ブロックLU分解CプログラムNO.4
2003/09/08 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをブロックガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
部分軸交換付で、軸交換は枢軸の右側だけを交換する。
ブロックサイズ(MB×MB)の単位に主力部の計算をするキャッシュ対策が組み込まれている
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 ;
// Function define
void BMULT(double *AA, double *BB, double *CC, int N, int MB) ;
//=================================================================C
int BLU4(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, KK1, NLG ;
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) {
KK1 = kk + MB ;
NLG = N - KK1 ;
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) )
if(NLG > 1) {
BMULT(&A[KK1][kk], &A[kk][KK1], &A[KK1][KK1], NLG, MB) ; }
}
//--- 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) ;
}
//=================================================================C
void BMULT(double A[][ND], double B[][ND], double C[][ND], int N, int MB)
//=================================================================C
// C=C-A*B by Block Multiplication with MB Block Size C
//-----------------------------------------------------------------C
// A[N][ND] R*8, In, 1st Input Matrix, Size=N*MB C
// B[MB][ND] R*8, In, 2nd Input Matrix, Size=MB*N C
// C[N][ND] R*8, I/O, Input/Output Matrix, Size=N*N C
// N I*4, In, Column Size of Matrix A C
// MB I*4, In, Block Size and Row Size of Matrix A C
//-----------------------------------------------------------------C
// Written by Yasunori Ushiro and Shouko Sakuma, 2003/09/08 C
// ( Hitachi Ltd. and Waseda University ) C
//=================================================================C
{ int i, j, k, ii, jj, ip, MC, MR ;
double AIK ;
// Outer Loop
for (ii=0; ii<N; ii=ii+MB) {
MC = MB ;
if(N-ii < MB) { MC = N - ii ; }
// Main Loop
for (jj=0; jj<N; jj=jj+MB) {
MR = MB ;
if(N-jj < MB) { MR = N - jj ; }
// C=C-A*B
for (i=0; i<MC; i++) {
ip = i + ii ;
for (k=0; k<MB; k++) {
AIK = A[ip][k] ;
for (j=0; j<MR; j++) {
C[ip][j+jj] = C[ip][j+jj] - AIK*B[k][j+jj] ; }
} }
} }
}