解の対称性を保つCG法プログラム

        2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )

1. 概要

 差分法による離散化で発生する規則的対称疎行列を係数とする連立一次方程式Ax=bをCG法
で数値解xを計算するFORTRAN及びCのソースプログラムを掲載する。
本プログラムは、解析対象場が2次元版では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(k.grad(u))=f, k=1, f=100を2次元差分法で離散化して疎行列を作成している。
 行列は5行の疎行列で、対称行列のため下三角部分を持ち左からA1,A2,A3(対角)順である。

(1) 結果xが対称なプログラム

 (a) 計算方法 :丸め誤差まで含めて対称な計算が行えるよう、ループを分割して計算順序を
  揃える工夫が施してある。コンパイラーでループ結合しなければ対称な計算解xが保障される。
 (b) プログラム : FORTRANCのソースプログラム

(2) 結果xが対称なプログラム(最適化依存)

 (a) 計算方法 : 丸め誤差まで含めて対称な計算が行えるよう、計算順序を揃える工夫が施して
  ある。但し、ループ分割して明示的に計算順序を保障したものではなく、同一式中の計算順序
  が保障されない最適化では、計算解xの対称性も保障されない。
 (b) プログラム : FORTRANCのソースプログラム

(3) 通常のCG法プログラム

 (a) 計算方法 : 丸め誤差のため、計算解xの対称性は保たれない。
 (b) プログラム : FORTRANCのソースプログラム

(4) テストプログラム

 (a) 計算方法 : 計算結果は、収束状況とy軸の中心の各x-格子点での解xの値を出力する。
  結果のy軸対称性を調べるため、y軸対称とならない結果を全て出力できる。
  FORTRANは内部表現の数値を16進で出力し、Cは10進数に変換して出力している。
 (b) プログラム : FORTRANCのソースプログラム

3. 3次元CG法プログラム

 x,y,z方向の分割数NX,NY,NZ入力し、3次元差分法で離散化した、連立一次方程式Ax=bの解xを
 CG法(共役勾配法)で求める。
 -div(k.grad(u))=f, k=1, f=100を3次元差分法で離散化して疎行列を作成している。
 行列は7行の疎行列で、対称行列のため下三角部分を持ち左からA1,A2,A3,A4(対角)順である。

(1) 結果xが対称なプログラム

 (a) 計算方法 :丸め誤差まで含めて対称な計算が行えるよう、ループを分割して計算順序を
  揃える工夫が施してある。コンパイラーでループ結合しなければ対称な計算解xが保障される。
 (b) プログラム : FORTRANCのソースプログラム

(2) 結果xが対称なプログラム(最適化依存)

 (a) 計算方法 : 丸め誤差まで含めて対称な計算が行えるよう、計算順序を揃える工夫が施して
  ある。但し、ループ分割して明示的に計算順序を保障したものではなく、同一式中の計算順序
  が保障されない最適化では、計算解xの対称性も保障されない。
 (b) プログラム : FORTRANCのソースプログラム

(3) 通常のCG法プログラム

 (a) 計算方法 : 丸め誤差のため、計算解xの対称性は保たれない。
 (b) プログラム : FORTRANCのソースプログラム

(4) テストプログラム

 (a) 計算方法 : 計算結果は、収束状況とy,z軸の中心の各x-格子点での解xの値を出力する。
  結果のy,z軸対称性を調べるため、y,z軸対称とならない結果を全て出力できる。
  FORTRANは内部表現の数値を16進で出力し、Cは10進数に変換して出力している。
 (b) プログラム : FORTRANCのソースプログラム