ブロックLU分解CプログラムNO.3

    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
int BLU3(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 ;
  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) {
    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) )
    for (i=MBK+1; i<N; i++) {
      for (k=kk; k<=MBK; k++) {
        AIK = A[i][k] ;
        for (j=MBK+1; j<N; j++) {
          A[i][j] = A[i][j] - AIK*A[k][j] ;  }
     }  }
   }
//--- 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) ;
  }