#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <conio.h>

#define MAX 20


//************ complex number handling ********************
typedef struct 
{
    double r,i;
	
}doublecomplex;

#define  COMPLEXMAG(x,y,z) z=sqrt( pow( x,2) + pow(y,2) );
#define COMPLEXMAG2(x,y,z) z=    ( pow( x,2) + pow(y,2) );


#define COMPLEXRECIP(ar, ai, xr, xi)              \
{                                                 \
           xr=  ar/( pow( ar,2) + pow(ai,2) );    \
           xi= -ai/( pow( ar,2) + pow(ai,2) );    \
}                                                 \

#define COMPLEXMUL(ar, ai, br, bi,cr,ci)          \
{                                                 \
           cr=  ar*br -ai*bi;                     \
           ci=  ar*bi +ai*br;                     \
}                                                 \

#define COMPLEXDIV(nr, ni, dr, di,cr,ci)          \
{                                                 \
           cr =  nr*dr +ni*di;                    \
           cr/= ( pow( dr,2) + pow(di,2) );       \
           ci =  ni*dr -nr*di;                    \
           ci/= ( pow( dr,2) + pow(di,2) );       \
}                                                 \

//*********************************************************


void dataentry(void);
void reserve(void);
void initialize(void);
void inputdim(void);
void printresults(void);
void erroranalysis(void);

void deter(long *nd, doublecomplex w[][MAX], double *ddetr,double *ddeti);


double ah;
double aq;
double aw;

doublecomplex aa,az,ba,c,temp;
doublecomplex a[MAX][MAX],b[MAX][MAX],x[MAX][MAX];
doublecomplex ac[MAX][MAX],det;

long i,ii,ik,imethod;
long j,jj,k,kk,n,nv;


char filename[5][FILENAME_MAX];
FILE *pFile[5];

//   Variables
//
//   n is the dimension of the input matrix.
//
//	 a is the input matrix.
//
//	 x is initialized as the identity matrix.  
//    It is transformed to the inverse matrix.
//
//   b is the product of matrix a and the inverse matrix x.
//     It should equal the identity matrix.
//

