/*************************************************************
*   Solving a complex homogeneous linear system by Gauss     *
*   Method with full pivoting.                               *
* ---------------------------------------------------------- *
* SAMPLE RUN:                                                *
*                        | -1  i -i |                        *
* Solve AX = 0 with  A = | -i -1  1 |                        *
*                        |  i  1 -1 |                        *
*                                                            *
* Basic Family of system solutions:                          *            
*  Solution 1                                                *
*   X1 = 0.000000 + 1.000000 I                               *
*   X2 = 1.000000 + 0.000000 I                               *
*   X3 = 0.000000 + 0.000000 I                               *
*                                                            *
*  Solution 2                                                *
*   X1 = 0.000000 - 1.000000 I                               *
*   X2 = 0.000000 + 0.000000 I                               *
*   X3 = 1.000000 + 0.000000 I                               *
*                                                            *
* ---------------------------------------------------------- *
* Ref.: "Algèbre, Algorithmes et programmes en Pascal        *
*        By Jean-Louis Jardrin, DUNOD Paris, 1988".          *
*                                                            *
*                         C++ Release By J-P Moreau, Paris.  *
*                                (www.jpmoreau.fr)           *
*************************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  20

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

  typedef COMPLEX Matc[NMAX][NMAX];

  int i,l,N,M0,R0;
  double eps;

  Matc  A, VX;


  double CABS(COMPLEX Z) {
    return sqrt(Z.R*Z.R + Z.I*Z.I);
  }

  void CADD(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) {
    Z->R = Z1.R + Z2.R;
    Z->I = Z1.I + Z2.I;
  }

  void CDIF(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) {
    Z->R = Z1.R - Z2.R;
    Z->I = Z1.I - Z2.I;
  }

  void CMUL(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) {
    Z->R = Z1.R*Z2.R - Z1.I*Z2.I;
    Z->I = Z1.R*Z2.I + Z1.I*Z2.R;
  }

  void CDIV(COMPLEX Z1, COMPLEX Z2, COMPLEX *Z) {
    double d; COMPLEX C;
    d = Z2.R*Z2.R+Z2.I*Z2.I;
    if (d>2.2e-16) { 
      C.R=Z2.R; C.I=-Z2.I;
      CMUL(Z1,C,Z);
      Z->R=Z->R/d; Z->I=Z->I/d;
    }
  }

  void CSwap(COMPLEX *Z1, COMPLEX *Z2) {
    COMPLEX  C;
    C.R=Z1->R; C.I=Z1->I; 
    Z1->R=Z2->R; Z1->I=Z2->I;
    Z2->R=C.R; Z2->I=C.I;
  }   

  /*********************************************************
  * Procedure RSHCGT solves the complex homogeneous linear *
  * system AX = 0 by Gauss method with full pivoting.      *
  * ------------------------------------------------------ *
  * INPUTS:                                                *
  *         eps:  required precision                       *
  *          N :  size of complex matrix A (NxN)           *
  * OUTPUTS:                                               *
  *          R0:  rank of A in output                      *
  *          M0:  dimension of system solutions space      *
  *               (if R0 <> N)                             *
  *          VX:  contains in output the M0 solution       *
  *               vectors (stored in columns 1..M0) when   *
  *               R0<>N, or the unique solution vector     *
  *               (stored in column 1) when R0=N.          * 
  *********************************************************/   
  void RSHCGT(double eps, int N, Matc A, int *R0, int *M0,
	            Matc VX)  {
    int i,j,k, k0,l,l0,n0;
    COMPLEX C0,C1, P0,Q0,S, Z0,Z1;
    int I0[NMAX];

    Z0.R=0.0; Z0.I=0.0;
    Z1.R=1.0; Z1.I=0.0;
    for (k=1; k<=N; k++)  I0[k]=k;
    *R0=0; 
	k=1;
    while (*R0==0 && k<=N) { 
      P0.R=A[k][k].R; P0.I=A[k][k].I; 
	  l0=k; k0=k;
      for (i=k; i<=N; i++)
        for (j=k; j<=N; j++)
		  if (CABS(A[i][j]) > CABS(P0)) {
            P0.R=A[i][j].R; P0.I=A[i][j].I; 
			l0=i; k0=j;
          }
	  if (CABS(P0)<eps) {
        *R0=k-1;
        *M0=N-(*R0);
      }
      else
		if (k==N) {
          *R0=N; *M0=0;
        }
        else {
          if (l0!=k)
            for (j=k; j<=N; j++)  CSwap(&A[k][j],&A[l0][j]);
          if (k0!=k)
            for (i=1; i<=N; i++)  CSwap(&A[i][k],&A[i][k0]);
          n0=I0[k]; I0[k]=I0[k0]; I0[k0]=n0;
        }
	  for (i=k+1; i<=N; i++) {
        CDIV(A[i][k],P0,&Q0);
        for (j=k+1; j<=N; j++) {
          C0.R=A[i][j].R; C0.I=A[i][j].I;
          CMUL(Q0,A[k][j],&C1);
          CDIF(C0,C1,&A[i][j]);
        }
      }
      k++;
    }

    if (*R0==N)
	  for (i=1; i<=N; i++) { 
	    VX[i][1].R=Z0.R;
        VX[i][1].I=Z0.I;
      }
	else {
	  for (l=1; l<=*M0; l++) {
        for (i=N; i>=*R0+1; i--)
		  if (i==*R0+l) {
		    VX[i][l].R=Z1.R;
            VX[i][l].I=Z1.I;
          }
          else {
		    VX[i][l].R=Z0.R;
            VX[i][l].I=Z0.I;
          }
		for (i=*R0; i>0; i--) {
          S.R=Z0.R; S.I=Z0.I;
          for (j=i+1; j<=N; j++) {
            C0.R=S.R; C0.I=S.I;
            CMUL(A[i][j],VX[j][l],&C1);
            CADD(C0,C1,&S);
          }
          CDIV(S,A[i][i],&C0);
          VX[i][l].R=-C0.R;
          VX[i][l].I=-C0.I;
        }
      }
      for (i=1; i<N; i++) {
        j=i+1;
        while (j<=N)
		  if (I0[j]==i) {
            for (l=1; l<=*M0; l++)  CSwap(&VX[j][l],&VX[i][l]);
            I0[j]=I0[i]; I0[i]=i;
            j=N+1;
          }
          else
            j++;
      }
    }

  } // RSHCGT()                          
  

