LU分解CプログラムNO.2
2007/05/30 東京工芸大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
特異性のチェックだけで軸交換はしていないプログラムである。
2. プログラム
#include <stdio.h>
#include <math.h>
//================================================================C
int GLU2(double *A, double *B, int N, int ND, double EPS)
//================================================================C
// Real-Dense LU Decomposition by Gauss Elimination C
// A ---> LU Decomposition without Partial Pivoting C
// with Singularity Check 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
// return I*4, Out, 0 : Normal Execution C
// 1 : Singular Stop C
// 2 ; Parameter Error C
//----------------------------------------------------------------C
// Written by Yasunori Ushiro , 2007/05/30 C
// ( Tokyo Polytechnic University ) C
// Ver.1 : No tuning Version ; Address is equal to Fortran C
//================================================================C
{ int i, j, k, IER, KND, JND ;
//----- Gauss Elimination -----------------------------
// Check Parameter
IER = 0 ;
if(ND < N)
{ IER = 2 ;
goto M100 ; }
// Gauss Elimination
for (k=0; k<N; k++) {
KND = k*ND ;
// Check Singularity
if(fabs(A[k+KND]) < EPS)
{ IER = 1 ;
goto M100 ; }
// A(*,k)=A(*,k)/A(k,k)
for (i=k+1; i<N; i++)
{ A[i+KND] = A[i+KND]/A[k+KND] ; }
// Main Elimination
for (j=k+1; j<N; j++) {
JND = j*ND ;
for (i=k+1; i<N; i++)
{ A[i+JND] = A[i+JND] - A[k+JND]*A[i+KND] ; }
}
}
//------ Solve Ax=b by Substitution with LUx=b ------
// Forward Substitution
for (j=0; j<N-1; j++) {
for (i=j+1; i<N; i++)
{ B[i] = B[i] - B[j]*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) ;
}