ブロックLU分解CプログラムNO.1
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
void BLU1(int N, int MB)
//=================================================================C
// Real-Dense LU Decomposition by Block Gauss Elimination C
// and Solve Ax=b by Substitution C
// A ---> LU Decomposition without Partial Pivoting C
//-----------------------------------------------------------------C
// A,B,ND are global define C
// N I*4, In, Matrix Order of A C
// MB I*4, In. Block Size 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, MBK ;
//----- Block Gauss Elimination Step ------------
// 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++) {
// Lower-Block Elimination
for (i=k+1; i<N; i++) {
// A(*,k)=A(*,k)/A(k,k)
A[i][k] = A[i][k]/A[k][k] ;
for (j=k+1; j<=MBK; j++) {
A[i][j] = A[i][j] - A[i][k]*A[k][j]; }
} }
// Block Process-2 ( Upper-Block Elimination )
for (i=kk+1; i<=MBK; i++) {
for (j=kk; j<i; j++) {
for (k=MBK+1; k<N; k++) {
A[i][k] = A[i][k] - A[j][k]*A[i][j] ; }
} }
// Block Process-3 ( Main Elimination by C=C-A*B) )
for (i=MBK+1; i<N; i++) {
for (k=kk; k<=MBK; k++) {
for (j=MBK+1; j<N; j++) {
A[i][j] = A[i][j] - A[i][k]*A[k][j] ; }
} }
}
//--- Solve LUx=b by Substitution -------------
// Forward Substitution
for (i=0; i<N; i++) {
for (j=0; j<i; j++) {
B[i] = B[i] - B[j]*A[i][j] ; }
}
// Backword Substitution
for (i=N-1; i>=0; i--) {
for (j=i+1; j<N; j++) {
B[i] = B[i] - A[i][j]*B[j] ; }
B[i] = B[i]/A[i][i] ;
}
}