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

1. 概要

2. プログラム

//=================================================================C
//  Test Program of CG for 2-Dimensional FDM  with symmetric       C
//    for symmetric solution on Y-axis                             C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2003/08/25                     C
//        ( Hitachi Ltd. and Waseda University )                   C
//=================================================================C
#include <stdio.h>
#include <math.h>
#include <time.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]
double A[3][NDY][NDX], B[NDY][NDX], X[NDY][NDX] ;
double R[NDY][NDX], P[NDY][NDX], Q[NDY][NDX] ;
FILE   *FT1 ;
// Function Definition
void SET2DS(int NX, int NY) ;
void OUT2D(int NX, int NY, int ID) ;
int CG2DM(int NX, int NY, double EPS, int *ITER, double *ERR) ;
int CG2DS(int NX, int NY, double EPS, int *ITER, double *ERR) ;
int CG2D2(int NX, int NY, double EPS, int *ITER, double *ERR) ;
// Main Program
void main ()
{ int     k, NX, NY, ITER, IERR, T1, T2 ;
  double  EPS, ERR, CPU, F1, F2, FLOP, FLOPS ; 
//  Open File
  FT1 = fopen("C-CG2DS.txt","w") ;
  EPS = 1.0e-8 ;
//  Size Input
  printf("Type In  NX,NY (<=300) \n") ; 
  scanf("%d %d",&NX,&NY) ;
  if(NX > NDX-1) NX = NDX-1 ;
  if(NY > NDY-1) NY = NDY-1 ;
  for (k=1; k<=3; k++) {
//  Set A,B,X
    printf("----- 2D FDM with CG -------------------------- \n") ;
    fprintf(FT1,"----- 2D FDM with CG -------------------------- \n") ;
    if(k == 1) {
      fprintf(FT1,"Modified symmetric solution CG for Y-axis \n") ; }
    else { if(k == 2) {
      fprintf(FT1,"Symmetric solution CG for Y-axis \n") ; }
     else {
      fprintf(FT1,"Unsymmetric solution CG for Y-axis \n") ; }
     }
    fprintf(FT1," CG Method for 2-Dimensional FDM,  NX=%d, NY=%d \n",NX,NY) ;
    fprintf(FT1," Type  Loop   Error    Time(s)   MFLOPS \n") ;
	SET2DS(NX, NY) ;
//  Solve Ax=b by CG
    ITER = NX*NY ;
    T1 = clock() ;
    F1 = 12 ;
    F2 = 19 ;
    if(k == 1) {
//   Modified symmetric solution for Y-axis
       IERR = CG2DM(NX, NY, EPS, &ITER, &ERR) ; }
     else { if(k == 2) {
//   Symmetric solution for Y-axis
       IERR = CG2DS(NX, NY, EPS, &ITER, &ERR) ; }
      else {
//   Unsymmetric solution for Y-axis
       IERR = CG2D2(NX, NY, EPS, &ITER, &ERR) ; }
     }
    T2 = clock() ;
//  Write
    CPU  = (double)(T2 - T1)/CLOCKS_PER_SEC ;
    FLOP = (NX-1)*(NY-1)*(F1+F2*ITER)*1.0e-6 ;
    FLOPS = 0.0 ;
    if(CPU > 0.0) FLOPS = FLOP/CPU ;
    fprintf(FT1,"%4d %4d %10.2e %8.2f %8.1f \n",
                  k,ITER,ERR,CPU,FLOPS) ;
    printf("%4d %4d %10.2e %8.2f %8.1f \n",k,ITER,ERR,CPU,FLOPS) ;
//  Output center-line
    OUT2D(NX, NY, 1) ;
   }
 }
//=================================================================C
void SET2DS(int NX, int NY)
//=================================================================C
//  Set A and B by 2 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
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro,   2003/08/21                     C
//        ( Hitachi Ltd. and Waseda University )                   C
//=================================================================C
{ int    i, j ;
  double DX, DY, HX, HY, F ;
//  Clear A,B,X
  for (j=0; j<=NY; j++) {
    for (i=0; i<=NX; i++) {
      A(i,j,1) = 0.0 ;
      A(i,j,2) = 0.0 ;
      A(i,j,3) = 0.0 ;
      B(i,j)   = 0.0 ;
      X(i,j)   = 0.0 ; }
    }
//  Initial Data
  DX = NX ;
  DY = NY ;
  HX = 1.0/DX ;
  HY = 1.0/DY ;
  F  = 100.0 ;
//  Set A,B
  for (j=1; j<=NY-1; j++) {
    for (i=1; i<=NX-1; i++) {
      A(i,j,1) = -HX*DY ;
      if(j == 1) A(i,j,1) = 0.0 ;
      A(i,j,2) = -HY*DX ;
      if(i == 1) A(i,j,2) = 0.0 ;
      A(i,j,3) = 2.0*(HY*DX+HX*DY) ;
      B(i,j)   = HX*HY*F ;  }
    }
 }
//=================================================================C
void OUT2D(int NX, int NY, int ID)
//-----------------------------------------------------------------C
{ int    i, j, NO, *a, *b ;
  double DX, XV ;
//  Output Y
  fprintf(FT1,"Output on Y-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)) ; }
//  Output Un-symmetric Solution
  if(ID == 1) {
    fprintf(FT1," ** Un-symmetrix Solution ** \n") ; 
    fprintf(FT1,"   i   j        X(i,j)           X(i,NY-j) \n") ;
    NO = 0 ;
    for (j=1; j<=NY/2-1; j++) {
     for (i=1; i<=NX/2-1; i++) {
      if(X(i,j) != X(i,NY-j)) {
        a = &X(i,j) ;  b = &X(i,NY-j) ;
//  Little Endian,  (Pentium)
        fprintf(FT1,"%4d %4d  %08X%08X  %08X%08X \n",i,j,
                    *(a+1),*a,*(b+1),*b) ;
//  Big Endian,  (SR8000) 
//      fprintf(FT1,"%4d %4d  %08X%08X  %08X%08X \n",i,j,
//                  *a,*(a+1),*b,*(b+1)) ;
        NO = NO + 1 ;  }
	 }  }
     printf("Un-symmetric solution numbers= %d \n",NO) ;
     fprintf(FT1,"Un-symmetric solution numbers= %d \n",NO) ;
   }
}