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