/********************************************************
* Inversion of a real square matrix by LU decomposition *
* with dynamic allocations                              *
*                                                       *
*                     C++ version by J-P Moreau, Paris  *
*                            (www.jpmoreau.fr)          *
* ----------------------------------------------------- *
* SAMPLE RUN:                                           *
*                                                       *
* Input file (inv_lu.dat):                              *
*                                                       *
*  4                                                    *
*  8  2    3  12                                        *
*  2  4    7   0.25                                     *
*  3  7    3   5                                        *
* 12  0.25 5   2                                        *
*                                                       *
* Output file (inv_lu.lst):                             *
*                                                       *
* ----------------------------------------------------- *
*  INVERSION OF A REAL SQUARE MATRIX:                   *
* ----------------------------------------------------- *
*  N=4                                                  *
*                                                       *
*  8.000000  2.000000  3.000000  12.00000               *
*  2.000000  4.000000  7.000000  0.250000               *
*  3.000000  7.000000  3.000000  5.000000               *
*  12.00000  0.250000  5.000000  2.000000               *
*                                                       *
*  Inverted matrix Y:                                   *
*                                                       *
* -0.040155 -0.085788  0.056557  0.110262               *
* -0.085788 -0.062018  0.202201  0.016978               *
*  0.056557  0.202201 -0.130316 -0.038826               *
*  0.110262  0.016978 -0.038826 -0.066631               *
*                                                       *
*  Verification A*Y = I:                                *
*                                                       *
*  1.000000 -0.000000  0.000000  0.000000               *
* -0.000000  1.000000 -0.000000  0.000000               *
*  0.000000  0.000000  1.000000 -0.000000               *
*  0.000000  0.000000 -0.000000  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 *);
  

  /******************************************                                     
  *       MULTIPLY TWO REAL MATRICES        *
  * --------------------------------------- *                                     
  * INPUTS:    A  MATRIX N*N                *                                     
  *            B  MATRIX N*N                *                                     
  *            N  INTEGER                   *                                     
  * --------------------------------------- *                                     
  * OUTPUT:    C  MATRIX N*N, PRODUCT A*B   *                                     
  *                                         *                                     
  ******************************************/
  void MatMult(REAL **A, REAL **B, REAL **C, int N)  { 
    REAL SUM;
    int  I,J,K;
    for (I=1; I<=N; I++) {
	  for (J=1; J<=N; J++)  {
        SUM=0.0;                                                                
        for (K=1; K<=N; K++)
          SUM += A[I][K]*B[K][J];                                               
        C[I][J]=SUM;                                                            
      }                                                                   
    }                                                                     
  }


int main()  {

  REAL **A;       // matrix n+1 x n+2
  REAL **A1;      // copy of matrix A
  REAL **Y;       // vector n+1 
  REAL *temp;     // vector n+1 
  int  *INDX;     // vector 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,"Inv_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 matrix A

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

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

  WriteHead(fp2,"  INVERSION OF A REAL SQUARE MATRIX:");
  fprintf(fp2,"  N=%d\n\n",n);

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

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

  //call solver if previous return code is ok
  //to obtain inverse of A one column at a time
  if (rc==0) {
	for (j=1; j<n+1; j++) {
      for (i=1; i<n+1; i++) temp[i]=Y[i][j];
      LUBKSB(A,n,INDX,temp);
      for (i=1; i<n+1; i++) Y[i][j]=temp[i];
    }
  }
  //the inverse matrix is now in matrix Y
  //the original matrix A is destroyed

  if (rc==1)
    fprintf(fp2,"\n  The system matrix is singular, no solution !\n");
  else  {
    fprintf(fp2,"\n  Inverted matrix Y:\n");
    WriteMat1 (fp2, n, n, Y);
  }

  // verify A1 x Y = I (result put in A)
  MatMult(A1,Y,A,n);
  // A should now contain identity matrix
  fprintf(fp2,"\n  Verification A*Y = I:\n");
  WriteMat1 (fp2, n, n, A);

  WriteEnd(fp2);
  fclose(fp2);
  printf("\n See results in file %s.\n\n",output);
  
  return 0;
} // main

// End of file inv_lu.cpp