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