以下に cg.c を掲げる。
1 /*
2 * cg.c --- CG 法のサンプル・プログラム
3 * コンパイル: gcc -O -o cg cg.c -lm
4 */
5
6 #include <stdio.h>
7 #include <math.h>
8 #include <stdlib.h>
9
10 typedef double *vector;
11 typedef vector *matrix;
12
13 vector new_vector(int n) { return malloc(sizeof(double) * n); }
14 void free_vector(vector v) { free(v); }
15
16 matrix new_matrix(int nrow, int ncol)
17 {
18 int i;
19 vector ap;
20 matrix a;
21
22 if ((a = malloc(nrow * sizeof(vector))) == NULL)
23 return NULL;
24 if ((ap = malloc(nrow * ncol * sizeof(double))) == NULL) {
25 free(a); return NULL;
26 }
27 for (i = 0; i < nrow; i++)
28 a[i] = ap + (i * ncol);
29 return a;
30 }
31
32 void free_matrix(matrix a)
33 {
34 free(a[0]);
35 free(a);
36 }
37
38 /* m行n列の零行列を作って返す */
39 matrix zeros(int m, int n)
40 {
41 int i, j;
42 matrix tmp;
43 tmp = new_matrix(m, n);
44 for (i = 0; i < m; i++)
45 for (j = 0; j < n; j++)
46 tmp[i][j] = 0.0;
47 return tmp;
48 }
49
50 /* n次元のベクトル z を零ベクトルにする */
51 void clear_vector(int n, vector z)
52 {
53 int i;
54 for (i = 0; i < n; i++)
55 z[i] = 0.0;
56 }
57
58 // n 次元ベクトルの内容を表示する
59 void print_vector(int n, vector x)
60 {
61 int i;
62 for (i = 0; i < n; i++) {
63 printf(" %g", x[i]);
64 if (i % 5 == 4)
65 printf("\n");
66 }
67 if (i % 5 != 0)
68 printf("\n");
69 }
70
71 // m行n列の行列の内容を表示する
72 void print_matrix(int m, int n, matrix a)
73 {
74 int i;
75 for (i = 0; i < m; i++)
76 print_vector(n, a[i]);
77 }
78
79 // ベクトル x, y の内積 (x,y) を計算する
80 double dotproduct(int n, vector x, vector y)
81 {
82 int i;
83 double sum;
84 sum = 0;
85 for (i = 0; i < n; i++)
86 sum += x[i] * y[i];
87 return sum;
88 }
89
90 // ベクトル x のノルム ||x|| を計算する
91 double norm(int n, vector x)
92 {
93 return sqrt(dotproduct(n, x, x));
94 }
95
96 // 行列 A, ベクトル b の積 A b を計算して、ベクトル Ab に書き込む
97 void multiply_mv(int m, int n, vector Ab, matrix A, vector b)
98 {
99 int i;
100 for (i = 0; i < m; i++)
101 Ab[i] = dotproduct(n, A[i], b);
102 }
103
104 // ベクトル x をベクトル y にコピーする
105 void copy_vector(int n, vector y, vector x)
106 {
107 int i;
108 for (i = 0; i < n; i++)
109 y[i] = x[i];
110 }
111
112 /* CG 法で n 元連立1次方程式 A x=b を解く。
113 * 呼び出すときの x は初期ベクトル, 戻ってきたときの x は解となる。
114 * ||残差|| < ε ||b|| となったら反復を停止する。
115 * 反復回数を maxiter に返す。
116 */
117
118 void cg(int n, matrix A, vector b, vector x, double eps, int *maxiter)
119 {
120 int i, k;
121 double eps_b, pAp, alpha, beta;
122 vector r, p, Ap;
123
124 r = new_vector(n);
125 p = new_vector(n);
126 Ap = new_vector(n);
127 /* eps_b := ε ||b|| */
128 eps_b = eps * norm(n, b);
129 /* r := A x */
130 multiply_mv(n, n, r, A, x);
131 /* r := b - A x */
132 for (i = 0; i < n; i++)
133 r[i] = b[i] - r[i];
134 /* p := r */
135 copy_vector(n, p, r);
136 /* 反復 */
137 for (k = 0; k < n; k++) {
138 /* Ap := A * p */
139 multiply_mv(n, n, Ap, A, p);
140 /* pAp := (p, Ap) */
141 pAp = dotproduct(n, p, Ap);
142 /* α = (r,p)/(p,Ap) */
143 alpha = dotproduct(n, r, p) / pAp;
144 /* x, r の更新 */
145 for (i = 0; i < n; i++) {
146 /* x = x + α*p */
147 x[i] += alpha * p[i];
148 /* r = r - α*A*p */
149 r[i] -= alpha * Ap[i];
150 }
151 /* ||r||<ε||b|| ならば反復を終了 */
152 if (norm(n, r) < eps_b)
153 break;
154 /* β := - (r, Ap) / (p, Ap) */
155 beta = - dotproduct(n, r, Ap) / pAp;
156 /* p := r + β*p */
157 for (i = 0; i < n; i++)
158 p[i] = r[i] + beta * p[i];
159 }
160 if (maxiter != NULL)
161 *maxiter = k;
162 }
163
164 int main()
165 {
166 int n, N, i, niter;
167 matrix a;
168 vector x, b;
169 double c1, c2;
170 /* 問題のサイズを決める */
171 printf("c1, c2, N=");
172 scanf("%lf %lf %d", &c1, &c2, &N);
173 n = N / 2;
174 N = 2 * n;
175 /* 行列 A, ベクトル x, b を記憶する変数の準備 */
176 a = zeros(N, N);
177 x = new_vector(N);
178 b = new_vector(N);
179 /* 係数行列を決める */
180 for (i = 0; i < N; i++) {
181 a[i][i] = 4.0;
182 if (i != 0)
183 a[i][i-1] = -1.0;
184 if (i != N-1)
185 a[i][i+1] = -1.0;
186 }
187 a[n-1][n] = 0.0;
188 a[n][n-1] = 0.0;
189 /* */
190 for (i = 0; i < n; i++) {
191 if (i != 0)
192 a[i][i-1] *= c1;
193 a[i][i] *= c1;
194 a[i][i+1] *= c1;
195 }
196 for (i = n; i < N; i++) {
197 a[i][i-1] *= c2;
198 a[i][i] *= c2;
199 if (i != N-1)
200 a[i][i+1] *= c2;
201 }
202 /* 係数行列を表示する */
203 if (N < 10) {
204 printf("a=\n");
205 print_matrix(N, N, a);
206 }
207 /* 解 x を決める */
208 for (i = 0; i < N; i++)
209 x[i] = i;
210 /* 解に対応する右辺 b を計算する */
211 multiply_mv(N, N, b, a, x);
212 /* b を表示する */
213 if (N < 10) {
214 printf("b=\n");
215 print_vector(N, b);
216 }
217 /* 初期ベクトル x を零ベクトルにする */
218 clear_vector(N, x);
219 /* CG 法で解く */
220 cg(N, a, b, x, 1e-12, &niter);
221 /* 解を表示する */
222 printf("x=\n");
223 print_vector(N, x);
224 printf("反復回数=%d\n", niter);
225
226 return 0;
227 }