2次元FDM用CG法(C)プログラム Y-対称NO.1
2003/08/26 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求める。
α=(r,r)/(p,Ap),β=new(r,r)/old(r,r)を使用した、演算量の少ないバージョン
係数行列Aと右辺ベクトルbが解析領域に対してyの中心で対称なら、計算解xも対称に
なる工夫をしたプログラムである。但し、Cコンパイラーが同一式中の加減算の順序を
変更するときには、対称性は保証されない。
2. プログラム
#include <stdio.h>
#include <math.h>
// Global Define
#define NDX 301
#define NDY 301
#define A(i,j,l) A[l-1][j][i]
#define B(i,j) B[j][i]
#define X(i,j) X[j][i]
#define R(i,j) R[j][i]
#define P(i,j) P[j][i]
#define Q(i,j) Q[j][i]
extern double A[3][NDY][NDX], B[NDY][NDX], X[NDY][NDX] ;
extern double R[NDY][NDX], P[NDY][NDX], Q[NDY][NDX] ;
extern FILE *FT1 ;
//=================================================================C
int CG2DS(int NX, int NY, double EPS, int *ITER, double *ERR)
//=================================================================C
// Solve Ax=b by CG NO.2 with 2 dimensional FDM C
// Alpha=(R,R)/(P,AP), Beta=new(R,R)/old(R,R) C
// for symmetric solution on Y-axis C
//-----------------------------------------------------------------C
// NX I*4, In, Grid Numbers on X-axis C
// NY I*4, In, Grid Numbers on Y-axis C
// EPS R*8, In, if ||r||/||b|| <= EPS --> return C
// ITER I*4, I/O, Number of Iteration C
// ERR R*8, Out, ERR=||r||/||b|| C
// return I*4, Out, IERR=0, Normal Return C
// =1, No Convergent C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/08/25 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j, k, kk, MNY, NY2, IERR ;
double BN, C0, C1, Alpha, Beta ;
// Even or Odd of NY
MNY = NY % 2 ;
if(MNY == 0) { NY2 = NY/2 - 1 ; }
else { NY2 = (NY - 1)/2 ; }
printf("MNY=%d NY2=%d\n",MNY,NY2) ;
// P=R=B-A*X
IERR = 0 ;
BN = 0.0 ;
C0 = 0.0 ;
for (j=1; j<=NY2; j++) {
k = NY - j ;
for (i=1; i<=NX-1; i++) {
BN = BN + B(i,j)*B(i,j) + B(i,k)*B(i,k) ;
R(i,j) = B(i,j) - A(i,j,1)*X(i,j-1) - A(i,j,2)*X(i-1,j)
- A(i,j,3)*X(i,j) - A(i+1,j,2)*X(i+1,j)
- A(i,j+1,1)*X(i,j+1) ;
R(i,k) = B(i,k) - A(i,k+1,1)*X(i,k+1) - A(i,k,2)*X(i-1,k)
- A(i,k,3)*X(i,k) - A(i+1,k,2)*X(i+1,k)
- A(i,k,1)*X(i,k-1) ;
C0 = C0 + R(i,j)*R(i,j) + R(i,k)*R(i,k) ;
P(i,j) = R(i,j) ;
P(i,k) = R(i,k) ; }
}
// Y-center line
if(MNY == 0) {
j = NY/2 ;
for (i=1; i<=NX-1; i++) {
BN = BN + B(i,j)*B(i,j) ;
R(i,j) = B(i,j) - A(i,j,1)*X(i,j-1) - A(i,j,2)*X(i-1,j)
- A(i,j,3)*X(i,j) - A(i+1,j,2)*X(i+1,j)
- A(i,j+1,1)*X(i,j+1) ;
C0 = C0 + R(i,j)*R(i,j) ;
P(i,j) = R(i,j) ; }
}
// Main Loop
for (kk=1; kk<=*ITER; kk++) {
// Q=A*P, C=(P,Q), Alpha=C0/(P,Q)
Alpha = 0.0 ;
for (j=1; j<=NY2; j++) {
k = NY - j ;
for (i=1; i<=NX-1; i++) {
Q(i,j) = A(i,j,1)*P(i,j-1) + A(i,j,2)*P(i-1,j)
+ A(i,j,3)*P(i,j) + A(i+1,j,2)*P(i+1,j)
+ A(i,j+1,1)*P(i,j+1) ;
Q(i,k) = A(i,k+1,1)*P(i,k+1) + A(i,k,2)*P(i-1,k)
+ A(i,k,3)*P(i,k) + A(i+1,k,2)*P(i+1,k)
+ A(i,k,1)*P(i,k-1) ;
Alpha = Alpha + P(i,j)*Q(i,j) + P(i,k)*Q(i,k) ; }
}
// Y-center line
if(MNY == 0) {
j = NY/2 ;
for (i=1; i<=NX-1; i++) {
Q(i,j) = A(i,j,1)*P(i,j-1) + A(i,j,2)*P(i-1,j)
+ A(i,j,3)*P(i,j) + A(i+1,j,2)*P(i+1,j)
+ A(i,j+1,1)*P(i,j+1) ;
Alpha = Alpha + P(i,j)*Q(i,j) ; }
}
Alpha = C0/Alpha ;
// X=X+Alpha*P, R=R-Alpha*Q
C1 = 0.0 ;
for (j=1; j<=NY-1; j++) {
for (i=1; i<=NX-1; i++) {
X(i,j) = X(i,j) + Alpha*P(i,j) ;
R(i,j) = R(i,j) - Alpha*Q(i,j) ;
C1 = C1 + R(i,j)*R(i,j) ; }
}
// if(ERR <= EPS) Convergent, Beta=C1/C0
*ERR = sqrt(C1/BN) ;
if(*ERR <= EPS) goto M100 ;
Beta = C1/C0 ;
C0 = C1 ;
// P=R+Beta*P
for (j=1; j<=NY-1; j++) {
for (i=1; i<=NX-1; i++) {
P(i,j) = R(i,j) + Beta*P(i,j) ; }
}
}
IERR = 1 ;
M100: *ITER = kk ;
return (IERR) ; }