DECLARE SUB ICGT (eps!, N%, it%)
DECLARE SUB MATMUL (N%)
DECLARE FUNCTION CABS! (Z!())
DECLARE SUB CDIV (Z1!(), Z2!(), Z!())
DECLARE SUB CMUL (Z1!(), Z2!(), Z!())
DECLARE SUB TSCGT (eps!, N%, it%, KP%(), LP%())
DECLARE SUB BSCGT (N%, KP%(), LP%())
'****************************************************************
'* calculate inverse of a complex square matrix by Gauss Method *
'* with full pivoting.                                          *
'* ------------------------------------------------------------ *
'* SAMPLE RUN:                                                  *
'* Calculate inverse of complex square matrix:                  *
'*                                                              *
'*       | (0, 1) ( 2,0) (-1,-1) |                              *
'*   A = | (1, 0) ( 1,1) ( 0,-3) |                              *
'*       | (0,-2) (-1,1) ( 3, 0) |                              *
'*                                                              *
'* Inverse matrix:                                              *
'* (  0.0000,  0.0000) (  0.3333,  0.0000) (  0.0000,  0.3333)  *
'* (  0.7500,  0.0000) ( -0.1667, -0.0833) (  0.3333,  0.0833)  *
'* (  0.2500, -0.2500) ( -0.0833,  0.2500) (  0.2500, -0.0833)  *
'*                                                              *
'* Product AM1*A:                                               *
'* (  1.0000,  0.0000) (  0.0000,  0.0000) (  0.0000,  0.0000)  *
'* (  0.0000,  0.0000) (  1.0000,  0.0000) (  0.0000,  0.0000)  *
'* (  0.0000,  0.0000) (  0.0000,  0.0000) (  1.0000,  0.0000)  *
'*                                                              *
'* ------------------------------------------------------------ *
'* Ref.: "Algebre, Algorithmes et programmes en Pascal          *
'*        By Jean-Louis Jardrin, DUNOD Paris, 1988".            *
'*                                                              *
'*                   Quick Basic Release By J-P Moreau, Paris.  *
'*                             (www.jpmoreau.fr)                *
'**************************************************************** 
'Program InvMatC

DEFINT I-N

OPTION BASE 1

DIM SHARED NMAX

  NMAX = 20

DIM SHARED A(NMAX, NMAX, 2), AM1(NMAX, NMAX, 2), MI(NMAX, NMAX, 2)

DIM SHARED C(NMAX, NMAX, 2), W(NMAX, 2), Z(NMAX, 2)


  N = 3  'size of complex matrix A

  A(1, 1, 1) = 0!: A(1, 1, 2) = 1!
  A(1, 2, 1) = 2!: A(1, 2, 2) = 0!
  A(1, 3, 1) = -1!: A(1, 3, 2) = -1!

  A(2, 1, 1) = 1!: A(2, 1, 2) = 0!
  A(2, 2, 1) = 1!: A(2, 2, 2) = 1!
  A(2, 3, 1) = 0!: A(2, 3, 2) = -3!

  A(3, 1, 1) = 0!: A(3, 1, 2) = -2!
  A(3, 2, 1) = -1!: A(3, 2, 2) = 1!
  A(3, 3, 1) = 3!: A(3, 3, 2) = 0!

  eps = 1E-10

  ICGT eps, N, it

  F$ = " (###.####,###.####)"

  CLS
  PRINT
  PRINT " Inverse matrix:"
  FOR i = 1 TO N
    FOR j = 1 TO N
      PRINT USING F$; AM1(i, j, 1); AM1(i, j, 2);
    NEXT j
    PRINT
  NEXT i

' Check AM1 * A = I
  MATMUL N

  PRINT
  PRINT " Product AM1*A:"
  FOR i = 1 TO N
    FOR j = 1 TO N
      PRINT USING F$; MI(i, j, 1); MI(i, j, 2);
    NEXT j
    PRINT
  NEXT i
  PRINT

