LU分解CプログラムNO.1
2003/08/18 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
実密行列を係数とする連立一次方程式Ax=bをガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
軸交換も特異性のチェックもしていない最も単純なプログラムである。
2. プログラム
#include <stdio.h>
#include <math.h>
//================================================================C
void GLU1(double *A, double *B, int N, int ND)
//================================================================C
// Real-Dense LU Decomposition by Gauss Elimination C
// A ---> LU Decomposition without Partial Pivoting 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
//----------------------------------------------------------------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, KND, JND ;
//----- Gauss Elimination -----------------------------
// Gauss Elimination
for (k=0; k<N; k++) {
KND = k*ND ;
// 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] ; }
}
}