DECLARE SUB CMUL (Z1!(), Z2!(), Z!())
DECLARE SUB RSHCGT (eps!, N%, R0 AS INTEGER, M0%)
DECLARE FUNCTION CABS! (Z!())
DECLARE SUB CDIV (Z1!(), Z2!(), Z!())
'**************************************************************
'*   Solving a complex homogeneous linear system by Gauss     *
'*   Method with full pivoting.                               *
'* ---------------------------------------------------------- *
'* SAMPLE RUN:                                                *
'*                        | -1  i -i |                        *
'* Solve AX = 0 with  A = | -i -1  1 |                        *
'*                        |  i  1 -1 |                        *
'*                                                            *
'* Basic Family of system solutions:                          *
'*  Solution  1                                               *
'*   X 1  = ( 0 , 1 )                                         *
'*   X 2  = ( 1 , 0 )                                         *
'*   X 3  = ( 0 , 0 )                                         *
'*                                                            *
'*  Solution  2                                               *
'*   X 1  = ( 0 ,-1 )                                         *
'*   X 2  = ( 0 , 0 )                                         *
'*   X 3  = ( 1 , 0 )                                         *
'*                                                            *
'* ---------------------------------------------------------- *
'* Ref.: "Algèbre, Algorithmes et programmes en Pascal        *
'*        By Jean-Louis Jardrin, DUNOD Paris, 1988".          *
'*                                                            *
'*                 Quick Basic Release By J-P Moreau, Paris.  *
'*                             (www.jpmoreau.fr)              *
'*************************************************************}
'Program TRSHCGT

DEFINT I-N

OPTION BASE 1

DIM SHARED NMAX
  NMAX = 20

DIM R0 AS INTEGER

'Complex matriczs qhared with main subroutine RSHCGT
DIM SHARED A(NMAX, NMAX, 2), VX(NMAX, NMAX, 2)

  N = 3   'size of complex matrix A

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

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

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

  eps = 1E-10

  RSHCGT eps, N, R0, M0

  CLS
  PRINT
  IF R0 = N THEN
    PRINT " System solution:"
    FOR i = 1 TO N
      PRINT "  X"; i; " = ("; VX(i, 1, 1); ","; VX(i, 1, 2); ")"
    NEXT i
  ELSE
    PRINT " Basic Family of system solutions:"
    FOR l = 1 TO M0
      PRINT " Solution "; l
      FOR i = 1 TO N
	PRINT "  X"; i; " = ("; VX(i, l, 1); ","; VX(i, l, 2); ")"
      NEXT i
    NEXT l
  END IF
  PRINT

END 'of main program

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

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

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

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

  SUB CSwap (Z1(), Z2())
  DIM T(2)
  T(1) = Z2(1): T(2) = Z2(2)
  Z2(1) = Z1(1): Z2(2) = Z1(2)
  Z1(1) = T(1): Z1(2) = T(2)
  END SUB

 '**********************************************************
 '* Procedure RSHCGT solves the complex homogeneous linear *
 '* system AX = 0 by Gauss method with full pivoting.      *
 '* ------------------------------------------------------ *
 '* INPUTS:                                                *
 '*         eps:  required precision                       *
 '*          N :  size of complex matrix A (NxN)           *
 '*          A :  Complex Matrix shared with main program  *
 '* OUTPUTS:                                               *
 '*          R0:  rank of A in output                      *
 '*          M0:  dimension of system solutions space      *
 '*               (if R0 <> N)                             *
 '*          VX:  contains in output the M0 solution       *
 '*               vectors (stored in columns 1..M0) when   *
 '*               R0<>N, or the unique solution vector     *
 '*               (stored in column 1) when R0=N.          *
