3次元FDM用CG法(C)プログラム Y-Z-対称NO.2

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

1. 概要

 3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求める。
α=(r,r)/(p,Ap),β=new(r,r)/old(r,r)を使用した、演算量の少ないバージョン
係数行列Aと右辺ベクトルbが解析領域に対してy及びzの中心で対称なら、計算解xも
対称になる工夫をしたプログラムである。Cコンパイラーが同一式中の加減算の順序を
変更する場合も、対称性を保証する。但し、ループ結合まで行う最適化の場合は結果の
対称性は保証されない。

2. プログラム

#include <stdio.h>
#include <math.h>
// Global Define
#define NDX  51
#define NDY  51
#define NDZ  51
#define A(i,j,k,l)  A[l-1][k][j][i]
#define B(i,j,k)    B[k][j][i]
#define X(i,j,k)    X[k][j][i]
#define R(i,j,k)    R[k][j][i]
#define P(i,j,k)    P[k][j][i]
#define Q(i,j,k)    Q[k][j][i]
extern double A[4][NDZ][NDY][NDX], B[NDZ][NDY][NDX], X[NDZ][NDY][NDX] ;
extern double R[NDZ][NDY][NDX], P[NDZ][NDY][NDX], Q[NDZ][NDY][NDX] ;
extern FILE *FT1 ;
//=================================================================C
int CG3DM(int NX, int NY, int NZ, double EPS, int *ITER, double *ERR)
//=================================================================C
//  Solve Ax=b by CG NO.1 with 3 dimensional FDM                   C
//    Alpha=(P,R)/(P,AP),   Beta=-(R,AP)/(P,AP)                    C
//    for symmetric solution on Y-axis                             C
//    Modified version of computing Q=A*P for symmetric solution   C
//-----------------------------------------------------------------C
//    NX           I*4, In,  Grid Numbers on X-axis                C
//    NY           I*4, In,  Grid Numbers on Y-axis                C
//    NZ           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, l, m, kk, NY2, NZ2, MNY, MNZ, IERR ;
  double  BN, C0, C1, Alpha, Beta ;
//  Even or Odd of NY,NZ
  MNY = NY % 2 ;
  if(MNY == 0) { NY2 = NY/2 - 1 ; }
  else { NY2 = (NY - 1)/2 ; }
  MNZ = NZ % 2 ;
  if(MNZ == 0) { NZ2 = NZ/2 - 1 ; }
  else { NZ2 = (NZ - 1)/2 ; }
//  P=R=B-A*X
  IERR = 0 ;
  BN   = 0.0 ;
  C0   = 0.0 ;
//   main points
  for (k=1; k<=NZ2; k++) {
    m = NZ - k ; 
    for (j=1; j<=NY2; j++) {
      l = NY - j ; 
//    Divided do loop for symmetric solution
      for (i=1; i<=NX-1; i++) {
        BN  = BN + B(i,j,k)*B(i,j,k) + B(i,j,m)*B(i,j,m)
                 + B(i,l,k)*B(i,l,k) + B(i,l,m)*B(i,l,m) ; 
        R(i,j,k) = B(i,j,k) - A(i,j,k,1)*X(i,j,k-1) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m+1,1)*X(i,j,m+1) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k,1)*X(i,l,k-1) ;
        R(i,l,m) = B(i,l,m) - A(i,l,m+1,1)*X(i,l,m+1) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,2)*X(i,j-1,k) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m,2)*X(i,j-1,m) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k+1,2)*X(i,l+1,k) ;
        R(i,l,m) = B(i,l,m) - A(i,l+1,m,2)*X(i,l+1,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,3)*X(i-1,j,k) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m,3)*X(i-1,j,m) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k,3)*X(i-1,l,k) ;
        R(i,l,m) = B(i,l,m) - A(i,l,m,3)*X(i-1,l,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,4)*X(i,j,k) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m,4)*X(i,j,m) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k,4)*X(i,l,k) ;
        R(i,l,m) = B(i,l,m) - A(i,l,m,4)*X(i,l,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i+1,j,k,3)*X(i+1,j,k) ;
        R(i,j,m) = B(i,j,m) - A(i+1,j,m,3)*X(i+1,j,m) ;
        R(i,l,k) = B(i,l,k) - A(i+1,l,k,3)*X(i+1,l,k) ;
        R(i,l,m) = B(i,l,m) - A(i+1,l,m,3)*X(i+1,l,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j+1,k,2)*X(i,j+1,k) ;
        R(i,j,m) = B(i,j,m) - A(i,j+1,m,2)*X(i,j+1,m) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k,2)*X(i,l-1,k) ;
        R(i,l,m) = B(i,l,m) - A(i,l,m,2)*X(i,l-1,m) ; }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k+1,1)*X(i,j,k+1) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m,1)*X(i,j,m-1) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k+1,1)*X(i,l,k+1) ;
        R(i,l,m) = B(i,l,m) - A(i,l,m,1)*X(i,l,m-1) ;
        C0 = C0 + R(i,j,k)*R(i,j,k) + R(i,j,m)*R(i,j,m) 
                + R(i,l,k)*R(i,l,k) + R(i,l,m)*R(i,l,m) ;
        P(i,j,k) = R(i,j,k) ;
        P(i,j,m) = R(i,j,m) ;
        P(i,l,k) = R(i,l,k) ;
        P(i,l,m) = R(i,l,m) ;  }
    } }
