3次元FDM用CG法(C)テストプログラム
2003/08/26 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
3次元差分法で離散化した疎行列を係数とする連立一次方程式Ax=bに対してCG法で
反復解xを求めるプログラムのテスト用プログラム。
計算対象場がy-zで対称なら、CG法で計算解xも対称性が保たれるかのテストプログラム。
2. プログラム
//=================================================================C
// Test Program of CG for 3-Dimensional FDM with symmetric C
// for symmetric solution on Y-Z-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 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 ID) ;
int CG3DM(int NX, int NY, int NZ, double EPS, int *ITER, double *ERR) ;
int CG3DS(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-CG3DS.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," Type Loop Error Time(s) MFLOPS \n") ;
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 3-Dimensional FDM, NX=%d, NY=%d, NZ=%d \n",
NX,NY,NZ) ;
printf(" CG Method for 3-Dimensional FDM, NX=%d, NY=%d, NZ=%d \n",
NX,NY,NZ) ;
SET3DS(NX, NY, NZ) ;
// Solve Ax=b by CG
ITER = NX*NY*NZ ;
T1 = clock() ;
F1 = 16 ;
F2 = 23 ;
if(k == 1) {
// Modified symmetric solution
IERR = CG3DM(NX, NY, NZ, EPS, &ITER, &ERR) ; }
else { if(k == 2) {
// Symmetric solution
IERR = CG3DS(NX, NY, NZ, EPS, &ITER, &ERR) ; }
else {
// Unsymmetric solution
IERR = CG3D2(NX, NY, NZ, EPS, &ITER, &ERR) ; }
}
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,"%3d %5d %10.2e %8.1f %8.1f \n",
k,ITER,ERR,CPU,FLOPS) ;
printf("%3d %5d %10.2e %8.1f %8.1f \n",k,ITER,ERR,CPU,FLOPS) ;
// Output center-line
OUT3D(NX, NY, NZ, 1) ;
}
}
//=================================================================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, int ID)
//-----------------------------------------------------------------C
{ int i, j, k, NO, *a, *b, *c, *d ;
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)) ; }
// Output Un-symmetric Solution
if(ID == 1) {
fprintf(FT1," ** Un-symmetrix Solution ** \n") ;
fprintf(FT1," i j k X(i,j,k) X(i,NY-j,k) "
" X(i,j,NZ-k) X(i,NY-j,NZ-k)\n") ;
NO = 0 ;
for (k=1; k<=NZ/2-1; k++) {
for (j=1; j<=NY/2-1; j++) {
for (i=1; i<=NX/2-1; i++) {
if(X(i,j,k) != X(i,NY-j,k) || X(i,j,k) != X(i,j,NZ-k)
|| X(i,j,k) != X(i,NY-j,NZ-k) ) {
a = &X(i,j,k) ; b = &X(i,NY-j,k) ;
c = &X(i,j,NZ-k) ; d = &X(i,NY-j,NZ-k) ;
// Little Endium (pentium)
fprintf(FT1,"%3d%3d%3d %08X%08X %08X%08X %08X%08X %08X%08X\n",
i,j,k,*(a+1),*a,*(b+1),*b,*(c+1),*c,*(d+1),*d ) ;
// Big Endium (SR8000)
// fprintf(FT1,"%3d%3d%3d %08X%08X %08X%08X %08X%08X %08X%08X\n",i,j,k,
// i,j,k,*a,*(a+1),*b,*(b+1),*c,*(c+1),*d,*(d+1) ) ;
NO = NO + 1 ; }
} } }
printf("Un-symmetric solution numbers= %d \n",NO) ;
fprintf(FT1,"Un-symmetric solution numbers= %d \n",NO) ;
}
}