/*********************************************************************
*     Solving a complex linear system AX = B by LU decomposition     *
*                                                                    *
*                                  C++ version by J-P Moreau, Paris  *
*                                        (www.jpmoreau.fr)           *
* ------------------------------------------------------------------ *
* SAMPLE RUN:                                                        *
*                                                                    *
* Input file (test_clu.dat):                                         *
*                                                                    *
* 15                                                                 *
* 0 0 0 0 6 0 1 0 1 0 4 0 2 0 1 0 2 0 2 0 0 0 3 0 0 0 5 0 0 1  1 1   *
* 6 0 2 0 5 0 2 0 4 0 5 0 5 0 2 0 1 0 2 0 3 0 1 0 5 0 1 0 3 2  2 2   *
* 6 0 2 0 5 0 6 0 3 0 6 0 5 0 0 0 0 0 1 0 3 0 0 0 4 0 0 0 5 3  3 3   *
* 5 0 4 0 3 0 1 0 4 0 4 0 6 0 4 0 6 0 1 0 4 0 2 0 0 0 5 0 3 4  4 4   *
* 3 0 4 0 6 0 4 0 2 0 0 0 4 0 6 0 5 0 5 0 1 0 1 0 6 0 5 0 4 5  5 5   *
* 4 0 2 0 0 0 4 0 1 0 6 0 1 0 1 0 3 0 6 0 1 0 1 0 1 0 1 0 4 6  6 6   *
* 1 0 1 0 4 0 5 0 5 0 0 0 0 0 0 0 0 0 5 0 5 0 2 0 4 0 3 0 5 7  7 7   *
* 2 0 5 0 5 0 0 0 0 0 0 0 0 0 1 0 4 0 1 0 3 0 1 0 2 0 5 0 2 8  8 8   *
* 4 0 4 0 6 0 2 0 3 0 2 0 0 0 5 0 5 0 3 0 5 0 2 0 3 0 4 0 6 9  9 9   * 
* 3 0 2 0 0 0 0 0 3 0 4 0 0 0 6 0 5 0 4 0 6 0 1 0 4 0 5 0 3 10 10 10 *
* 1 0 1 0 6 0 2 0 0 0 6 0 0 0 2 0 0 0 0 0 6 0 3 0 4 0 0 0 4 11 11 11 *
* 6 0 3 0 2 0 3 0 6 0 2 0 0 0 2 0 3 0 0 0 5 0 5 0 1 0 6 0 3 12 12 12 *
* 2 0 1 0 5 0 4 0 0 0 4 0 0 0 6 0 3 0 2 0 3 0 1 0 1 0 4 0 6 13 13 13 *
* 3 0 0 0 5 0 1 0 1 0 6 0 0 0 4 0 3 0 6 0 1 0 0 0 1 0 0 0 2 14 14 14 *
* 0 0 3 0 3 0 3 0 4 0 5 0 0 0 3 0 1 0 2 0 4 0 6 0 4 0 1 0 4 15 15 15 *
*                                                                    *
* Output file (test_clu.lst):                                        *
*                                                                    *
*  COMPLEX LINEAR SYSTEM TO BE SOLVED:  N=15                         *
*  --/--                                                             *
*                                                                    *
*  System solution:                                                  *
*                                                                    *
*  X1=   -0.057831     0.092402                                      *
*  X2=   -0.885623     1.415035                                      *
*  X3=   -0.106641     0.170390                                      *
*  X4=   -0.421776     0.673908                                      *
*  X5=    0.063331    -0.101190                                      *
*  X6=   -0.098265     0.157007                                      *
*  X7=    0.217855    -0.348086                                      *
*  X8=   -0.698291     1.115719                                      *
*  X9=    0.869955    -1.390002                                      * 
*  X10=   -0.338608     0.541023                                     *
*  X11=   -0.463159     0.740028                                     *
*  X12=    0.118877    -0.189941                                     *
*  X13=    0.514975    -0.822819                                     *
*  X14=    0.013271    -0.021204                                     *
*  X15=    0.731170    -1.168252                                     *
*                                                                    *
*********************************************************************/
#include <stdio.h>
#include <math.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 lu.cpp called by main
int  LUDCMP(CMAT, int, IMAT, int *);
void LUBKSB(CMAT, int, IMAT, CVEC);


  void ReadVec1 (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;
    }
  }

  void WriteVec1(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");
  }

int main()  {

  CMAT A;         // matrice n+1 x n+2
  CVEC B;         // vecteur n+1 
  CVEC temp;      // vecteur n+1 
  IMAT INDX;      // vecteur 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;

  strcpy(s,"cmat15");
  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

  fprintf(fp2,"\n  COMPLEX 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].R = temp[j].R;
      A[i][j].I = temp[j].I;
    }
    B[i].R = temp[n+1].R;
    B[i].I = temp[n+1].I;
    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 %12.6f\n",i, B[i].R, B[i].I);
  }
  fclose(fp2);
  printf("\n See results in file %s.\n\n",output);
  
  return 0;
} // main

// End of file test_clu.cpp