/************************************************************************************
*           Inversion of a complex square matrix by LU decomposition                *
*           with static allocations                                                 *
*                                                                                   *
*                                                 C++ version by J-P Moreau, Paris  *
*                                                        (www.jpmoreau.fr)          *
* --------------------------------------------------------------------------------- *
* SAMPLE RUN:                                                                       *
*                                                                                   *
* Input file (inv_clu.dat):                                                         *
*                                                                                   *
*  4                                                                                *
*  8  0  2    0  3  0  12    1                                                      *
*  2  0  4    0  7  0   0.25 2                                                      *
*  3  0  7    0  3  0   5    3                                                      *
* 12  0  0.25 0  5  0   2    4                                                      *
*                                                                                   *
* Output file (inv_clu.lst):                                                        *
*                                                                                   *
* INVERSION OF A COMPLEX SQUARE MATRIX:                                             *
*                                                                                   *
* N=4                                                                               *
*                                                                                   *
*  8.00   0.00   2.00   0.00   3.00   0.00  12.00   1.00                            *
*  2.00   0.00   4.00   0.00   7.00   0.00   0.25   2.00                            *
*  3.00   0.00   7.00   0.00   3.00   0.00   5.00   3.00                            * 
* 12.00   0.00   0.25   0.00   5.00   0.00   2.00   4.00                            *
*                                                                                   *
* Inverted matrix Y:                                                                *
*                                                                                   *
*1 -0.03022 -0.04162    -0.08426 -0.00641     0.05306  0.01466     0.10426  0.02515 *
*2 -0.07421 -0.04847    -0.06024 -0.00746     0.19813  0.01707     0.00998  0.02929 *
*3  0.05443  0.00890     0.20187  0.00137    -0.12957 -0.00313    -0.03754 -0.00538 *
*4  0.10431  0.02491     0.01606  0.00384    -0.03673 -0.00877    -0.06304 -0.01505 *
*                                                                                   *
*                                                                                   *
* Verification A*Y = I:                                                             *
*                                                                                   *
*1  1.00000  0.00000    -0.00000 -0.00000     0.00000  0.00000     0.00000  0.00000 *
*2  0.00000 -0.00000     1.00000  0.00000     0.00000  0.00000     0.00000  0.00000 *
*3 -0.00000 -0.00000    -0.00000 -0.00000     1.00000  0.00000     0.00000  0.00000 *
*4  0.00000  0.00000     0.00000  0.00000    -0.00000  0.00000     1.00000  0.00000 *
*                                                                                   *
* Uses: clu1.cpp                                                                    *
************************************************************************************/
#include <stdio.h>
#include <string.h>

#define SIZE 25

  //complex number
  typedef struct {
     double R,I;  //algebraic form
  } COMPLEX; 

  typedef COMPLEX CMAT[SIZE][SIZE];
  typedef COMPLEX CVEC[SIZE];
  typedef int IMAT[SIZE];


//Functions of clu1.cpp called by main
int  LUDCMP(CMAT, int, IMAT, int *);
void LUBKSB(CMAT, int, IMAT, CVEC);
void CMUL(COMPLEX, COMPLEX, COMPLEX *);

// read n elements of a complex vector from an input text file
  void ReadCVec (FILE *fp, int n, CVEC x)  {
    int i;
    COMPLEX tmp;

    for (i = 1; i < n+1; i++) {
      fscanf(fp,"%lf %lf", &tmp.R, &tmp.I);
      x[i].R = tmp.R; x[i].I = tmp.I;
    }
  }

// write n elements of a complex vector to an output text file
  void WriteCVec(FILE *fp, int n, CVEC x)  {
    int i;

    for (i = 1; i < n+1; i++)
      fprintf(fp," %6.2f %6.2f", x[i].R, x[i].I);
    fprintf(fp,"\n");
  }

// write m x n elements of a complex matrix to an output text file
  int WriteCMat (FILE *fp, int m, int n, CMAT a) {
    int i, j;

    if (fprintf (fp,"\n") <= 0) return (-1);

    for (i = 1; i < m+1; i++) {
      for (j = 1; j < n+1; j++)
        if (j==1) 
          fprintf(fp," %d %8.5f %8.5f",i, a[i][j].R, a[i][j].I);
        else
          fprintf(fp,"    %8.5f %8.5f", a[i][j].R, a[i][j].I);
      if (fprintf (fp,"\n") <= 0) return (-1);
	}
    if (fprintf (fp,"\n") <= 0) return (-1);

    return (0);
  }   

  /******************************************                                     
  *      MULTIPLY TWO COMPLEX MATRICES      *
  * --------------------------------------- *                                     
  * INPUTS:    A  MATRIX N*N                *                                     
  *            B  MATRIX N*N                *                                     
  *            N  INTEGER                   *                                     
  * --------------------------------------- *                                     
  * OUTPUT:    C  MATRIX N*N, PRODUCT A*B   *                                     
  *                                         *                                     
  ******************************************/
  void MatCMul(CMAT A, CMAT B, CMAT C, int N)  { 
    COMPLEX SUM, tmp;
    int  I,J,K;
    for (I=1; I<=N; I++) {
	  for (J=1; J<=N; J++)  {
        SUM.R=0.0; SUM.I=0.0;                                                                
        for (K=1; K<=N; K++) {
          CMUL(A[I][K],B[K][J],&tmp);
          SUM.R = SUM.R + tmp.R;
          SUM.I = SUM.I + tmp.I;
		}  
        C[I][J].R=SUM.R;
		C[I][J].I=SUM.I;
      }                                                                   
    }                                                                     
  }


int main()  {

  CMAT A;        // complex matrix n+1 x n+2
  CMAT A1;       // copy of matrix A
  CMAT Y;        // complex matrix n+1 x n+1 
  CVEC temp;     // vector n+1 
  IMAT INDX;     // vector n+1 

// NOTA: index zero not used here. 

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

  printf("\n Input file name (without <.dat>): ");
  scanf("%s", s);
  strcpy(s,"inv_clu");
  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

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

  for (i=1; i<n+1; i++)  {
    ReadCVec(fp1,n,temp);
    for (j=1; j<n+1; j++) {
	  A[i][j].R = temp[j].R;
      A[i][j].I = temp[j].I;
      A1[i][j].R = temp[j].R;
      A1[i][j].I = temp[j].I;
      Y[i][j].R = 0.0;
      Y[i][j].I = 0.0;
    }
	Y[i][i].R = 1.0;
    Y[i][i].I = 0.0;
    WriteCVec(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");
    WriteCMat (fp2, n, n, Y);
  }

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

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

// End of file inv_clu.cpp