END 'of main program

 '**************************************************************
 '* Function BSCGT calculates the solution of upper triangular *
 '* system [S(n-1)] by the back substitution method and deduces*
 '* from it the solution of system [S]. The call must be made  *
 '* only after a call to TSCGT and the matrix of [S] must be   *
 '* regular.                                                   *
 '* ---------------------------------------------------------- *
 '* INPUTS:                                                    *
 '*         N : size of matrix C                               *
 '*         C:  contains the upper triangular matrix and the   *
 '*             multipliers used during the triangularization  *
 '*             process (in output of TSCGT).     .            *
 '*         W : contains the right-side coefficients           *
 '*         KP: contains the column exchanges.                 *
 '*         LP: contains the line exchanges.                   *
 '* OUTPUT:                                                    *
 '*         Z : system solution complex vector.                *
 '**************************************************************
 'Note: C, W, Z are shared with main program and ICGT
  SUB BSCGT (N, KP(), LP())
  DIM C0(2), C1(2), S(2), Z0(2)'Complex numbers
  DIM TMP(2), TMP1(2)

    Z0(1) = 0!: Z0(2) = 0!
    FOR k = 1 TO N - 1
      l0 = LP(k)
      IF l0 <> k THEN
        'CSwap(W(k),W(l0))
        C0(1) = W(l0, 1): C0(2) = W(l0, 2)
        W(l0, 1) = W(k, 1): W(l0, 2) = W(k, 2)
        W(k, 1) = C0(1): W(k, 2) = C0(2)
      END IF
      FOR i = k + 1 TO N
        C0(1) = W(i, 1): C0(2) = W(i, 2)
        'CMUL(C(i,k),W(k),C1)
        TMP(1) = C(i, k, 1): TMP(2) = C(i, k, 2)
        TMP1(1) = W(k, 1): TMP1(2) = W(k, 2)
        CMUL TMP(), TMP1(), C1()
        'CDIF(C0,C1,W(i))
        W(i, 1) = C0(1) - C1(1)
        W(i, 2) = C0(2) - C1(2)
      NEXT i
    NEXT k
    'CDIV(W(N),C(N,N),Z(N))
    TMP(1) = W(N, 1): TMP(2) = W(N, 2)
    TMP1(1) = C(N, N, 1): TMP1(2) = C(N, N, 2)
    CDIV TMP(), TMP1(), C1()
    Z(N, 1) = C1(1): Z(N, 2) = C1(2)
    FOR i = N - 1 TO 1 STEP -1
      S(1) = Z0(1): S(2) = Z0(2)
      FOR j = i + 1 TO N
        C0(1) = S(1): C0(2) = S(2)
        'CMUL(C(i,j),Z(j),C1)
        TMP(1) = C(i, j, 1): TMP(2) = C(i, j, 2)
        TMP1(1) = Z(j, 1): TMP1(2) = Z(j, 2)
        CMUL TMP(), TMP1(), C1()
        'CADD(C0,C1,S)
        S(1) = C0(1) + C1(1): S(2) = C0(2) + C1(2)
      NEXT j
      'CDIF(W(i),S,C0)
      C0(1) = W(i, 1) - S(1): C0(2) = W(i, 2) - S(2)
      'CDIV(C0,C(i,i),Z(i))
      TMP(1) = C(i, i, 1): TMP(2) = C(i, i, 2)
      CDIV C0(), TMP(), TMP1()
      Z(i, 1) = TMP1(1): Z(i, 2) = TMP1(2)
    NEXT i
    FOR k = N - 1 TO 1 STEP -1
      k0 = KP(k)
      IF k0 <> k THEN
        'CSwap(Z(k),Z(k0))
        TMP(1) = Z(k0, 1): TMP(2) = Z(k0, 2)
        Z(k0, 1) = Z(k, 1): Z(k0, 2) = Z(k, 2)
        Z(k, 1) = TMP(1): Z(k, 2) = TMP(2)
      END IF
    NEXT k
  END SUB  'BSCGT

  FUNCTION CABS (Z())
    CABS = SQR(Z(1) ^ 2 + Z(2) ^ 2)
  END FUNCTION

  SUB CADD (Z1(), Z2(), Z())
    Z(1) = Z1(1) + Z2(1)
    Z(2) = Z1(2) + Z2(2)
  END SUB

  SUB CDIF (Z1(), Z2(), Z())
    Z(1) = Z1(1) - Z2(1)
    Z(2) = Z1(2) - Z2(2)
  END SUB

  SUB CDIV (Z1(), Z2(), Z())
    DIM C1(2)
    D = Z2(1) ^ 2 + Z2(2) ^ 2
    IF D < 2E-12 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

 '****************************************************************
 '* calculate inverse of a complex square matrix by Gauss Method *
 '* with full pivoting.                                          *
 '* ------------------------------------------------------------ *
 '* INPUTS:                                                      *
 '*        eps : required precision                              *
 '*         N  : size of complex matrix A                        *
 '*         A  : input complex matrix (NxN)                      *
 '* OUTPUTS:                                                     *
 '*         it : =0, matrix is singular; =1, matris is regular   *
 '*        AM1 : inverse complex matrix (NxN).                   *
 '****************************************************************
 ' A, AM1, C, W, Z are shared with main program and BSCGT
  SUB ICGT (eps, N, it)
  DIM Z0(2), Z1(2)
  DIM KP(NMAX), LP(NMAX)

    Z0(1) = 0!: Z0(2) = 0!
    Z1(1) = 1!: Z1(2) = 0!

    TSCGT eps, N, it, KP(), LP()

    IF it = 1 THEN
      FOR l = 1 TO N
        FOR i = 1 TO N
          IF i = l THEN
            W(i, 1) = Z1(1): W(i, 2) = Z1(2)
          ELSE
            W(i, 1) = Z0(1): W(i, 2) = Z0(2)
          END IF
        NEXT i

        BSCGT N, KP(), LP()

        FOR i = 1 TO N
          AM1(i, l, 1) = Z(i, 1)
          AM1(i, l, 2) = Z(i, 2)
        NEXT i
      NEXT l
    END IF
  END SUB 'ICGT    

