Jacobi

/*

        jacobi

        imax 반복횟수

        b[i] 상수행렬

        a[i][j] 계수행렬

        residual 잔차 잔여

*/

#include <stdio.h>

#include <math.h>

#define nd 20

double x[nd], x0[nd], a[nd][nd], b[nd];

int imax = 1000, n, iter;

double epsilon = 0.00001, resavg;

void jacobi();

int main() {

        int i, j;

        printf("방정식 수");

        scanf("%d", &n);

        printf("계수행렬 a(i,j) 을 입력하시오");

        for (i = 1; i <= n; i = i + 1)

        {

                for (j = 1; j <= n; j = j + 1) {

                        printf("a(%d,j) ", i);

                        scanf("%lf", &a[i][j]);

                }

        }

        printf("상수행렬 b[i]를 입력하시오 \n");

        printf("b[i] i=1,,,n=");

        for (i = 1; i <= n; i = i + 1) scanf("%lf", &b[i]);

        printf("가정값 x0[i] 입력하시오 \n");

 printf("x0[i] i=1,,,n=");

        for (i = 1; i <= n; i = i + 1) scanf("%lf", &x0[i]);

        printf("\n INPUT DATA 출력");

        for (i = 1; i <= n; i = i + 1)

        {

                printf("\n");

                for (j = 1; j <= n; j = j + 1)

                {

                        printf("%10.5f", a[i][j]);

                }

                printf("%10.5f", b[i]);

        }

        printf("\n");

        jacobi();

        printf("\n 계산결과 출력");

        printf("\n iter = %d, 평균 residual = %f", iter, resavg);

        for (i = 1; i <= n; i++) {

                printf("\n x(%d)=%10.6f", i, x[i]);

        }

}

void jacobi() {

        double ahat[nd][nd], bhat[nd], res[nd];

        double sum, ressum;

        int i, j;

for (i = 1; i <= n; i++) {

                for (j = 1; j <= n; j++) {

                        ahat[i][j] = a[i][j] / a[i][i];

                }

                bhat[i] = b[i] / a[i][i];

                ahat[i][i] = 0;

        }

        for (iter = 1; iter <= imax; ++iter) {

                for (i = 1; i <= n; ++i) {

                        x[i] = bhat[i];

                        for (j = 1; j <= n; ++j) {

                                x[i] = x[i] - ahat[i][j] * x0[j];

                        }

                }

                ressum = 0;

                for (i = 1; i <= n; i++) {

                        sum = x[i] - bhat[i];

                        for (j = 1; j <= n; ++j) {

                                sum = sum + ahat[i][j] * x[j];

                        }

                        res[i] = fabs(sum);

                        ressum = ressum + res[i];

                }

                resavg = ressum / n;

                if (resavg < epsilon) {

                        goto END;

                }

                for (i = 1; i <= n; ++i) {

                        x0[i] = x[i];

                }

        }

        printf("\n 수렴하지 않는다");

END:;

}