//   Y center line
  if(MNY == 0) {
    for (k=1; k<=NZ2; k++) {
      m = NZ - k ;
      j = NY/2 ;
//    Divided do loop for symmetric solution
      for (i=1; i<=NX-1; i++) {
        BN = BN + B(i,j,k)*B(i,j,k) + B(i,j,m)*B(i,j,m) ;
        R(i,j,k) = B(i,j,k) - A(i,j,k,1)*X(i,j,k-1) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m+1,1)*X(i,j,m+1) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,2)*X(i,j-1,k) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m,2)*X(i,j-1,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,3)*X(i-1,j,k) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m,3)*X(i-1,j,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,4)*X(i,j,k) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m,4)*X(i,j,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i+1,j,k,3)*X(i+1,j,k) ;
        R(i,j,m) = B(i,j,m) - A(i+1,j,m,3)*X(i+1,j,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j+1,k,2)*X(i,j+1,k) ;
        R(i,j,m) = B(i,j,m) - A(i,j+1,m,2)*X(i,j+1,m) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k+1,1)*X(i,j,k+1) ;
        R(i,j,m) = B(i,j,m) - A(i,j,m,1)*X(i,j,m-1) ;
        C0  = C0 + R(i,j,k)*R(i,j,k) + R(i,j,m)*R(i,j,m) ;
        P(i,j,k) = R(i,j,k) ;
        P(i,j,m) = R(i,j,m) ;  }
    }  }
//   Z center line
  if(MNZ == 0 ) {
    k = NZ/2 ;
    for (j=1; j<=NY2; j++) {
	  l = NY - j ;
//    Divided do loop for symmetric solution
      for (i=1; i<=NX-1; i++) {
        BN = BN + B(i,j,k)*B(i,j,k) + B(i,l,k)*B(i,l,k) ;
        R(i,j,k) = B(i,j,k) - A(i,j,k,1)*X(i,j,k-1) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k,1)*X(i,l,k-1) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,2)*X(i,j-1,k) ;
        R(i,l,k) = B(i,l,k) - A(i,l+1,k,2)*X(i,l+1,k) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,3)*X(i-1,j,k) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k,3)*X(i-1,l,k) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k,4)*X(i,j,k) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k,4)*X(i,l,k) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i+1,j,k,3)*X(i+1,j,k) ;
        R(i,l,k) = B(i,l,k) - A(i+1,l,k,3)*X(i+1,l,k) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j+1,k,2)*X(i,j+1,k) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k,2)*X(i,l-1,k) ;  }
      for (i=1; i<=NX-1; i++) {
        R(i,j,k) = B(i,j,k) - A(i,j,k+1,1)*X(i,j,k+1) ;
        R(i,l,k) = B(i,l,k) - A(i,l,k+1,1)*X(i,l,k+1) ;
        C0 = C0 + R(i,j,k)*R(i,j,k) + R(i,l,k)*R(i,l,k) ;
        P(i,j,k) = R(i,j,k) ;
        P(i,l,k) = R(i,l,k) ;  }
    }  }
 //   Y and Z center line
  if(MNY == 0 & MNZ ==0) {
    k = NZ/2 ;
    j = NY/2 ;
    for (i=1; i<=NX-1; i++) {
      BN       = BN + B(i,j,k)*B(i,j,k) ;
      R(i,j,k) = B(i,j,k) - A(i,j,k,1)*X(i,j,k-1)
               - A(i,j,k,2)*X(i,j-1,k) - A(i,j,k,3)*X(i-1,j,k)
               - A(i,j,k,4)*X(i,j,k) - A(i+1,j,k,3)*X(i+1,j,k)
               - A(i,j+1,k,2)*X(i,j+1,k) - A(i,j,k+1,1)*X(i,j,k+1) ;
      C0       = C0 + R(i,j,k)*R(i,j,k) ;
      P(i,j,k) = R(i,j,k) ;  }
    }