void main()
{
		 printf("\n\n Inversion of a Complex Matrix via Gaussian Elimination with Partial Pivoting. \n");
		 printf("\n Version 1.0,  May 17, 2002 \n");
		 printf("\n by Tom Irvine ");
	     printf("\n Email:  tomirvine@aol.com \n");
		 
		 inputdim();
		 initialize();
		 dataentry();
	     reserve();

		 deter(&n,a,&det.r,&det.i);
		
		 if((fabs(det.r)<= 1.0e-200) && (fabs(det.i)<= 1.0e-200) )
		 {
		    printf("\n\n Error:  determinant of input matrix  = (%12.5e , %12.5e)\n",det.r,det.i);
			exit(1);
		 }
		 else
		 {
			 printf("\n\n Determinant of input matrix  = (%12.5e , %12.5e)\n",det.r,det.i);
		 }

/****** begin forward elimination  *********/

         for(i=0;i<n;++i)
         {

			 COMPLEXMAG(a[i][i].r,a[i][i].i,aq)

	          for(jj=i+1;jj<n;++jj)
	          {

			    COMPLEXMAG(a[jj][i].r,a[jj][i].i,ah)

	            if(aq<ah)
	            {
			        COMPLEXMAG(a[jj][i].r,a[jj][i].i,aq)

					for(kk=0;kk<n;++kk)
					{
						temp.r=a[jj][kk].r;
						temp.i=a[jj][kk].i;

						a[jj][kk].r=a[i][kk].r;
						a[jj][kk].i=a[i][kk].i;

						a[i][kk].r=temp.r;		
						a[i][kk].i=temp.i;						
					}
					for(nv=0; nv<n; nv++)
					{
						temp.r=x[jj][nv].r;
						temp.i=x[jj][nv].i;						
						
						x[jj][nv].r=x[i][nv].r;
						x[jj][nv].i=x[i][nv].i;				
						
						
						x[i][nv].r=temp.r;
						x[i][nv].i=temp.i;
					}	 
	            }
			  }
	          ik=i+1;
	          for(ii=ik;ii<n;++ii)
	          {
				  COMPLEXMAG(a[i][i].r,a[i][i].i,aw)

				  if( fabs(aw)<1.0e-200)
				  {
					  printf("\n Error:  zero in diagonal.\n");
					  exit(1);
				  }

				  COMPLEXDIV(a[ii][i].r,a[ii][i].i,a[i][i].r,a[i][i].i,aa.r,aa.i)

				  for(nv=0; nv<n; nv++)
				  {

                     COMPLEXDIV( a[ii][i].r,a[ii][i].i,a[i][i].r,a[i][i].i,az.r,az.i) 
                     COMPLEXMUL(az.r,az.i,x[i][nv].r,x[i][nv].i,ba.r,ba.i )


	                 x[ii][nv].r+=-ba.r;
	                 x[ii][nv].i+=-ba.i;


				  }
	              for(j=i;j<n;++j)
	              {

					 COMPLEXMUL(a[i][j].r,a[i][j].i,aa.r,aa.i,ba.r,ba.i );

 		                 a[ii][j].r+=-ba.r;
						 a[ii][j].i+=-ba.i;

		             if(( fabs(a[ii][j].r) <1.0e-10) && ( fabs(a[ii][j].i) <1.0e-10) )
		             {
		                 a[ii][j].r=0.;
						 a[ii][j].i=0.;
		             }
	              }
	          }
         }
		
		

/********* back substitution ***********************************/

		 for(nv=0; nv<n; nv++)
		 { 

			   COMPLEXDIV( x[n-1][nv].r,x[n-1][nv].i,a[n-1][n-1].r,a[n-1][n-1].i,ba.r,ba.i)
			   
			   x[n-1][nv].r=ba.r;
			   x[n-1][nv].i=ba.i;
			  

			   for(i=n-2;i>=0;i=i-1)
			   {              
				  for(k=i+1;k<n;++k)
	              {

					  COMPLEXMUL(x[k][nv].r, x[k][nv].i, a[i][k].r, a[i][k].i,c.r,c.i)
	                  
					  x[i][nv].r+=-c.r;
					  x[i][nv].i+=-c.i;

	              }
	             
                  COMPLEXDIV(x[i][nv].r,x[i][nv].i, a[i][i].r,a[i][i].i, ba.r,ba.i)

                  x[i][nv].r = ba.r;
				  x[i][nv].i = ba.i;
			   }
		 }

/************* end of gaussian elimination ***********************/

	     printresults();

         deter(&n,x,&det.r,&det.i);
		 printf("\n\n Determinant of inverse matrix  = (%12.5e , %12.5e) \n",det.r,det.i);

       	 erroranalysis();

	     deter(&n,b,&det.r,&det.i);
		 printf("\n Determinant of error matrix  = (%12.5e , %12.5e) \n\n",det.r,det.i);

		 if(imethod==1)
		 {
			 printf("\n Input data written to file:    inv_data.in ");
		 }
		 printf("\n Output data written to file:   inv_data.out");
         printf("\n Error matrix written to file:  inv_err.out \n");

		 printf("\n Calculation complete. \n Press any key to exit. \n\n");
		 getch();

/********************** end of main program **************************/

}
void deter(long *nd, doublecomplex w[][MAX], double *ddetr,double *ddeti)
{
	/********* calculate determinant *********/

//  The determinant of an upper triangular matrix A is the product of 
//  the diagonal elements. 

//  The determinant of a general square matrix B is defined as the determinant 
//  of that upper triangular matrix that has been obtained from B by swapping 
//  rows and adding them. 

//  One can arrive to an upper triangular matrix from a given general 
//  matrix B in several different ways, and the resulting upper triangular 
//  matrix is not unique. It can, nevertheless, be shown that the determinant 
//  does not depend on the various choices. 


	     long i,j,k;
		 doublecomplex aa,dd[MAX][MAX];

		 for(i=0; i<*nd; i++)
         {	
	        for(j=0; j<*nd; j++)
	        {
	          dd[i][j].r=w[i][j].r;
	          dd[i][j].i=w[i][j].i;
	        }
         }


	 

		 for( i=0; i< *nd; i++)
		 {
			for( k=i+1; k < *nd; k++)
			{
                if(( fabs(dd[i][i].r) < 1.0e-50) && ( fabs(dd[i][i].i) < 1.0e-50) )
				{
					*ddetr = 0;
					*ddeti = 0;
					return;
				}

				COMPLEXDIV(dd[k][i].r,dd[k][i].i,dd[i][i].r,dd[i][i].i,aa.r,aa.i )

				for( j = 0; j < *nd; j++)
				{	
                   COMPLEXMUL( aa.r,aa.i,dd[i][j].r,dd[i][j].i,ba.r,ba.i)


					dd[k][j].r-=ba.r;
					dd[k][j].i-=ba.i;
				}	
			}
		}

		*ddetr=dd[0][0].r;
		*ddeti=dd[0][0].i;

        for(i=1; i< *nd; i++)
		{
	
            COMPLEXMUL( *ddetr,*ddeti,dd[i][i].r,dd[i][i].i,ba.r,ba.i)

			*ddetr=ba.r;
			*ddeti=ba.i;
		}
/*
		printf("\n\n");

        for(i=0; i< *nd; i++)
		{
			for(j=0; j< *nd; j++)
			{
				printf(" %ld %ld  dd = ( %lf , %lf )\n",i,j,dd[i][j].r,dd[i][j].i);
			}
        }

		printf("\n\n");
*/

}

