3次元FDM用CG法(C)テストプログラム

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

1. 概要

 3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求めるプログラムのテスト用プログラム
CG法は比較的直交性に強いバージョンと演算量の少ないバージョンの両方をテストする。

2. プログラム

//=================================================================C
//  Test Program of CG for 3-Dimensional FDM  with symmetric       C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2003/08/21                     C
//        ( Hitachi Ltd. and Waseda University )                   C
//=================================================================C
#include <stdio.h>
#include <math.h>
#include <time.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]
double A[4][NDZ][NDY][NDX], B[NDZ][NDY][NDX], X[NDZ][NDY][NDX] ;
double R[NDZ][NDY][NDX], P[NDZ][NDY][NDX], Q[NDZ][NDY][NDX] ;
FILE   *FT1 ;
// Function Definition
void SET3DS(int NX, int NY, int NZ) ;
void OUT3D(int NX, int NY, int NZ) ;
int CG3D1(int NX, int NY, int NZ, double EPS, int *ITER, double *ERR) ;
int CG3D2(int NX, int NY, int NZ, double EPS, int *ITER, double *ERR) ;
// Main Program
void main ()
{ int     k, NX, NY, NZ, ITER, IERR, T1, T2 ;
  double  EPS, ERR, CPU, F1, F2, FLOP, FLOPS ;
//  Open File
  FT1 = fopen("C-CG3D.txt","w") ;
  EPS = 1.0e-8 ;
//  Size Input
  printf("Type In  NX,NY,NZ (<=300) \n") ;
  scanf("%d %d %d",&NX,&NY,&NZ) ;
  if(NX > NDX-1) NX = NDX-1 ;
  if(NY > NDY-1) NY = NDY-1 ;
  if(NZ > NDZ-1) NZ = NDZ-1 ;
  fprintf(FT1,"  CG Method for 3-Dimensional FDM,  NX=%d, NY=%d, NZ=%d \n",
          NX,NY,NZ) ;
  fprintf(FT1," Type  Loop   Error   Time(s)   MFLOPS \n") ;
  for (k=1; k<=2; k++) {
//  Set A,B,X
    SET3DS(NX, NY, NZ) ;
//  Solve Ax=b by CG
    ITER = NX*NY*NZ ;
    T1 = clock() ;
    if(k == 1) {
       IERR = CG3D1(NX, NY, NZ, EPS, &ITER, &ERR) ;
       F1 = 14 ;
       F2 = 27 ;  }
     else {
       IERR = CG3D2(NX, NY, NZ, EPS, &ITER, &ERR) ;
       F1 = 16 ;
       F2 = 23 ; }
    T2 = clock() ;
//  Write
    CPU  = (double)(T2 - T1)/CLOCKS_PER_SEC ;
    FLOP = (NX-1)*(NY-1)*(NZ-1)*(F1+F2*ITER)*1.0e-6 ;
    FLOPS = 0.0 ;
    if(CPU > 0.0) FLOPS = FLOP/CPU ;
    fprintf(FT1,"%4d %4d %10.2e %8.1f %8.1f \n",
                  k,ITER,ERR,CPU,FLOPS) ;
    printf("%4d %4d %10.2e %8.1f %8.1f \n",k,ITER,ERR,CPU,FLOPS) ;
//  Output center-line
    OUT3D(NX, NY, NZ) ;
   }
 }
//=================================================================C
void SET3DS(int NX, int NY, int NZ)
//=================================================================C
//  Set A and B by 3 Dimensional FDM  with symmetric               C
//    -div(K.grad(X)) = f,   K=1.0, f=10,0 and 1x1 square          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
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2003/08/21                     C
//        ( Hitachi Ltd. and Waseda University )                   C
//=================================================================C
{ int    i, j, k ;
  double DX, DY, DZ, HX, HY, HZ, F ;
//  Clear A,B,X
  for (k=0; k<=NZ; k++) {
   for (j=0; j<=NY; j++) {
    for (i=0; i<=NX; i++) {
      A(i,j,k,1) = 0.0 ;
      A(i,j,k,2) = 0.0 ;
      A(i,j,k,3) = 0.0 ;
      A(i,j,k,4) = 0.0 ;
      B(i,j,k)   = 0.0 ;
      X(i,j,k)   = 0.0 ; }
   } }
//  Initial Data
  DX = NX ;
  DY = NY ;
  DZ = NZ ;
  HX = 1.0/DX ;
  HY = 1.0/DY ;
  HZ = 1.0/DZ ;
  F  = 100.0 ;
//  Set A,B
  for (k=1; k<=NZ-1; k++) {
   for (j=1; j<=NY-1; j++) {
    for (i=1; i<=NX-1; i++) {
      A(i,j,k,1) = -HX*HY*DZ ;
      if(k == 1) A(i,j,k,1) = 0.0 ;
      A(i,j,k,2) = -HZ*HX*DY ;
      if(j == 1) A(i,j,k,2) = 0.0 ;
      A(i,j,k,3) = -HY*HZ*DX ;
      if(i == 1) A(i,j,k,3) = 0.0 ;
      A(i,j,k,4) = 2.0*(HY*HZ*DX+HZ*HX*DY+HX*HY*DZ) ;
      B(i,j,k)   = HX*HY*HZ*F ;  }
   } }
 }
//=================================================================C
void OUT3D(int NX, int NY, int NZ)
//-----------------------------------------------------------------C
{ int    i ;
  double DX, XV ;
//  Output Y
  fprintf(FT1,"Output on YZ-center line \n","  X-axis  Solution(X) \n") ;
  DX = NX ;
  for (i=0; i<=NX; i++) {
    XV = i/DX ;
    fprintf(FT1,"%10.6f %10.6f \n",XV,X(i,NY/2,NZ/2)) ; }
  }