LU分解CプログラムNO.4
2003/08/18 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
部分軸交換付きで軸交換は対象行のうち右側の行だけ交換するプログラム。
2. プログラム
#include <stdio.h>
#include <math.h>
//================================================================C
int GLU4(double *A, double *B, int N, int ND, double EPS, int *IP)
//================================================================C
// Real-Dense LU Decomposition by Gauss Elimination C
// A ---> LU Decomposition with Partial Pivoting C
// Type-2 Partial Pivoting (Changing partial rows) C
//----------------------------------------------------------------C
// A(ND*N) R*8, I/O, Coefficient Matrix C
// B(N) R*8, I/O, Righ-hand Vector C
// N I*4, In, Matrix Size of A C
// ND I*4, In, Array Size of A ( ND >= N ) C
// EPS R*8, In, Value for Singularity Check C
// IP(N) I*4, Out, Pivot Number C
// return I*4, Out, 0 : Normal Execution C
// 1 : Singular Stop C
// 2 ; Parameter Error C
//----------------------------------------------------------------C
// Written by Yasunori Ushiro , 2003/08/18 C
// ( Waseda University and Hitachi Ltd. ) C
// Ver.1 : No tuning Version ; Address is equal to Fortran C
//================================================================C
{ int i, j, k, KPIV, IER, KND, JND ;
double PIV, AKJ, DPIV, BW ;
//----- Gauss Elimination -----------------------------
// Check Parameter
IER = 0 ;
if(ND < N)
{ IER = 2 ;
goto M100 ; }
// Gauss Elimination
for (k=0; k<N; k++) {
KND = k*ND ;
// Search Maximum Value in k's column
KPIV = k ;
PIV = fabs(A[k+KND]) ;
for (i=k+1; i<N; i++) {
if(fabs(A[i+KND]) >= PIV)
{ KPIV = i ;
PIV = fabs(A[i+KND]) ; }
}
// Check Singularity
if(PIV < EPS)
{ IER = 1 ;
goto M100 ; }
IP[k] = KPIV ;
// Change A(k,k) <--> A(KPIV,k)
if(KPIV != k)
{ AKJ = A[k+KND] ;
A[k+KND] = A[KPIV+KND] ;
A[KPIV+KND] = AKJ ; }
// Pivot Value
DPIV = 1.0/A[k+KND] ;
A[k+KND] = DPIV ;
// A(*,k)=A(*,k)/A(k,k)
for (i=k+1; i<N; i++)
{ A[i+KND] = A[i+KND]*DPIV ; }
// Main Elimination
for (j=k+1; j<N; j++) {
JND = j*ND ;
// Change A(k,j) <--> A(KPIV,j)
AKJ = A[KPIV+JND] ;
A[KPIV+JND] = A[k+JND] ;
A[k+JND] = AKJ ;
// Gauss Elimination
for (i=k+1; i<N; i++)
{ A[i+JND] = A[i+JND] - AKJ*A[i+KND] ; }
}
}
//------ Solve Ax=b by Substitution with LUx=b ------
// Interchange Entries of B
for (j=0; j<N-1; j++) {
k = IP[j] ;
BW = B[k] ;
B[k] = B[j] ;
B[j] = BW ;
// Forward Substitution
for (i=j+1; i<N; i++)
{ B[i] = B[i] - BW*A[i+j*ND] ; }
}
// Back Substitution
for (j=N-1; j>=0; j--) {
B[j] = B[j]*A[j+j*ND] ;
for (i=0; i<=j-1; i++)
{ B[i] = B[i] - A[i+j*ND]*B[j] ; }
}
M100: return (IER) ;
}