void dataentry(void)
{
	     double anum;

	     if(imethod==1)
		 {
			 strcpy(filename[0],"inv_data.in");
		     pFile[0]=fopen(filename[0],"w");
 
			 for(i=0; i< n; i++)
			 {
				 for(j=0; j<n; j++)
				 {
					 printf("\n Enter a[%ld][%ld].r \n",i+1,j+1);
					 scanf("%lf",&anum);
					 a[i][j].r=anum;


					 printf("\n Enter a[%ld][%ld].i \n",i+1,j+1);
					 scanf("%lf",&anum);
					 a[i][j].i=anum;
					 

					 if( j <(n-1) )
					 {
					    fprintf(pFile[0],"(%12.5e , %12.5e )\t",a[i][j].r,a[i][j].i);
					 }
					 else
					 {
			            fprintf(pFile[0],"(%12.5e , %12.5e )\n",a[i][j].r,a[i][j].i);
					 }
				 }
			 }
		 }
		 if(imethod==2)
		 {
			 printf( "\n Input filename.  Format is free, but no header lines allowed. \n");
			              
			 scanf("%s",filename[1]);
			 pFile[1]=fopen(filename[1], "rb");

			 if(pFile[1] == NULL)
			 {
				 printf("Input file error \n");
				 exit(1);
			 }

			 for(i=0; i< n; i++)
			 {
				 for(j=0; j<n; j++)
				 {
					 fscanf(pFile[1],"%lf",&anum);
					 a[i][j].r=anum;

					 fscanf(pFile[1],"%lf",&anum);
					 a[i][j].i=anum;
				 }	
			 }
			 
			 if(n<=6)
		     {
				printf("\n\n  Input Matrix \n\n");

				for(i=0; i<n; i++)
				{  
					for(j=0; j<n; j++)
					{              
				 		if( j <(n-1) )
						{
							printf("(%12.5e , %12.5e)\t",a[i][j].r,a[i][j].i);
						}
						else
						{
							printf("(%12.5e , %12.5e) \n",a[i][j].r,a[i][j].i);
						}
					}
				}
		    }
		 }
		 if(imethod !=1 && imethod !=2 ){exit(1);}
}
void reserve(void)
{
/****  copy data into reserve matrices  **/

         for(i=0;i<n;++i)
         {	
	        for(j=0;j<n;++j)
	        {
	          ac[i][j].r=a[i][j].r;
	          ac[i][j].i=a[i][j].i;
	        }
         }

}
void initialize(void)
{
	     for(i=0; i<MAX; i++)
	     {
			  for(j=0; j<MAX; j++)
			  {
				 a[i][j].r=0.;
				 ac[i][j].r=0.;
			     x[i][j].r=0.;
				 b[i][j].r=0.;

				 a[i][j].i=0.;
				 ac[i][j].i=0.;
			     x[i][j].i=0.;
				 b[i][j].i=0.;
			  }
			  x[i][i].r=1.;
			  x[i][i].i=0.;
	     }

}
void inputdim(void)
{
	     printf("\n The input matrix has dimension n x n. \n\n Enter n \n");
		 scanf("%ld",&n);

		 if(n>20)
		 {
			 printf("\n Error:  maximum dimension = 20 \n");
			 exit(1);
		 }
		 if(n<2)
		 {
			 printf("\n Error:  minimum dimension = 2 \n");
			 exit(1);
		 }

		 printf("\n\n Select data entry method:  1=keyboard  2=input file \n");

		 scanf("%ld",&imethod);
}
void printresults(void)
{
	     strcpy(filename[2],"inv_data.out");
		 pFile[2]=fopen(filename[2],"w");

         if(n<=6)
		 {
			 printf("\n\n Inverse Matrix \n\n");

			 for(i=0; i<n; i++)
		     {  
			   for(j=0; j<n; j++)
			   {              
				 	 if( j <(n-1) )
					 {
					    printf("(%12.5e , %12.5e)\t",x[i][j].r,x[i][j].i);
					 }
					 else
					 {
			            printf("(%12.5e , %12.5e) \n",x[i][j],x[i][j].i);
					 }

			   }
		     }

		 }
	  
	     for(i=0; i<n; i++)
		 {  
			   for(j=0; j<n; j++)
			   {              
				 	 if( j <(n-1) )
					 {
					    fprintf(pFile[2],"(%12.5e , %12.5e)\t",x[i][j].r,x[i][j].i);
					 }
					 else
					 {
			            fprintf(pFile[2],"(%12.5e , %12.5e)\n",x[i][j].r,x[i][j].i);
					 }

			   }
		 }
       

}
void erroranalysis(void)
{
	  for(i=0; i<n; i++)
	  {  
          for(j=0; j<n; j++)
		  {    
             for(k=0; k<n; k++)
		     {
				COMPLEXMUL(ac[i][k].r,ac[i][k].i,x[k][j].r,x[k][j].i,ba.r,ba.i)

                b[i][j].r+=ba.r;
                b[i][j].i+=ba.i;
		     }
          }
	  }

	  	
	     strcpy(filename[3],"inv_err.out");
		 pFile[3]=fopen(filename[3],"w");

         if(n<=6)
		 {
			 printf("\n Error Analysis: \n Input Matrix x Inverse Matrix (should equal identity matrix). \n\n");

			 for(i=0; i<n; i++)
		     {  
			   for(j=0; j<n; j++)
			   {              
				 	 if( j <(n-1) )
					 {
					    printf("(%12.5e , %12.5e)\t",b[i][j].r,b[i][j].i);
					 }
					 else
					 {
			            printf("(%12.5e , %12.5e)\n",b[i][j].r,b[i][j].i);
					 }

			   }
		     }

		 }
 
	     for(i=0; i<n; i++)
		 {  
			   for(j=0; j<n; j++)
			   {              
				 	 if( j <(n-1) )
					 {
					    fprintf(pFile[3],"(%12.5e , %12.5e)\t",b[i][j].r,b[i][j].i);
					 }
					 else
					 {
			            fprintf(pFile[3],"(%12.5e , %12.5e)\n",b[i][j].r,b[i][j].i);
					 }

			   }
		 }
        

}