ブロックLU分解CプログラムNO.3
2003/09/08 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをブロックガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
部分軸交換付で、軸交換は枢軸の右側だけを交換する。
2. プログラム
#include <stdio.h>
#include <math.h>
// Global define
#define ND 2000
extern double A[ND][ND], B[ND] ;
extern double IP[ND] ;
extern FILE *FT1 ;
//=================================================================C
int BLU3(int N, int MB, double EPS)
//=================================================================C
// Real-Dense LU Decomposition by Block Gauss Elimination C
// and Solve Ax=b by Substitution C
// A ---> LU Decomposition with Partial Pivoting C
// Type-2 Partial Pivoting (Changing partial rows) C
//-----------------------------------------------------------------C
// A,B,IP,ND are global define C
// N I*4, In, Matrix Order of A C
// MB I*4, In. Block Size C
// EPS R*8, In, Value for Singularity check C
// (return) I*4, Out, 0 : Normal Execution C
// 1 : Singular Stop C
// 2 ; Parameter Error C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro , 2003/09/08 C
// ( Hitachi Ltd. and Waseda University ) C
// Ver.1 No tuning Version C
//=================================================================C
{ int i, j, k, kk, jj, KPIV, IER, MBK, ME ;
double PIV, DPIV, AKJ, AIK, AIJ, BW, S ;
//----- Block Gauss Elimination Step ------------
// Check Parameter
IER = 0 ;
if(N > ND || MB > N)
{ IER = 2 ;
goto M100 ; }
// Gauss Elimination
for (kk=0; kk<N; kk=kk+MB) {
MBK = kk + MB - 1;
if(MBK > N-1) { MBK = N - 1 ; }
// Block Process-1
for (k=kk; k<=MBK; k++) {
// Search Maximum Value in k's column
KPIV = k ;
PIV = fabs(A[k][k]) ;
for (i=k+1; i<N; i++) {
if(fabs(A[i][k]) > PIV) {
KPIV = i ;
PIV = fabs(A[i][k]) ; }
}
// Check Singularity
if(PIV < EPS) {
IER = 1 ;
goto M100 ; }
IP[k] = KPIV ;
// Change A(k,*) <--> A(KPIV,*)
if(KPIV != k) {
for (j=kk; j<N; j++) {
AKJ = A[k][j] ;
A[k][j] = A[KPIV][j] ;
A[KPIV][j] = AKJ ; }
}
// Pivot Value
DPIV = 1.0/A[k][k] ;
A[k][k] = DPIV ;
// Lower-Block Elimination
for (i=k+1; i<N; i++) {
// A(*,k)=A(*,k)/A(k,k)
AIK = A[i][k]*DPIV ;
A[i][k] = AIK ;
for (j=k+1; j<=MBK; j++) {
A[i][j] = A[i][j] - AIK*A[k][j]; }
}
}
// Block Process-2 ( Upper-Block Elimination )
for (i=kk+1; i<=MBK; i++) {
for (j=kk; j<i; j++) {
AIJ = A[i][j] ;
for (k=MBK+1; k<N; k++) {
A[i][k] = A[i][k] - AIJ*A[j][k] ; }
} }
// Block Process-3 ( Main Elimination by C=C-A*B) )
for (i=MBK+1; i<N; i++) {
for (k=kk; k<=MBK; k++) {
AIK = A[i][k] ;
for (j=MBK+1; j<N; j++) {
A[i][j] = A[i][j] - AIK*A[k][j] ; }
} }
}
//--- Solve LUx=b by Substitution -------------
// Interchange Entries of B
for (jj=0; jj<N-1; jj=jj+MB) {
ME = jj + MB ;
if(ME > N-1) { ME = N - 1 ; }
for (i=jj; i<ME; i++) {
k = IP[i] ;
BW = B[k] ;
B[k] = B[i] ;
B[i] = BW ; }
// Forward Substitution
for (i=jj+1; i<=ME; i++) {
S = 0.0 ;
for (j=jj; j<i; j++) {
S = S + B[j]*A[i][j] ; }
B[i] = B[i] - S ;
}
for (i=ME+1; i<N; i++) {
S = 0.0 ;
for (j=jj; j<ME; j++) {
S = S + B[j]*A[i][j] ; }
B[i] = B[i] - S ;
}
}
// Backword Substitution
for (i=N-1; i>=0; i--) {
S = 0.0 ;
for (j=i+1; j<N; j++) {
S = S + A[i][j]*B[j] ; }
B[i] = (B[i] - S)*A[i][i] ;
}
M100: ;
return (IER) ;
}