야코비(Jacobi) Method C언어로 구현하기

1.  
2. /* Jacobi Method */
3.  
4. #include <stdio.h>
5. #include <conio.h>
6. #include <stdlib.h>
7. #include <math.h>
8.  
9. int it=0;
10. double **a;
11. double *b;
12. double *s;
13. #define         EPS             1e-5            // 최대 오차
14. #define         ITOR    500             // 최대 반복 횟수
15.  
16. int n; // 행렬의 갯수
17.  
18. FILE *fp;
19.  
20.  
21. double** matrix_malloc(int n)
22. {
23.  
24.         double **a;
25.  
26.         int i,j;
27.         a=(double **)malloc(n*sizeof(double));
28.         for(i=0;i<n;i++)
29.         a[i]=(double *)malloc(n*sizeof(double));
30.          for(i=0;i<n;i++){
31.                  for(j=0;j<n;j++){
32.                          a[i][j] = 0.0;
33.                  }
34.          }
35.  
36.         return a;
37.  
38. }
39. double * b_malloc(int n)
40. {
41.         double *b;
42.  
43.         b = (double*)malloc(n*sizeof(double));
44.  
45.         return b;
46.  
47. }
48.  
49. double * s_malloc(int n)
50. {
51.         double *s;
52.  
53.         s = (double*)malloc(n*sizeof(double));
54.  
55.         return s;
56.  
57. }
58.  
59.  
60. void matrix_free(int n)
61. {
62.  
63.         int i;
64.          
65.         for(i=0;i<n;i++)
66.                 free(a[i]);
67.  
68.          free(a);
69.  
70. }
71.  
72. void b_free()
73. {
74.         free(b);
75. }
76.  
77.  
78.  
79. void s_free()
80. {
81.         free(s);
82. }
83.  
84. void Jacobi_Method(void)
85. {
86.         int i,j;
87.     double a_norm, x_norm, w;
88.     double *x, *x1;
89.         int it =0;
90.  
91.         x = s_malloc(n);
92.         x1 = s_malloc(n);
93.  
94.     for(i=0;i<n;i++)
95.             x[i] = 0.0;
96.                
97.          printf("Iteration");
98.         for(i=0;i<n;i++)
99.                 printf("    x%d    ",i);
100.         printf("\n");
101.  
102.  
103.  
104.           do{
105.                 it++;
106.                 for(j=0;j<n;j++)
107.                                         x1[j]=x[j];
108.  
109.                 a_norm = 0.0;
110.                 x_norm = 0.0;
111.  
112.              
113.                                 for(i=0;i<n;i++) {
114.                                         w = b[i];
115.                                         for(j=0;j<n;j++)
116.                                         {
117.                                                 if(j!=i)w -= a[i][j]*x1[j];
118.                                         }
119.  
120.                         w /= a[i][i];
121.                         a_norm += fabs(x[i]-w);
122.                         x_norm += fabs(w);
123.                         x[i] = w;
124.                                 }
125.                 printf("    %-2d ",it);
126.                 for(i=0;i<n;i++)
127.                         printf(" %9.5lf", x[i]);
128.                 printf("\n");
129.                 if( it > ITOR){
130.                         printf("Iteration= %d    Check Error!!!!", it);
131.                         exit(1);
132.                 }
133.         } while(a_norm/x_norm > EPS);
134.  
135.          free(x);
136.          free(x1);
137.  
138. }
139.  
140.  
141. int main(void)
142. {
143.         int i, j;
144.  
145.  
146.        
147.         if((fp=fopen("matrix.txt", "r"))==NULL){
148.                  printf("File error!!\n");
149.                 exit(1);
150.         }
151.  
152.         fscanf(fp, "%d", &n);
153.  
154.         a = matrix_malloc(n);
155.         b = b_malloc(n);
156.         s = s_malloc(n);
157.  
158. // 파일 입력
159.         for(i=0;i<n;i++){
160.                  for(j=0;j<n;j++){
161.  
162.                          fscanf(fp,"%lf",&a[i][j]);
163.                          printf("%.2f ",a[i][j]);
164.                  }
165.                 fscanf(fp,"%lf",&b[i]);
166.         printf("%.2f \n",b[i]);
167.         }
168.  
169.  
170.         Jacobi_Method();
171.  
172.  
173.         matrix_free(n);
174.  
175.         b_free();
176.         s_free();
177.         fclose(fp);
178.        
179.         return 0;
180. }