//================================================================C
// Test Program of Solver for Dense-Real Linear Equation Ax=b     C
//   by Gauss Elimination with LU Decomposition                   C
//----------------------------------------------------------------C
//    Written by Yasunori Ushiro ,  2007/05/30                    C
//        ( Tokyo Polytechnic University )                        C
//     後 保範（東京工芸大学）                                    C
//================================================================C
#include <stdio.h> 
#include <math.h>
#include <time.h>
// Global Define
#define NND  1000
double A[NND*NND], B[NND], WA[NND*NND], WB[NND] ;
int    IP[NND] ; 
FILE   *FT1 ;
// Function Definition
void DRMAT(double *AA, double *BB, int N, int ND, int IDA,
          int IDB, double SIG) ;
void  GLU1(double *AA, double *BB, int N, int ND) ;
int  GLU2(double *AA, double *BB, int N, int ND, double EPS) ;
int  GLU3(double *AA, double *BB, int N, int ND, double EPS, int *IP) ;
int  GLU4(double *AA, double *BB, int N, int ND, double EPS, int *IP) ;
// Main Program
void main ()
{ int    i, j, k, IDA, IDB, IER, T1, T2, N, ND, JND ; 
  double EPS, SIG, FLOP, CPU, PRC, SUM ; 
//  Open File
  ND  = NND ;
  FT1 = fopen("C-GLU.txt","w") ;
  T1  = clock() ;                    // Start Time
  EPS = 1.0e-14 ;
//  Size Input
  printf("Type In  Matrix Size(N<=1000) \n") ; 
  scanf("%d",&N) ;
  if(N > ND) { N=ND ; }
  fprintf(FT1,"  Solve Ax=b by Gauss Elimination  N=%d \n",N) ;
  fprintf(FT1,"Type. IER   Precision  Time(s)  MFLOPS \n") ;
//  Set WA,WB
  SIG = 0.15 ;
  IDA = 1 ;
  IDB = 0 ;
  DRMAT(WA, WB, N, ND, IDA, IDB, SIG) ;
//  Loop for Random Test Case
  for (k=1; k<=4; k++) {
//  Set A,B
    for (j=0; j<N; j++)
	 { B[j] = WB[j] ;
       JND  = j*ND ;
       for (i=0; i<N; i++)
        { A[i+JND] = WA[i+JND] ; }
	  }
//  LU Decomposition and Solve Ax=b
    T1 = clock() ;
    if(k <= 2) {
//  No Pivoting
      if(k <=1 ) {
        IER = 0 ;
        GLU1(A, B, N, ND) ; }
      else {
        IER = GLU2(A, B, N, ND, EPS) ; }
    }
    else {
//  Partial Pivoting
      if(k <=3 ) {
        IER = GLU3(A, B, N, ND, EPS, IP) ; }
      else {
        IER = GLU4(A, B, N, ND, EPS, IP) ; }
    }		
    T2 = clock() ;
//  Output
    CPU = (double)(T2 - T1)/CLOCKS_PER_SEC ; 
    FLOP = 0.0 ;
    if(CPU > 0.0) { FLOP = 2.0*N/3.0*N*N*1.0e-6/CPU ; }
    PRC = 0.0 ;
    SUM = 0.0 ;
    for (i=0; i<N; i++)
     { PRC = PRC + (B[i] - 1.0)*(B[i] - 1.0) ;
       SUM = SUM + B[i]*B[i] ;  }
    PRC = sqrt(PRC/SUM) ;
    fprintf(FT1,"%3d %3d %10.2e %8.1f %8.1f \n",
                  k,IER,PRC,CPU,FLOP) ;
    printf("%3d %3d %10.2e %8.1f %8.1f \n",k,IER,PRC,CPU,FLOP) ;
  } }