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

    2003/09/08 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------

1. 概要

 実密行列を係数とする連立一次方程式Ax=bをブロックガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
部分軸交換付で、軸交換は枢軸の右側だけを交換する。
ブロックサイズ(MB×MB)の単位に主力部の計算をするキャッシュ対策が組み込まれている

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 ;
// Function define
void BMULT(double *AA, double *BB, double *CC, int N, int MB) ;
//=================================================================C
int BLU4(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, KK1, NLG ;
  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) {
    KK1 = kk + MB ;
    NLG = N - KK1 ;
    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) )
      if(NLG > 1) {
        BMULT(&A[KK1][kk], &A[kk][KK1], &A[KK1][KK1], NLG, MB) ;  }
   }
//--- 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) ;
  }
//=================================================================C
void BMULT(double A[][ND], double B[][ND], double C[][ND], int N, int MB)
//=================================================================C
//  C=C-A*B by Block Multiplication with MB Block Size             C
//-----------------------------------------------------------------C
//    A[N][ND]  R*8, In,  1st Input Matrix, Size=N*MB              C
//    B[MB][ND] R*8, In,  2nd Input Matrix, Size=MB*N              C
//    C[N][ND]  R*8, I/O, Input/Output Matrix, Size=N*N            C
//    N         I*4, In,  Column Size of Matrix A                  C
//    MB        I*4, In,  Block Size and Row Size of Matrix A      C
//-----------------------------------------------------------------C
//    Written by Yasunori Ushiro and Shouko Sakuma,  2003/09/08    C
//        ( Hitachi Ltd. and Waseda University )                   C
//=================================================================C
{ int     i, j, k, ii, jj, ip, MC, MR ;
  double  AIK ; 
//  Outer Loop
  for (ii=0; ii<N; ii=ii+MB) {
    MC = MB ;
    if(N-ii < MB) { MC = N - ii ; }
//  Main Loop
    for (jj=0; jj<N; jj=jj+MB) {
      MR = MB ;
      if(N-jj < MB) { MR = N - jj ; }
//   C=C-A*B
      for (i=0; i<MC; i++) {
        ip = i + ii ;
        for (k=0; k<MB; k++) {
          AIK = A[ip][k] ;
          for (j=0; j<MR; j++) { 
            C[ip][j+jj] = C[ip][j+jj] - AIK*B[k][j+jj] ; }
       } }
    } }
  }