DECLARE FUNCTION CABS# (X#, Y#)
DECLARE SUB CDIV (Z1#(), Z2#(), Z#())
DECLARE SUB CMUL (Z1#(), Z2#(), Z#())
'**************************************************************************************
'*             Inverse of a complex square matrix A by LU decomposition               *
'*                                                                                    *
'*                                                Basic version by J-P Moreau, Paris  *
'*                                                        (www.jpmoreau.fr)           *
'* ---------------------------------------------------------------------------------- *
'* Reference:                                                                         *
'*                                                                                    *
'*  "Numerical Recipes by W.H. Press, B. P. Flannery, S.A. Teukolsky                  *
'*   and W.T. Vetterling, Cambridge University Press, 1986."                          *
'* ---------------------------------------------------------------------------------- *
'* SAMPLE RUN:                                                                        *
'*                                                                                    *
'* Input file (test_clu.dat):                                                         *
'*                                                                                    *
'* 4                                                                                  *
'*  8  0  2    0  3  0  12     1                                                      *
'*  2  0  4    0  7  0   0.25  2                                                      *
'*  3  0  7    0  3  0   5     3                                                      *
'* 12  0  0.25 0  5  0   2     4                                                      *
'*                                                                                    *
'* Output file (inv_clu.lst):                                                         *
'* ----------------------------------------------------                               *
'*  INVERSION OF A COMPLEX MATRIX BY LU DECOMPOSITION:                                *
'* ----------------------------------------------------                               *
'*                                                                                    *
'*  N= 4                                                                              *
'*                                                                                    *
'*  8.00   0.00    2.00   0.00    3.00   0.00   12.00   1.00                          *
'*  2.00   0.00    4.00   0.00    7.00   0.00    0.25   2.00                          *
'*  3.00   0.00    7.00   0.00    3.00   0.00    5.00   3.00                          *
'* 12.00   0.00    0.25   0.00    5.00   0.00    2.00   4.00                          *
'*                                                                                    *
'*  Inverted Matrix:                                                                  *
'*                                                                                    *
'* -0.03022  -0.04162   -0.08426  -0.00641    0.05306   0.01466    0.10426   0.02515  *
'* -0.07421  -0.04847   -0.06024  -0.00746    0.19813   0.01707    0.00998   0.02929  *
'*  0.05443   0.00890    0.20187   0.00137   -0.12957  -0.00313   -0.03754  -0.00538  *
'*  0.10431   0.02491    0.01606   0.00384   -0.03673  -0.00877   -0.06304  -0.01505  *
'*                                                                                    *
'*  Verification A*Y = I:                                                             *
'*                                                                                    *
'*  1.00000   0.00000   -0.00000  -0.00000    0.00000   0.00000   -0.00000   0.00000  *
'*  0.00000  -0.00000    1.00000   0.00000    0.00000   0.00000   -0.00000   0.00000  *
'* -0.00000  -0.00000   -0.00000  -0.00000    1.00000   0.00000    0.00000   0.00000  *
'*  0.00000   0.00000    0.00000   0.00000   -0.00000   0.00000    1.00000   0.00000  *
'* ----------------------------------------------------                               *
'*  End of file inv_clu.lst                                                           *
'*                                                                                    *
'* ---------------------------------------------------------------------------------- *
'* Subroutines used:                                                                  *
'*   1000 LUDCMP: LU decomposition of a complex square matrix                         *
'*   2000 LUBKSB: Solve linear system AX = B, A being the output of LUDCMP.           *
'*                or inversion of A complex matrix.                                   *
'*                                                                                    *
'**************************************************************************************
 DEFINT I-N
 DEFDBL A-H, O-Z

  CLS : PRINT
  INPUT " Data file name (without .dat): ", nom$
  in$ = nom$ + ".dat"
  ou$ = nom$ + ".lst"

  OPEN in$ FOR INPUT AS #1
  OPEN ou$ FOR OUTPUT AS #2

  INPUT #1, N     'size of complex matrix

  n1 = N + 1
  N2 = n1 * n1
  DIM AR(N2), AI(N2), BR(n1), BI(n1), temp(n1, 2), INDX(n1)
  DIM A1R(N2), A1I(N2), YR(N2), YI(N2)
  DIM SUM(2), TMP(2), TMP1(2), TMP2(2)

'------------------------------------------------------------------
'NOTA: Complex Matrix A of size N+1 x N+1 is stored in txo vectors,
'      line by. AR one for real part and AI for imaginary part.
'      The element (i,j) of the matrix is at location i(n+1) + j in
'      the corresponding vector. The access is faster.
'------------------------------------------------------------------

  PRINT #2, "----------------------------------------------------"
  PRINT #2, " INVERSION OF A COMPLEX MATRIX BY LU DECOMPOSITION: "
  PRINT #2, "----------------------------------------------------"
  PRINT #2,
  PRINT #2, "  N="; N
  PRINT #2,

  FOR I = 1 TO N
    FOR J = 1 TO N
      INPUT #1, temp(J, 1), temp(J, 2)  'read a line
    NEXT J
    FOR J = 1 TO N
      AR(I * n1 + J) = temp(J, 1)
      AI(I * n1 + J) = temp(J, 2)
      A1R(I * n1 + J) = temp(J, 1)  'copy of A for later verfication
      A1I(I * n1 + J) = temp(J, 2)
      YR(I * n1 + J) = 0#
      YI(I * n1 + J) = 0#
    NEXT J
    FOR J = 1 TO N
      PRINT #2, USING "  ###.## ###.##"; temp(J, 1); temp(J, 2);
    NEXT J
    PRINT #2,
    YR(I * n1 + I) = 1#
  NEXT I
  CLOSE #1

  GOSUB 1000  'LUDCMP(A,n,INDX,D,CODE)

'call solver if previous return code is ok
'to obtain inverse of A one column at a time
  IF CODE = 0 THEN
    FOR JJ = 1 TO N   'Caution, J is used in 2000
      FOR I = 1 TO N
	BR(I) = YR(I * n1 + JJ)
	BI(I) = YI(I * n1 + JJ)
      NEXT I
      GOSUB 2000  'LUBKSB(A,n,INDX,B);
      FOR I = 1 TO N
	YR(I * n1 + JJ) = BR(I)
	YI(I * n1 + JJ) = BI(I)
      NEXT I
    NEXT JJ
  END IF
' inverse of matrix A is now in matrix Y
' the original matrix A is destroyed

  IF CODE = 1 THEN
    PRINT #2, " The complex matrix is singular!"
  ELSE
    PRINT #2,
    PRINT #2, "  Inverted Matrix:"
    PRINT #2,
    FOR I = 1 TO N
      FOR J = 1 TO N
	PRINT #2, USING "  ###.##### ###.#####"; YR(I * n1 + J); YI(I * n1 + J);
      NEXT J
      PRINT #2,
    NEXT I
  END IF

' verify A1 x Y = I (result put in A)
' Multiply A1 by Y and put in A
  FOR I = 1 TO N
    FOR J = 1 TO N
      SUM(1) = 0#: SUM(2) = 0#
      FOR K = 1 TO N
	TMP(1) = A1R(I * n1 + K): TMP(2) = A1I(I * n1 + K)
	TMP1(1) = YR(K * n1 + J): TMP1(2) = YI(K * n1 + J)
	CMUL TMP(), TMP1(), TMP2()
	SUM(1) = SUM(1) + TMP2(1)
	SUM(2) = SUM(2) + TMP2(2)
	AR(I * n1 + J) = SUM(1)
	AI(I * n1 + J) = SUM(2)
      NEXT K
    NEXT J
  NEXT I

' A should now contain identity matrix
  PRINT #2,
  PRINT #2, "  Verification A*Y = I:"
  PRINT #2,
  FOR I = 1 TO N
    FOR J = 1 TO N
      PRINT #2, USING "  ###.##### ###.#####"; AR(I * n1 + J); AI(I * n1 + J);
    NEXT J
    PRINT #2,
  NEXT I

  PRINT #2, "----------------------------------------------------"
  PRINT #2, " End of file "; ou$
  CLOSE #2
  PRINT
  PRINT " Results in file "; ou$
  PRINT

END 'of main program


'***************************************************************
'* 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.   *
'***************************************************************
1000 'SUBROUTINE LUDCMP FOR COMPLEX MATRIX
  NMAX = 100: TINY = 2.2E-16
  DIM VV(N, 2), AMAX(2), DUM(2)

  D = 1: CODE = 0

  FOR I = 1 TO N
    AMAX(1) = 0#: AMAX(2) = 0#
    FOR J = 1 TO N
      IF CABS(AR(I * n1 + J), AI(I * n1 + J)) > CABS(AMAX(1), AMAX(2)) THEN
	AMAX(1) = AR(I * n1 + J)
	AMAX(2) = AI(I * n1 + J)
      END IF
    NEXT J

    IF CABS(AMAX(1), AMAX(2)) < TINY THEN
      CODE = 1
      RETURN
    END IF

'   VV(I) = 1.0 / AMAX
    TMP(1) = 1#: TMP(2) = 0#
    CDIV TMP(), AMAX(), DUM()
    VV(I, 1) = DUM(1): VV(I, 2) = DUM(2)
  NEXT I  'i loop

  FOR J = 1 TO N
    FOR I = 1 TO J - 1
      SUM(1) = AR(I * n1 + J)
      SUM(2) = AI(I * n1 + J)
      FOR K = 1 TO I - 1
'       SUM = SUM - A(I * N1 + K) * A(K * N1 + J)
	TMP(1) = AR(I * n1 + K): TMP(2) = AI(I * n1 + K)
	TMP1(1) = AR(K * n1 + J): TMP1(2) = AI(K * n1 + J)
	CMUL TMP(), TMP1(), DUM()
	SUM(1) = SUM(1) - DUM(1)
	SUM(2) = SUM(2) - DUM(2)
      NEXT K
      AR(I * n1 + J) = SUM(1)
      AI(I * n1 + J) = SUM(2)
    NEXT I  'i loop
    AMAX(1) = 0#: AMAX(2) = 0#
    FOR I = J TO N
      SUM(1) = AR(I * n1 + J)
      SUM(2) = AI(I * n1 + J)
      FOR K = 1 TO J - 1
'       SUM = SUM - A(I * N1 + K) * A(K * N1 + J)
	TMP(1) = AR(I * n1 + K): TMP(2) = AI(I * n1 + K)
	TMP1(1) = AR(K * n1 + J): TMP1(2) = AI(K * n1 + J)
	CMUL TMP(), TMP1(), DUM()
	SUM(1) = SUM(1) - DUM(1)
	SUM(2) = SUM(2) - DUM(2)
      NEXT K
      AR(I * n1 + J) = SUM(1)
      AI(I * n1 + J) = SUM(2)
'     DUM = VV(I) * ABS(SUM)
      DUM(1) = VV(I, 1) * CABS(SUM(1), SUM(2))
      DUM(2) = VV(I, 2) * CABS(SUM(1), SUM(2))
      IF CABS(DUM(1), DUM(2)) >= CABS(AMAX(1), AMAX(2)) THEN
	IMAX = I
	AMAX(1) = DUM(1): AMAX(2) = DUM(2)
      END IF
    NEXT I  'i loop

    IF J <> IMAX THEN
      FOR K = 1 TO N
	DUM(1) = AR(IMAX * n1 + K)
	DUM(2) = AI(IMAX * n1 + K)
	AR(IMAX * n1 + K) = AR(J * n1 + K)
	AI(IMAX * n1 + K) = AI(J * n1 + K)
	AR(J * n1 + K) = DUM(1)
	AI(J * n1 + K) = DUM(2)
      NEXT K  'k loop
      D = -D
      VV(IMAX, 1) = VV(J, 1)
      VV(IMAX, 2) = VV(J, 2)
    END IF

    INDX(J) = IMAX
    IF CABS(AR(J * n1 + J), AI(J * n1 + J)) < TINY THEN
      AR(J * n1 + J) = TINY
      AI(J * n1 + J) = 0#
    END IF

    IF J <> N THEN
'     DUM = 1# / A(J * N1 + J)
      TMP(1) = 1#: TMP(2) = 0#
      TMP1(1) = AR(J * n1 + J): TMP1(2) = AI(J * n1 + J)
      CDIV TMP(), TMP1(), DUM()
      FOR I = J + 1 TO N
'       A(I * N1 + J) = A(I * N1 + J) * DUM
	TMP(1) = AR(I * n1 + J): TMP(2) = AI(I * n1 + J)
	CMUL TMP(), DUM(), TMP1()
	AR(I * n1 + J) = TMP1(1)
	AI(I * n1 + J) = TMP1(2)
      NEXT I
    END IF
  NEXT J  'j loop

RETURN  '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.                     *
'******************************************************************
2000  'SUBROUTINE LUBKSB
  II = 0: n1 = N + 1
  FOR I = 1 TO N
    LL = INDX(I)
    SUM(1) = BR(LL): SUM(2) = BI(LL)
    BR(LL) = BR(I): BI(LL) = BI(I)
    IF II <> 0 THEN
      FOR J = II TO I - 1
'       SUM = SUM - A(I * N1 + J) * B(J)
	TMP(1) = AR(I * n1 + J): TMP(2) = AI(I * n1 + J)
	TMP1(1) = BR(J): TMP1(2) = BI(J)
	CMUL TMP(), TMP1(), DUM()
	SUM(1) = SUM(1) - DUM(1)
	SUM(2) = SUM(2) - DUM(2)
      NEXT J
    ELSEIF CABS(SUM(1), SUM(2)) <> 0 THEN
      II = I
    END IF
    BR(I) = SUM(1): BI(I) = SUM(2)
  NEXT I  'i loop

  FOR I = N TO 1 STEP -1
    SUM(1) = BR(I): SUM(2) = BI(I)
    IF I < N THEN
      FOR J = I + 1 TO N
'       SUM = SUM - A(I * N1 + J) * B(J)
	TMP(1) = AR(I * n1 + J): TMP(2) = AI(I * n1 + J)
	TMP1(1) = BR(J): TMP1(2) = BI(J)
	CMUL TMP(), TMP1(), DUM()
	SUM(1) = SUM(1) - DUM(1)
	SUM(2) = SUM(2) - DUM(2)
      NEXT J
    END IF
'   B(I) = SUM / A(I * N1 + I)
    TMP(1) = AR(I * n1 + I): TMP(2) = AI(I * n1 + I)
    CDIV SUM(), TMP(), DUM()
    BR(I) = DUM(1): BI(I) = DUM(2)
  NEXT I  'i loop
RETURN

  FUNCTION CABS (X, Y)
    CABS = SQR(X * X + Y * Y)
  END FUNCTION

  SUB CDIV (Z1(), Z2(), Z())
    DIM C1(2)
    D = Z2(1) ^ 2 + Z2(2) ^ 2
    IF D < 2.2E-16 THEN
      PRINT #2, " Complex Divide by zero!"
    ELSE
      C1(1) = Z2(1): C1(2) = -Z2(2)
      CMUL Z1(), C1(), Z()
      Z(1) = Z(1) / D: Z(2) = Z(2) / D
    END IF
  END SUB

  SUB CMUL (Z1(), Z2(), Z())
    Z(1) = Z1(1) * Z2(1) - Z1(2) * Z2(2)
    Z(2) = Z1(1) * Z2(2) + Z1(2) * Z2(1)
  END SUB