//  Main Loop
  for (kk=1; kk<=*ITER; kk++) {
//   Q=A*P, C=(P,Q), Alpha=C0/(P,Q)
    Alpha = 0.0 ;
    for (k=1; k<=NZ2; k++) {
      m = NZ - k ;
      for (j=1; j<=NY2; j++) {
        l = NY - j ;
//    Divided do loop for symmetric solution
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = A(i,j,k,1)*P(i,j,k-1) + A(i,j,k,2)*P(i,j-1,k) ;
          Q(i,j,m) = A(i,j,m+1,1)*P(i,j,m+1) + A(i,j,k,2)*P(i,j-1,m) ;
          Q(i,l,k) = A(i,l,k,1)*P(i,l,k-1) + A(i,l+1,k,2)*P(i,l+1,k) ;
          Q(i,l,m) = A(i,l,m+1,1)*P(i,l,m+1) + A(i,l+1,m,2)*P(i,l+1,m) ; }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k,3)*P(i-1,j,k) ;
          Q(i,j,m) = Q(i,j,m) + A(i,j,m,3)*P(i-1,j,m) ;
          Q(i,l,k) = Q(i,l,k) + A(i,l,k,3)*P(i-1,l,k) ;
          Q(i,l,m) = Q(i,l,m) + A(i,l,m,3)*P(i-1,l,m) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k,4)*P(i,j,k) ;
          Q(i,j,m) = Q(i,j,m) + A(i,j,k,4)*P(i,j,m) ;
          Q(i,l,k) = Q(i,l,k) + A(i,l,k,4)*P(i,l,k) ;
          Q(i,l,m) = Q(i,l,m) + A(i,l,m,4)*P(i,l,m) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i+1,j,k,3)*P(i+1,j,k) ;
          Q(i,j,m) = Q(i,j,m) + A(i+1,j,m,3)*P(i+1,j,m) ;
          Q(i,l,k) = Q(i,l,k) + A(i+1,l,k,3)*P(i+1,l,k) ;
          Q(i,l,m) = Q(i,l,m) + A(i+1,l,m,3)*P(i+1,l,m) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j+1,k,2)*P(i,j+1,k) ;
          Q(i,j,m) = Q(i,j,m) + A(i,j+1,k,2)*P(i,j+1,m) ;
          Q(i,l,k) = Q(i,l,k) + A(i,l,k,2)*P(i,l-1,k) ;
          Q(i,l,m) = Q(i,l,m) + A(i,l,m,2)*P(i,l-1,m) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k+1,1)*P(i,j,k+1) ;
          Q(i,j,m) = Q(i,j,m) + A(i,j,m,1)*P(i,j,m-1) ;
          Q(i,l,k) = Q(i,l,k) + A(i,l,k+1,1)*P(i,l,k+1) ;
          Q(i,l,m) = Q(i,l,m) + A(i,l,m,1)*P(i,l,m-1) ;
          Alpha = Alpha + P(i,j,k)*Q(i,j,k) + P(i,j,m)*Q(i,j,m)
                        + P(i,l,k)*Q(i,l,k) + P(i,l,m)*Q(i,l,m) ; }
     } }
