/******************************************************
* Solving a linear system AX = B by LU decomposition  *
*                                                     *
*                   C++ version by J-P Moreau, Paris  *
*                         (www.jpmoreau.fr)           *
* --------------------------------------------------- *
* SAMPLE RUN:                                         *
*                                                     *
* Input file (test_lu.dat):                           *
*                                                     *
*  4                                                  *
*  8  2    3  12     25.0                             *
*  2  4    7   0.25  13.25                            *
*  3  7    3   5     18.0                             *
* 12  0.25 5   2     19.25                            *
*                                                     *
* Output file (test_lu.lst):                          *
*                                                     *
* --------------------------------------------------- *
*  LINEAR SYSTEM TO BE SOLVED:                        *
* --------------------------------------------------- *
*  N=4                                                *
*                                                     *
*  8.000000  2.000000  3.000000  12.00000  25.00000   *
*  2.000000  4.000000  7.000000  0.250000  13.25000   *
*  3.000000  7.000000  3.000000  5.000000  18.00000   *
*  12.00000  0.250000  5.000000  2.000000  19.25000   *
*                                                     *
*  System solution:                                   *
*                                                     *
*  X1=  1.000000                                      *
*  X2=  1.000000                                      *
*  X3=  1.000000                                      *
*  X4=  1.000000                                      *
* --------------------------------------------------- *
* Uses: basis_r.cpp, lu.cpp, wmblock.cpp              *
******************************************************/
#include <basis.h>
#include <vmblock.h>

//Functions of lu.cpp called by main
int  LUDCMP(REAL **, int, int *, int *);
void LUBKSB(REAL **, int, int *, REAL *);
  

int main()  {

  REAL **A;       // matrice n+1 x n+2
  REAL *B;        // vecteur n+1 
  REAL *temp;     // vecteur n+1 
  int  *INDX;     // vecteur n+1 
  void *vmblock = NULL; 

// NOTA: index zero not used here. 

  int d=0, i, j, n, rc=0;
  char input[20], output[20], s[20];
  FILE *fp1, *fp2;

  strcpy(s,"Test_lu");
  strcpy(input,s);
  strcat(input,".dat");
  strcpy(output,s);
  strcat(output,".lst");

  fp1 = fopen(input,"r");
  fp2 = fopen(output,"w");

  fscanf(fp1,"%d",&n);     // size of given system AX=B

  //dynamic allocations of real matrix A, real vector B, 
  //temporary vector temp and integer vector INDX.
  vmblock = vminit();  
  A       = (REAL **)vmalloc(vmblock, MATRIX,  n+1, n+1);  
  B       = (REAL *) vmalloc(vmblock, VEKTOR,  n+1, 0);
  temp    = (REAL *) vmalloc(vmblock, VEKTOR,  n+2, 0);
  INDX    = (int *)  vmalloc(vmblock, VEKTOR,  n+1, 0);

  if (! vmcomplete(vmblock))
  {
    LogError ("No Memory", 0, __FILE__, __LINE__);
    return 1;
  }

  WriteHead(fp2,"  LINEAR SYSTEM TO BE SOLVED:");
  fprintf(fp2,"  N=%d\n\n",n);

  for (i=1; i<n+1; i++)  {
    ReadVec1(fp1,n+1,temp);
    for (j=1; j<n+1; j++) A[i][j] = temp[j];
    B[i] = temp[n+1];
    WriteVec1(fp2,n+1,temp);
  }
  fclose(fp1);

  //Call LU decomposition routine
  rc = LUDCMP(A,n,INDX,&d);

  //if result is ok (rc=0), call appropriate solver
  if (rc==0) LUBKSB(A,n,INDX,B);

  if (rc==1)
    fprintf(fp2,"\n  The system matrix is singular, no solution !\n");
  else  {
    fprintf(fp2,"\n  System solution:\n\n");
    for (i=1; i<n+1; i++)  
      fprintf(fp2,"  X%d=%12.6f\n",i,B[i]);
  }
  WriteEnd(fp2);
  fclose(fp2);
  printf("\n See results in file %s.\n\n",output);
  
  return 0;
} // main

// End of file test_lu.cpp