解の対称性を保つCG法プログラム
2007/06/04 後 保範 (Ushiro Yasunori、東京工芸大学)
1. 概要
差分法による離散化で発生する規則的対称疎行列を係数とする連立一次方程式Ax=bをCG法
で数値解xを計算するFORTRAN及びCのソースプログラムを掲載する。
本プログラムは、解析対象場が2次元版ではxとyの中心で対称のとき、3次元版ではｘ,yとzの中心
で対称のとき、計算解xをそれぞれのもとで対称にするものである。比較のために丸め誤差の
ため、計算解xが対称にならない通常のプログラムも掲載している。
解析対象場がその軸で対称とは、FDMで離散化後の係数行列A及び右辺ベクトルがその軸で対称
になっていることである。
計算解xが対称になっているかどうかを確かめるため、結果xを2次元の場合はy軸対称性のテスト
結果を、3次元の場合はy,z両軸で対称性のテスト結果を出力する。
対称性のテスト結果は、非対称となる格子点の組数及びその詳細を出力する。
FORTRANは16進数で内部表現の数値をそのまま出力している。一方、Cは10進数に変換した結果
を出力している。
本プログラムは内容を理解してもらうことを目的とし、計算時間及び使用メモリ量の効率化
は目的としていない。配列の格納順はCの場合も＃define文を使用しFORTRANに合わせてある。
2. 2次元CG法プログラム
x及びy方向の分割数NX,NY入力し、2次元差分法で離散化した、連立一次方程式Ax=bの解xを
CG法(共役勾配法)で求める。
-div(ｋ.grad(u))=f, k=1, f=100を2次元差分法で離散化して疎行列を作成している。
行列は5行の疎行列で、対称行列のため下三角部分を持ち左からA1,A2,A3(対角)順である。
(1) 結果xが対称なプログラム
(a) 計算方法 ：丸め誤差まで含めて対称な計算が行えるよう、q=Apの計算においてペアーとなる
項を括弧で括る工夫がしてある。コンパイラーで括弧を外さなければ対称な計算解xが保障される。
(b) プログラム ：
FORTRAN(CG2DS.f),と
C(CG2DS.c)のソースプログラム
(2) 通常のCG法プログラム
(a) 計算方法 ： 丸め誤差のため、計算解xの対称性は保たれない。
(b) プログラム ：
FORTRAN(CG2D2.f),と
C(CG2D2.c)のソースプログラム
(3) テストプログラム
(a) 計算方法 ： 計算結果は、収束状況とy軸の中心の各x-格子点での解xの値を出力する。
結果のy軸対称性を調べるため、y軸対称とならない結果を全て出力できる。
FORTRANは内部表現の数値を16進で出力し、Cは10進数に変換して出力している。
(b) プログラム ：
FORTRAN(TCGD2S.f),と
C(TCGD2S.c)のソースプログラム
3. 3次元CG法プログラム
x,y,z方向の分割数NX,NY,NZ入力し、3次元差分法で離散化した、連立一次方程式Ax=bの解xを
CG法(共役勾配法)で求める。
-div(ｋ.grad(u))=f, k=1, f=100を3次元差分法で離散化して疎行列を作成している。
行列は7行の疎行列で、対称行列のため下三角部分を持ち左からA1,A2,A3,A4(対角)順である。
(1) 結果xが対称なプログラム
(a) 計算方法 ：丸め誤差まで含めて対称な計算が行えるよう、q=Apの計算においてペアーとなる
項を括弧で括る工夫がしてある。コンパイラーで括弧を外さなければ対称な計算解xが保障される。
(b) プログラム ：
FORTRAN(CG3DS.f),と
C(CG3DS.c)のソースプログラム
(2) 通常のCG法プログラム
(a) 計算方法 ： 丸め誤差のため、計算解xの対称性は保たれない。
(b) プログラム ：
FORTRAN(CG3D2.f),と
C(CG3D2.c)のソースプログラム
(3) テストプログラム
(a) 計算方法 ： 計算結果は、収束状況とy,z軸の中心の各x-格子点での解xの値を出力する。
結果のy,z軸対称性を調べるため、y,z軸対称とならない結果を全て出力できる。
FORTRANは内部表現の数値を16進で出力し、Cは10進数に変換して出力している。
(b) プログラム ：
FORTRAN(TCGD3S.f), と
C(TCGD3S.c)のソースプログラム