DECLARE SUB RSLCGTC (eps!, dta!, M%, N%, it%)
DECLARE FUNCTION CABS! (Z!())
DECLARE SUB CADD (Z1!(), Z2!(), Z!())
DECLARE SUB CDIV (Z1!(), Z2!(), Z!())
DECLARE SUB CMUL (Z1!(), Z2!(), Z!())
DECLARE SUB BSCGT (N%, KP%(), LP%())
DECLARE SUB TSCGT (eps!, N%, it%, KP%(), LP%())
'********************************************************************
'* Solve a Complex Linear System By Gauss Method with full pivoting *
'* and correction process.                                          *
'* ---------------------------------------------------------------- *
'* SAMPLE RUN:                                                      *
'*  Solve Complex Linear System AX = B with:                        *
'*                                                                  *
'*      ( 47,-15) ( 62,5) (  0,-72) (61, 20)         ( 629,988)     *
'*  A = (  6, 14) (-17,3) (-102,91) ( 7,-12) and B = (-180,825)     *
'*      ( 13, 55) ( 32,8) (  41, 7) (25,  1)         ( 877,441)     *
'*      (111,25)  ( 40,0) ( 12,-82) (58,-30)         ( 734,-88)     *
'*                                                                  *
'* System solutions:                                                *
'*  ( -1.0000,  3.0000)                                             *
'*  (  2.0000, 10.0000)                                             *
'*  (  7.0000, -5.0000)                                             *
'*  ( 17.0000,  6.0000)                                             *
'*                                                                  *
'* ---------------------------------------------------------------- *
'* Reference:  "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 CSysLin

DEFINT I-N

OPTION BASE 1

DIM SHARED NMAX
  NMAX = 20

' Matc = Array(1..NMAX,1..NMAX) of Complex
' Vecc = Array(1..NMAX) of Complex
' Veci = Array(1..NMAX) of Integer

' I,it,M,N: Integer
' dta, eps: Real

' Common complex matrices and vectors
DIM SHARED A(NMAX, NMAX, 2)
DIM SHARED B(NMAX, 2), X(NMAX, 2)
DIM SHARED C(NMAX, NMAX, 2), W(NMAX, 2), Z(NMAX, 2)

  N = 4  'Size of linear system

  A(1, 1, 1) = 47!: A(1, 1, 2) = -15!
  A(1, 2, 1) = 62!: A(1, 2, 2) = 5!
  A(1, 3, 1) = 0!: A(1, 3, 2) = -72!
  A(1, 4, 1) = 61!: A(1, 4, 2) = 20!

  A(2, 1, 1) = 6!: A(2, 1, 2) = 14!
  A(2, 2, 1) = -17!: A(2, 2, 2) = 3!
  A(2, 3, 1) = -102!: A(2, 3, 2) = 91!
  A(2, 4, 1) = 7!: A(2, 4, 2) = -12!

  A(3, 1, 1) = 13!: A(3, 1, 2) = -55!
  A(3, 2, 1) = 32!: A(3, 2, 2) = 8!
  A(3, 3, 1) = 41!: A(3, 3, 2) = 7!
  A(3, 4, 1) = 25!: A(3, 4, 2) = 1!

  A(4, 1, 1) = 111!: A(4, 1, 2) = 25!
  A(4, 2, 1) = 40!: A(4, 2, 2) = 0!
  A(4, 3, 1) = 12!: A(4, 3, 2) = -82!
  A(4, 4, 1) = 58!: A(4, 4, 2) = -30!

  B(1, 1) = 629!: B(1, 2) = 988!
  B(2, 1) = -180!: B(2, 2) = 825!
  B(3, 1) = 877!: B(3, 2) = 441!
  B(4, 1) = 734!: B(4, 2) = -88!

  eps = 1E-10: dta = 1E-10: M = 1

  RSLCGTC eps, dta, M, N, it

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

  CLS
  PRINT
  PRINT " System solutions:"
  FOR J = 1 TO N
    PRINT USING F$; X(J, 1); X(J, 2)
  NEXT J
  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 other subroutines. 
  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

 '*************************************************************
 '* Solve a Complex Linear System AX = B By Gauss Method with *
 '* full pivoting and a correction process.                   *
 '* --------------------------------------------------------- *
 '* INPUTS:                                                   *
 '*         eps, dta : absolute and relative precisions       *
 '*         M : maximum number of iterations                  *
 '*         N : size of linear system                         *
 '*         A : complex matrix                                *
 '*         B : right-hand side (complex vector)              *
 '* OUTPUTS:                                                  *
 '*         it: flag, =-1 if no convergence, =0 if matrix A   *
 '*             is singular, =1 convergence ok.               *
 '*         X : contains system solution (X1,X2,...,Xn)       *
 '*                                                           *
 '*************************************************************
 'Note: A, B, X are shared with main program
  SUB RSLCGTC (eps, dta, M, N, it)
  DIM C0(2), C1(2), S(2), TMP(2), TMP1(2), Z0(2)  'Complex numbers
  DIM KP(NMAX), LP(NMAX)

    Z0(1) = 0!: Z0(2) = 0!
    
    TSCGT eps, N, it, KP(), LP()

    IF it = 1 THEN

      FOR I = 1 TO N
	W(I, 1) = B(I, 1): W(I, 2) = B(I, 2)
      NEXT I


      BSCGT N, KP(), LP()

      FOR I = 1 TO N
	X(I, 1) = Z(I, 1)
	X(I, 2) = Z(I, 2)
      NEXT I

      PRINT " Vector X after first BSCGT:"
      FOR I = 1 TO N
	PRINT " "; X(I, 1); " "; X(I, 2)
      NEXT I

      INPUT " ", R$

      it = -1: L = 1
      WHILE it = -1 AND L <= M
	phi1 = 0!
	FOR I = 1 TO N
	  TMP(1) = X(I, 1): TMP(2) = X(I, 2)
	  IF CABS(TMP()) > phi1 THEN phi1 = CABS(TMP())
	NEXT I
	FOR I = 1 TO N
	  S(1) = Z0(1): S(2) = Z0(2)
	  FOR J = 1 TO N
	    C0(1) = S(1): C0(2) = S(2)
	    'CMUL(A(I,J),X(J),C1)
	    TMP(1) = A(I, J, 1): TMP(2) = A(I, J, 2)
	    TMP1(1) = X(J, 1): TMP1(2) = X(J, 2)
	    CMUL TMP(), TMP1(), C1()
	    CADD C0(), C1(), S()
	  NEXT J
	  'CDIF(B(I),S,W(I))
	  W(I, 1) = B(I, 1) - S(1)
	  W(I, 2) = B(I, 2) - S(2)
	NEXT I

	BSCGT N, KP(), LP()

	FOR I = 1 TO N
	  C0(1) = X(I, 1): C0(2) = X(I, 2)
	  'CADD(C0,Z(I),X(I))
	  X(I, 1) = C0(1) + Z(I, 1)
	  X(I, 2) = C0(2) + Z(I, 2)
	NEXT I
	phi2 = 0!
	FOR I = 1 TO N
	  TMP(1) = Z(I, 1): TMP(2) = Z(I, 2)
	  IF CABS(TMP()) > phi2 THEN phi2 = CABS(TMP())
	NEXT I
	IF (phi2 / phi1) < dta THEN
	  it = 1
	ELSE
	  L = L + 1
	END IF
      WEND
    END IF
  END SUB 'SLCGTC

 '*****************************************************************
 '* 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.                    *
 '*****************************************************************
 'NOTE: C is shared with main program and other Subroutines.
  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