SUB MATMUL (N)
'*******************************************                                    
'*  MULTIPLICATION OF TWO SQUARE COMPLEX   *                                    
'*  MATRICES                               *
'* --------------------------------------- *                                    
'* INPUTS:  AM1  MATRIX N*N                *                                    
'*            A  MATRIX N*N                *                                    
'*            N  INTEGER                   *                                    
'* --------------------------------------- *                                    
'* OUTPUTS:  MI  MATRIX N*N PRODUCT AM1*A  *                                    
'*                                         *
'*******************************************
  DIM SUM(2), PROD(2), TMP(2), TMP1(2)'Complex numbers
                                         
  FOR i = 1 TO N
    FOR j = 1 TO N
      SUM(1) = 0!: SUM(2) = 0!
      FOR k = 1 TO N
        'SUM=SUM+AM1(I,K)*A(K,J)
        TMP(1) = AM1(i, k, 1): TMP(2) = AM1(i, k, 2)
        TMP1(1) = A(k, j, 1): TMP1(2) = A(k, j, 2)
        CMUL TMP(), TMP1(), PROD()
        'CADD(SUM,PROD,SUM)
        SUM(1) = SUM(1) + PROD(1)
        SUM(2) = SUM(2) + PROD(2)
      NEXT k
      MI(i, j, 1) = SUM(1)
      MI(i, j, 2) = SUM(2)
    NEXT j
  NEXT i
END SUB

 '*****************************************************************
 '* 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.                    *
 '*****************************************************************
  SUB TSCGT (eps, N, it, KP(), LP())
  DIM C0(2), C1(2), P0(2), TMP(2), TMP1(2)  'Complex numbers

    FOR i = 1 TO N
      FOR j = 1 TO N
        C(i, j, 1) = A(i, j, 1)
        C(i, j, 2) = A(i, j, 2)
      NEXT j
    NEXT i

    it = 1: k = 1
    WHILE it = 1 AND k < N
      P0(1) = C(k, k, 1): P0(2) = C(k, k, 2)
      l0 = k: k0 = k
      FOR i = k TO N
        FOR j = k TO N
          C0(1) = C(i, j, 1): C0(2) = C(i, j, 2)
          IF CABS(C0()) > CABS(P0()) THEN
            P0(1) = C(i, j, 1): P0(2) = C(i, j, 2)
            l0 = i: k0 = j
          END IF
        NEXT j
      NEXT i
      LP(k) = l0: KP(k) = k0
      IF CABS(P0()) < eps THEN
        it = 0
      ELSE
        IF l0 <> k THEN
          FOR j = k TO N
            'CSwap(C(k,j),C(l0,j))
            C0(1) = C(l0, j, 1): C0(2) = C(l0, j, 2)
            C(l0, j, 1) = C(k, j, 1): C(l0, j, 2) = C(k, j, 2)
            C(k, j, 1) = C0(1): C(k, j, 2) = C0(2)
          NEXT j
        END IF
        IF k0 <> k THEN
          FOR i = 1 TO N
            'CSwap(C(i,k),C(i,k0))
            C0(1) = C(i, k0, 1): C0(2) = C(i, k0, 2)
            C(i, k0, 1) = C(i, k, 1): C(i, k0, 2) = C(i, k, 2)
            C(i, k, 1) = C0(1): C(i, k, 2) = C0(2)
          NEXT i
        END IF
        FOR i = k + 1 TO N
          C0(1) = C(i, k, 1): C0(2) = C(i, k, 2)
          'CDIV(C0,P0,C(i,k))
          CDIV C0(), P0(), TMP()
          C(i, k, 1) = TMP(1): C(i, k, 2) = TMP(2)
          FOR j = k + 1 TO N
            C0(1) = C(i, j, 1): C0(2) = C(i, j, 2)
            'CMUL(C(i,k),C(k,j),C1)
            TMP(1) = C(i, k, 1): TMP(2) = C(i, k, 2)
            TMP1(1) = C(k, j, 1): TMP1(2) = C(k, j, 2)
            CMUL TMP(), TMP1(), C1()
            'CDIF(C0,C1,C(i,j))
            C(i, j, 1) = C0(1) - C1(1)
            C(i, j, 2) = C0(2) - C1(2)
          NEXT j
        NEXT i
        k = k + 1
      END IF
    WEND
    C0(1) = C(N, N, 1): C0(2) = C(N, N, 2)
    IF it = 1 AND CABS(C0()) < eps THEN it = 0
  END SUB  'TSCGT