CのブロックLU分解法テストプログラム
2003/09/08 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをブロックガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求めるプログラムのテストプログラム。
4種類のLU分解プログラムを一回にテストする。
1000次元までの乱数データの行列で、結果の精度、計算速度及びMFLOPSを
各LU分解法ごとに出力する。計算する次元数(N)及びブロック長(MB)を画面より与える。
2. プログラム
//================================================================C
// Test Program of Solver for Dense-Real Linear Equation Ax=b C
// by Block-Gauss Elimination with LU Decomposition C
//----------------------------------------------------------------C
// Written by Yasunori Ushiro, 2003/09/04 C
// ( Waseda University and Hitachi Ltd. ) C
//================================================================C
#include <stdio.h>
#include <math.h>
#include <time.h>
// Global Define
#define ND 2000
double A[ND][ND], B[ND], WA[ND][ND], WB[ND] ;
int IP[ND] ;
FILE *FT1 ;
// Function Definition
void DMATB(int N, int IDA, int IDB, double SIG) ;
int BLU1(int N, int MB) ;
int BLU2(int N, int MB, double EPS) ;
int BLU3(int N, int MB, double EPS) ;
int BLU4(int N, int MB, double EPS) ;
// Main Program
void main ()
{ int i, j, k, N, MB, IDA, IDB, IER, T1, T2 ;
double EPS, SIG, FLOP, CPU, PRC, SUM ;
// Open File
FT1 = fopen("C-BLU.txt","w") ;
T1 = clock() ; // Start Time
EPS = 1.0e-14 ;
// Size Input
printf("Type In N(<=%d), MB(Block Size) \n",ND) ;
scanf("%d %d",&N,&MB) ;
if(N > ND) { N = ND ; }
if(MB > N) { MB = N ; }
fprintf(FT1," Solve Ax=b by Block-Gauss Elimination N=%d, MB=%d \n",N,MB) ;
fprintf(FT1,"Type. IER Precision Time(s) MFLOPS \n") ;
// Set WA,WB
SIG = 0.15 ;
IDA = 1 ;
IDB = 0 ;
DMATB(N, IDA, IDB, SIG) ;
// Loop for Random Test Case
for (k=1; k<=4; k++) {
// Set A,B
for (i=0; i<N; i++)
{ B[i] = WB[i] ;
for (j=0; j<N; j++)
{ A[i][j] = WA[i][j] ; }
}
// Block-LU Decomposition and Solve Ax=b
T1 = clock() ;
if(k <= 2) {
// No Pivoting
if(k <=1 ) {
IER = 0 ;
BLU1(N, MB) ; }
else {
IER = BLU2(N, MB, EPS) ; }
}
else {
// Partial Pivoting
if(k <=3 ) {
IER = BLU3(N, MB, EPS) ;
}
else {
IER = BLU4(N, MB, EPS) ;
}
}
T2 = clock() ;
// Output
CPU = (double)(T2 - T1)/CLOCKS_PER_SEC ;
FLOP = 0.0 ;
if(CPU > 0.0) { FLOP = 2.0*N/3.0*N*N*1.0e-6/CPU ; }
PRC = 0.0 ;
SUM = 0.0 ;
for (i=0; i<N; i++)
{ PRC = PRC + (B[i] - 1.0)*(B[i] - 1.0) ;
SUM = SUM + B[i]*B[i] ; }
PRC = sqrt(PRC/SUM) ;
fprintf(FT1,"%3d %3d %10.2e %8.1f %8.1f \n",
k,IER,PRC,CPU,FLOP) ;
printf("%3d %3d %10.2e %8.1f %8.1f \n",k,IER,PRC,CPU,FLOP) ;
} }