//    Y center line
    if(MNY == 0) {
      for (k=1; k<=NZ2; k++) {
        m = NZ - k ;
        j = NY/2 ;
//    Divided do loop for symmetric solution
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = A(i,j,k,1)*P(i,j,k-1) + A(i,j,k,2)*P(i,j-1,k) ;
          Q(i,j,m) = A(i,j,m+1,1)*P(i,j,m+1) + A(i,j,k,2)*P(i,j-1,m) ; }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k,3)*P(i-1,j,k) ;
          Q(i,j,m) = Q(i,j,m) + A(i,j,m,3)*P(i-1,j,m) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k,4)*P(i,j,k) ;
          Q(i,j,m) = Q(i,j,m) + A(i,j,k,4)*P(i,j,m) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i+1,j,k,3)*P(i+1,j,k) ;
          Q(i,j,m) = Q(i,j,m) + A(i+1,j,m,3)*P(i+1,j,m) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j+1,k,2)*P(i,j+1,k) ;
          Q(i,j,m) = Q(i,j,m) + A(i,j+1,k,2)*P(i,j+1,m) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k+1,1)*P(i,j,k+1) ;
          Q(i,j,m) = Q(i,j,m) + A(i,j,m,1)*P(i,j,m-1) ;
          Alpha = Alpha + P(i,j,k)*Q(i,j,k) + P(i,j,m)*Q(i,j,m) ; }
      }  }
//    Z center line
    if(MNZ == 0) {
	  k = NZ/2 ;
      for (j=1; j<=NY2; j++) {
        l = NY - j ;
//    Divided do loop for symmetric solution
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = A(i,j,k,1)*P(i,j,k-1) + A(i,j,k,2)*P(i,j-1,k) ;
          Q(i,l,k) = A(i,l,k,1)*P(i,l,k-1) + A(i,l+1,k,2)*P(i,l+1,k) ; }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k,3)*P(i-1,j,k) ;
          Q(i,l,k) = Q(i,l,k) + A(i,l,k,3)*P(i-1,l,k) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k,4)*P(i,j,k) ;
          Q(i,l,k) = Q(i,l,k) + A(i,l,k,4)*P(i,l,k) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i+1,j,k,3)*P(i+1,j,k) ;
          Q(i,l,k) = Q(i,l,k) + A(i+1,l,k,3)*P(i+1,l,k) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j+1,k,2)*P(i,j+1,k) ;
          Q(i,l,k) = Q(i,l,k) + A(i,l,k,2)*P(i,l-1,k) ;  }
        for (i=1; i<=NX-1; i++) {
          Q(i,j,k) = Q(i,j,k) + A(i,j,k+1,1)*P(i,j,k+1) ;
          Q(i,l,k) = Q(i,l,k) + A(i,l,k+1,1)*P(i,l,k+1) ;
          Alpha = Alpha + P(i,j,k)*Q(i,j,k) + P(i,l,k)*Q(i,l,k) ;  }
        for (i=1; i<=NX-1; i++) {
       }
        for (i=1; i<=NX-1; i++) {
       }
	  }  }
//    Y-Z center line
    if(MNY == 0 & MNZ == 0) {
      k = NZ/2 ;
      j = NY/2 ;
      for (i=1; i<=NX-1; i++) {
        Q(i,j,k) = A(i,j,k,1)*P(i,j,k-1) + A(i,j,k,2)*P(i,j-1,k)
                 + A(i,j,k,3)*P(i-1,j,k) + A(i,j,k,4)*P(i,j,k)
                 + A(i+1,j,k,3)*P(i+1,j,k) + A(i,j+1,k,2)*P(i,j+1,k)
                 + A(i,j,k+1,1)*P(i,j,k+1) ;
        Alpha    = Alpha + P(i,j,k)*Q(i,j,k) ;  }
	  }
    Alpha = C0/Alpha ;
//   X=X+Alpha*P,  R=R-Alpha*Q
    C1   = 0.0 ;
    for (k=1; k<=NZ-1; k++) {
     for (j=1; j<=NY-1; j++) {
      for (i=1; i<=NX-1; i++) {
        X(i,j,k) = X(i,j,k) + Alpha*P(i,j,k) ;
        R(i,j,k) = R(i,j,k) - Alpha*Q(i,j,k) ;
        C1       = C1 + R(i,j,k)*R(i,j,k) ;  }
     } }
//   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 (k=1; k<=NZ-1; k++) {
     for (j=1; j<=NY-1; j++) {
      for (i=1; i<=NX-1; i++) {
        P(i,j,k) = R(i,j,k) + Beta*P(i,j,k) ;  }
     } }
   }
    IERR = 1 ;
//
 M100: *ITER = kk ;
 return (IERR) ; }