''*               Shared with mainprogram.                 *
 '**********************************************************
  SUB RSHCGT (eps, N, R0 AS INTEGER, M0)
  DIM C0(2), C1(2), P0(2), Q0(2), S(2), TMP(2), TMP1(2)'Complex numbers
  DIM Z0(2), Z1(2)
  DIM I0(NMAX)

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

    FOR K = 1 TO N
      I0(K) = K
    NEXT K
    R0 = 0: K = 1
    WHILE R0 = 0 AND K <= N
      P0(1) = A(K, K, 1): P0(2) = A(K, K, 2)
      l0 = K: k0 = K
      FOR i = K TO N
	FOR j = K TO N
	  TMP(1) = A(i, j, 1): TMP(2) = A(i, j, 2)
	  IF CABS(TMP()) > CABS(P0()) THEN
	    P0(1) = A(i, j, 1): P0(2) = A(i, j, 2)
	    l0 = i: k0 = j
	  END IF
	NEXT j
      NEXT i
      IF CABS(P0()) < eps THEN
	R0 = K - 1
	M0 = N - R0
      ELSE
	IF K = N THEN
	  R0 = N: M0 = 0
	ELSE
	  IF l0 <> K THEN
	    FOR j = K TO N
	      'CSwap(A(k,j),A(l0,j))
	      TMP(1) = A(l0, j, 1): TMP(2) = A(l0, j, 2)
	      A(l0, j, 1) = A(K, j, 1): A(l0, j, 2) = A(K, j, 2)
	      A(K, j, 1) = TMP(1): A(K, j, 2) = TMP(2)
	    NEXT j
	  END IF
	  IF k0 <> K THEN
	    FOR i = 1 TO N
	      'CSwap(A(i,k),A(i,k0))
	      TMP(1) = A(i, k0, 1): TMP(2) = A(i, k0, 2)
	      A(i, k0, 1) = A(i, K, 1): A(i, k0, 2) = A(i, K, 2)
	      A(i, K, 1) = TMP(1): A(i, K, 2) = TMP(2)
	    NEXT i
	  END IF
	  N0 = I0(K): I0(K) = I0(k0): I0(k0) = N0
	END IF
      END IF

      FOR i = K + 1 TO N
	'Q0 = A(i,k)/P0
	TMP(1) = A(i, K, 1): TMP(2) = A(i, K, 2)
	CDIV TMP(), P0(), Q0()
	FOR j = K + 1 TO N
	  C0(1) = A(i, j, 1): C0(2) = A(i, j, 2)
	  TMP1(1) = A(K, j, 1): TMP1(2) = A(K, j, 2)
	  'C1 = Q0 * A(k,j)
	  CMUL Q0(), TMP1(), C1()
	  A(i, j, 1) = C0(1) - C1(1)
	  A(i, j, 2) = C0(2) - C1(2)
	NEXT j
      NEXT i
      K = K + 1
    WEND

    IF R0 = N THEN
      FOR i = 1 TO N
	VX(i, 1, 1) = Z0(1): VX(i, 1, 2) = Z0(2)
      NEXT i
    ELSE
      FOR l = 1 TO M0
	FOR i = N TO R0 + 1 STEP -1
	  IF i = R0 + l THEN
	    VX(i, l, 1) = Z1(1): VX(i, l, 2) = Z1(2)
	  ELSE
	    VX(i, l, 1) = Z0(1): VX(i, l, 2) = Z0(2)
	  END IF
	NEXT i
	FOR i = R0 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)
	    'C1 = A(i,j) * VX(j,l)
	    TMP(1) = A(i, j, 1): TMP(2) = A(i, j, 2)
	    TMP1(1) = VX(j, l, 1): TMP1(2) = VX(j, l, 2)
	    CMUL TMP(), TMP1(), C1()
	    S(1) = C0(1) + C1(1)
	    S(2) = C0(2) + C1(2)
	  NEXT j
	  'C0 = S/A(i,i)
	  TMP(1) = A(i, i, 1): TMP(2) = A(i, i, 2)
	  CDIV S(), TMP(), C0()
	  VX(i, l, 1) = -C0(1)
	  VX(i, l, 2) = -C0(2)
	NEXT i
      NEXT l
      FOR i = 1 TO N - 1
	j = i + 1
	WHILE j <= N
	  IF I0(j) = i THEN
	    FOR l = 1 TO M0
	      'CSwap(VX(j,l),VX(i,l))
	      TMP(1) = VX(i, l, 1): TMP(2) = VX(i, l, 2)
	      VX(i, l, 1) = VX(j, l, 1): VX(i, l, 2) = VX(j, l, 2)
	      VX(j, l, 1) = TMP(1): VX(j, l, 2) = TMP(2)
	    NEXT l
	    I0(j) = I0(i): I0(i) = i
	    j = N + 1
	  ELSE
	    j = j + 1
	  END IF
	WEND
      NEXT i
    END IF
  END SUB  'RSHCGT