/******************************************************
*    LU decomposition routines used by test_lu.cpp    *
*    with static allocations                          *
*                                                     *
*                  C++ version by J-P Moreau, Paris   *
*                  for a complex linear system.       *
*                       (www.jpmoreau.fr)             * 
* --------------------------------------------------- *
* Reference:                                          *
*                                                     *
* "Numerical Recipes by W.H. Press, B. P. Flannery,   *
*  S.A. Teukolsky and W.T. Vetterling, Cambridge      *
*  University Press, 1986".                           *
*                                                     *
******************************************************/
#include <stdio.h>
#include <math.h>

#define SIZE 25
#define TINY 2.5e-16

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

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


  double CABS(COMPLEX Z) {
    return sqrt(Z.R*Z.R + Z.I*Z.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;
    }
  }


/**************************************************************
* Given an N x N matrix A, this routine replaces it by the LU *
* decomposition of a rowwise permutation of itself. A and N   *
* are input. INDX is an output vector which records the row   *
* permutation effected by the partial pivoting; D is output   *
* as -1 or 1, depending on whether the number of row inter-   *
* changes was even or odd, respectively. This routine is used *
* in combination with LUBKSB to solve linear equations or to  *
* invert a matrix. Return code is 1, if matrix is singular.   *
**************************************************************/
int LUDCMP(CMAT A, int n, IMAT INDX, int *d)  {

  COMPLEX AMAX,DUM, SUM, TMP;
  int  I,IMAX,J,K;
  CVEC VV;

  *d=1; 

  for  (I=1; I<n+1; I++)  {
    AMAX.R=0.0; AMAX.I=0.0; 
    for (J=1; J<n+1; J++)  
	  if (CABS(A[I][J]) > CABS(AMAX)) {
	    AMAX.R=A[I][J].R;
        AMAX.I=A[I][J].I;
      } 


    if(CABS(AMAX) < TINY)  return 1;
 // VV[I] = 1.0 / AMAX;
    TMP.R=1.0; TMP.I=0.0;
	CDIV(TMP, AMAX, &VV[I]);
  } // i loop 

  for (J=1; J<n+1;J++)  {

    for (I=1; I<J; I++)  { 
      SUM.R = A[I][J].R;
      SUM.I = A[I][J].I;
	  for (K=1; K<I; K++) {  
//      SUM = SUM - A[I][K]*A[K][J]
        CMUL(A[I][K], A[K][J], &TMP);
		SUM.R = SUM.R - TMP.R; 
		SUM.I = SUM.I - TMP.I;
      }
      A[I][J].R = SUM.R;
      A[I][J].I = SUM.I;
	} // i loop 
    
	AMAX.R = 0.0; AMAX.I=0.0;

	for (I=J; I<n+1; I++)  {
      SUM = A[I][J];
      for  (K=1; K<J; K++) {  
//      SUM = SUM - A[I][K]*A[K][J]
        CMUL(A[I][K], A[K][J], &TMP);
		SUM.R = SUM.R - TMP.R; 
		SUM.I = SUM.I - TMP.I;
      }
      A[I][J].R = SUM.R;
      A[I][J].I = SUM.I;
//    DUM = VV[I]*ABS(SUM)
      DUM.R=VV[I].R*CABS(SUM);
      DUM.I=VV[I].I*CABS(SUM);

      if (CABS(DUM) >= CABS(AMAX)) {
        IMAX = I;
        AMAX.R = DUM.R; AMAX.I = DUM.I; 
      }
	} // i loop
   
	if (J != IMAX)  {
	  for (K=1; K<n+1; K++)  {
        DUM.R = A[IMAX][K].R;
        DUM.I = A[IMAX][K].I; 
        A[IMAX][K].R = A[J][K].R;
        A[IMAX][K].I = A[J][K].I;
        A[J][K].R = DUM.R;
        A[J][K].I = DUM.I;
	  } // k loop 
      *d = -*d;
      VV[IMAX].R = VV[J].R;
      VV[IMAX].I = VV[J].I;
    }

    INDX[J] = IMAX;

    if (CABS(A[J][J]) < TINY) {
	  A[J][J].R = TINY;
      A[J][J].I = 0.0;
    }

    if (J != n)  {
//    DUM = 1.0 / A[J][J]
      TMP.R=1.0; TMP.I=0.0;
	  CDIV(TMP, A[J][J], &DUM);
      for (I=J+1; I<n+1; I++) {  
//      A[I][J] *= DUM;
        CMUL(A[I][J], DUM, &TMP);
        A[I][J].R = TMP.R;
        A[I][J].I = TMP.I;
      }
    } 
  } // j loop 

  return 0;

} // subroutine LUDCMP 


/*****************************************************************
 * Solves the set of N linear equations A . X = B.  Here A is     *
 * input, not as the matrix A but rather as its LU decomposition, *
 * determined by the routine LUDCMP. INDX is input as the permuta-*
 * tion vector returned by LUDCMP. B is input as the right-hand   *
 * side vector B, and returns with the solution vector X. A, N and*
 * INDX are not modified by this routine and can be used for suc- *
 * cessive calls with different right-hand sides. This routine is *
 * also efficient for plain matrix inversion.                     *
 *****************************************************************/
void LUBKSB(CMAT A, int n, IMAT INDX, CVEC B)  {
  COMPLEX SUM, TMP;
  int  I,II,J,LL;

  II = 0; 

  for (I=1; I<n+1; I++)  {
    LL = INDX[I];
    SUM.R = B[LL].R; SUM.I = B[LL].I;
    B[LL].R = B[I].R; B[LL].I = B[I].I;
    if (II != 0) 
	  for (J=II; J<I; J++) {  
//      SUM = SUM - A[I][J]*B[J]
        CMUL(A[I][J],B[J],&TMP);
		SUM.R -= TMP.R;
		SUM.I -= TMP.I;
      }
    else if (CABS(SUM) != 0.0)  II = I;
    B[I].R = SUM.R; B[I].I = SUM.I;
  } // i loop

  for (I=n; I>0; I--)  {
    SUM.R = B[I].R; SUM.I = B[I].I;
    if (I < n)  {
	  for (J=I+1; J<n+1; J++) {
//      SUM = SUM - A[I][J]*B[J]
        CMUL(A[I][J],B[J],&TMP);
		SUM.R -= TMP.R;
		SUM.I -= TMP.I;
      }
    }
//  B[I] = SUM / A[I][I]
	CDIV(SUM, A[I][I], &B[I]);
  } // i loop 


} // LUBKSB    

// end of file clu1.cpp