#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 ,  2007/05/30                    C
//        ( Tokyo Polytechnic University )                        C
//     後 保範（東京工芸大学）                                    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) ; 
  }