#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 ,  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, 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] ; }
  }
}