2次元FDM用CG法(C)テストプログラム
2003/08/24 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求めるプログラムのテスト用プログラム
CG法は比較的直交性に強いバージョンと演算量の少ないバージョンの両方をテストする。
2. プログラム
//=================================================================C
// Test Program of CG for 2-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 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 CG2D1(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-CG2D.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 ;
fprintf(FT1," CG Method for 2-Dimensional FDM, NX=%d, NY=%d \n",NX,NY) ;
fprintf(FT1," Type Loop Error Time(s) MFLOPS \n") ;
for (k=1; k<=2; k++) {
// Set A,B,X
SET2DS(NX, NY) ;
// Solve Ax=b by CG
ITER = NX*NY ;
T1 = clock() ;
if(k == 1) {
IERR = CG2D1(NX, NY, EPS, &ITER, &ERR) ;
F1 = 10 ;
F2 = 23 ; }
else {
IERR = CG2D2(NX, NY, EPS, &ITER, &ERR) ;
F1 = 12 ;
F2 = 19 ; }
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.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
OUT2D(NX, NY) ;
}
}
//=================================================================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)
//-----------------------------------------------------------------C
{ int i ;
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)) ; }
}