void main() {

  N=3;   // size of complex matrix A

  A[1][1].R=-1.0; A[1][1].I= 0.0;
  A[1][2].R= 0.0; A[1][2].I= 1.0;
  A[1][3].R= 0.0; A[1][3].I=-1.0;

  A[2][1].R= 0.0; A[2][1].I=-1.0;
  A[2][2].R=-1.0; A[2][2].I= 0.0;
  A[2][3].R= 1.0; A[2][3].I= 0.0;

  A[3][1].R= 0.0; A[3][1].I= 1.0;
  A[3][2].R= 1.0; A[3][2].I= 0.0;
  A[3][3].R=-1.0; A[3][3].I= 0.0;

  eps=1e-10;

  RSHCGT(eps,N,A,&R0,&M0,VX);


  if (R0==N) {
    printf("\n System solution:\n");
    for (i=1; i<=N; i++)  printf("  X%d = %f\n", i, VX[i][1].R);
  }
  else {
    printf("\n Basic Family of system solutions:\n");
    for (l= 1; l<=M0; l++) {
      printf("  Solution %d\n", l);
      for (i=1; i<=N; i++) {
        printf("   X%d = %f", i, VX[i][l].R);
        if (VX[i][l].I >= 0.0) printf(" + ");
        else printf(" - ");
        printf("%f I\n", fabs(VX[i][l].I));
      }
    }
  }       
  printf("\n");
}

// end of file rshcgt.cpp