/************************************************************
*         Determinant of a complex square matrix            *
*          By Gauss Method with full pivoting               *
* --------------------------------------------------------- *
* SAMPLE RUN:                                               *
* Calculate the determinant of complex matrix:              *
*  ( 47,-15) ( 62,5) (  0,-72) (61, 20)                     *
*  (  6, 14) (-17,3) (-102,91) ( 7,-12)                     *
*  ( 13, 55) ( 32,8) (  41, 7) (25,  1)                     *
*  (111,25)  ( 40,0) ( 12,-82) (58,-30)                     *
*                                                           *
*  Det =  17416564.000000 - 10598320.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 N;
  double eps;
  Matc A;
  COMPLEX det;


  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<1e-12)
      printf(" Complex Divide by zero!\n");
    else {
      C.R=Z2.R; C.I=-Z2.I;
      CMUL(Z1,C,Z);
      Z->R=Z->R/d; Z->I=Z->I/d;
    }
  }
 
  /****************************************************************
  * TSCGT procedure implements the triangularization algorithm of *
  * Gauss with full pivoting at each step for a complex matrix, A *
  * and saves the made transformations in KP and LP.              *
  * ------------------------------------------------------------- *
  * INPUTS:                                                       *
  *          N:   size of complex matrix A                        *
  *          A:   complex matrix of size N x N                    *
  * OUTPUTS;                                                      *
  *          it:  =0 if A is singular, else =1.                   *
  *           C:  contains the upper triangular matrix and the    *
  *               multipliers used during the process.            *
  *          KP:  contains the column exchanges.                  *
  *          LP:  contains the line exchanges.                    *
  ****************************************************************/
  void TSCGT(double eps, int N, Matc A, int *it, Matc C, int *KP, int *LP) {
  int i,j,k,k0,l0;
  COMPLEX C0,C1,P0,T0;

    for (i=1; i<=N; i++)
	  for(j=1; j<=N; j++) {
        C[i][j].R=A[i][j].R;
        C[i][j].I=A[i][j].I;
      }

    *it=1; k=1;
    while (*it==1 && k<=N) {
      P0.R=C[k][k].R; P0.I=C[k][k].I;   
      l0=k; k0=k;
      for (i=k; i<=N; i++)
        for (j=1; j<=N; j++)
		  if (CABS(C[i][j]) > CABS(P0)) {
            P0.R=C[i][j].R; P0.I=C[i][j].I;
            l0=i; k0=j;
          }
      LP[k]=l0; KP[k]=k0;
      if (CABS(P0) < eps)
        *it=0;
      else {
        if (l0!=k)
		  for (j=k; j<=N; j++) {
            T0.R=C[k][j].R; T0.I=C[k][j].I; 
            C[k][j].R=C[l0][j].R;
            C[k][j].I=C[l0][j].I;
            C[l0][j].R=T0.R;
            C[l0][j].I=T0.I;
          }
        if (k0!=k)
		  for (i=1; i<=N; i++) {
            T0.R=C[i][k].R; T0.I=C[i][k].I;
            C[i][k].R=C[i][k0].R;
            C[i][k].I=C[i][k0].I;
            C[i][k0].R=T0.R;
            C[i][k0].I=T0.I;
          }
		for (i=k+1; i<=N; i++) {
          C0.R=C[i][k].R; C0.I=C[i][k].I; 
		  CDIV(C0,P0,&C[i][k]);
          for (j=k+1; j<=N; j++) {
            C0.R=C[i][j].R; C0.I=C[i][j].I;
            CMUL(C[i][k],C[k][j],&C1);
            CDIF(C0,C1,&C[i][j]);
          }
        }
        k++;
      }
    }
    if (*it==1 && CABS(C[N][N]) < eps)  *it=0;
  } //TSCGT()


  /***************************************************************
  * The DCGT procedure calculates the complex determinant of a   *
  * complex square matrix by the Gauss method with full pivoting *
  * ------------------------------------------------------------ *
  * INPUTS:                                                      *
  *         eps:  required precision                             *
  *          N :  size of matrix A                               *
  *          A :  complex matrix of size N x N                   *
  * OUTPUT:                                                      *
  *         det:  complex determinant.                           *                  
  ***************************************************************/
  void DCGT(double eps, int N, Matc A, COMPLEX *det) {
    int it,k,l;
    COMPLEX C0,Z0,Z1;
    Matc C;
    int KP[NMAX], LP[NMAX];
    Z0.R=0.0; Z0.I=0.0;
    Z1.R=1.0; Z1.I=0.0;
    TSCGT(eps,N,A,&it,C,KP,LP);

    if (it==0) {
      det->R=Z0.R; det->I=Z0.I;
    }
    else {
      det->R=Z1.R; det->I=Z1.I;
      for (k=1; k<=N; k++) {
        C0.R=det->R; C0.I=det->I;  
		CMUL(C0,C[k][k],det);
      }
      l=0;
      for (k=1; k<N; k++) {
        if (LP[k]!=k) l++;
        if (KP[k]!=k) l++;
      }
	  if ((l%2)!=0) { 
        det->R=-det->R; det->I=-det->I;
      }
    }
  } // DCGT()


void main() { 


/* Example #1
  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;

  ( Det = -4.0 + 0 I )

  Example #2  */

  N=4;

  A[1][1].R=47.0; A[1][1].I=-15.0;
  A[1][2].R=62.0; A[1][2].I=  5.0;
  A[1][3].R= 0.0; A[1][3].I=-72.0;
  A[1][4].R=61.0; A[1][4].I= 20.0;

  A[2][1].R=   6.0; A[2][1].I= 14.0;
  A[2][2].R=- 17.0; A[2][2].I=  3.0;
  A[2][3].R=-102.0; A[2][3].I= 91.0;
  A[2][4].R=   7.0; A[2][4].I=-12.0;

  A[3][1].R= 13.0; A[3][1].I=-55.0;
  A[3][2].R= 32.0; A[3][2].I=  8.0;
  A[3][3].R= 41.0; A[3][3].I=  7.0;
  A[3][4].R= 25.0; A[3][4].I=  1.0;

  A[4][1].R=111.0; A[4][1].I= 25.0;
  A[4][2].R= 40.0; A[4][2].I=  0.0;
  A[4][3].R= 12.0; A[4][3].I=-82.0;
  A[4][4].R= 58.0; A[4][4].I=-30.0;

  eps=1e-10;

  DCGT(eps, N, A, &det);

  printf("\n  Det = %f", det.R);
  if (det.I>=0.0) printf(" + ");
  else printf(" - ");
  printf("%f I\n\n", fabs(det.I));

}

// end of file cdetmat.cpp