DECLARE FUNCTION CABS! (Z!())
DECLARE SUB CDIF (Z1!(), Z2!(), Z!())
DECLARE SUB MATPRINT2 (titre$, N%, M%)
DECLARE SUB MATREAD (N%)
DECLARE SUB MATPRINT (titre$, N%)
DECLARE SUB MATREAD1 (N%, M%)
DECLARE SUB MATPRINT1 (titre$, N%, M%)
DECLARE SUB CMATINV (N%, M%, DET!())
DECLARE SUB MATMUL1 (N%, M%)
DECLARE SUB MATMUL (N%)
DECLARE SUB MATPRINT3 (titre$, N%)
DECLARE SUB CADD (Z1!(), Z2!(), Z!())
DECLARE SUB CMUL (Z1!(), Z2!(), Z!())
DECLARE SUB CDIV (Z1!(), Z2!(), Z!())
'***********************************************************************************************
'*            SOLVING A COMPLEX LINEAR MATRIX SYSTEM AX=B BY GAUSS-JORDAN METHOD               *
'* ------------------------------------------------------------------------------------------- *
'*                                                  Quick Basic Version By J-P Moreau, Paris   *
'*                                                              (www.jpmoreau.fr)              *
'* SAMPLE RUN:                                                                                 *
'* INPUTS  from example file csysmat.dat,                                                      *
'* OUTPUTS to file csysmat.lst.                                                                *
'* ------------------------------------------------------------------------------------------- *
'* SAMPLE RUN:                                                                                 *
'* Input file                                                                                  *
'* 4                                                                                           *
'* 47.0 -15.0  62.0 5.0    0.0 -72.0 61.0  20.0                                                *
'*  6.0  14.0 -17.0 3.0 -102.0  91.0  7.0 -12.0                                                *
'* 13.0 -55.0  32.0 8.0   41.0   7.0 25.0   1.0                                                *
'*111.0  25.0  40.0 0.0   12.0 -82.0 58.0 -30.0                                                *
'* 1                                                                                           *
'* 629.0 988.0                                                                                 *
'*-180.0 825.0                                                                                 *
'* 877.0 441.0                                                                                 *
'* 734.0 -88.0                                                                                 *
'*                                                                                             *
'* Output file                                                                                 *
'*-------------------------------------------------------                                      *
'*  Solving a complex linear matrix system AX = B                                              *
'*      (by GAUSS-JORDAN method)\n\n                                                           *
'*-------------------------------------------------------                                      *
'*  N =  4                                                                                     *
'*                                                                                             *
'*  MATRIX A:                                                                                  *
'*                                                                                             *
'* ( 47 ,-15 ) ( 62 , 5 ) ( 0 ,-72 ) ( 61 , 20 )                                               *
'* ( 6 , 14 ) (-17 , 3 ) (-102 , 91 ) ( 7 ,-12 )                                               *
'* ( 13 ,-55 ) ( 32 , 8 ) ( 41 , 7 ) ( 25 , 1 )                                                *
'* ( 111 , 25 ) ( 40 , 0 ) ( 12 ,-82 ) ( 58 ,-30 )                                             *
'*                                                                                             *
'*  M =  1                                                                                     *
'*                                                                                             *
'*  MATRIX B :                                                                                 *
'*                                                                                             *
'* ( 629 , 988 )                                                                               *
'* (-180 , 825 )                                                                               *
'* ( 877 , 441 )                                                                               *
'* ( 734 ,-88 )                                                                                *
'*                                                                                             *
'*  INVERSE MATRIX:                                                                            ****
'*                                                                                                *
'* ( 3.940673E-03 ,-7.98777E-03 ) (-8.374093E-03 ,-4.057823E-03 ) (-9.490407E-03 , 1.793796E-02 ) *
'* ( 2.100789E-04 ,-1.660345E-03 )                                                                *
'* ( 3.678096E-02 , 2.182006E-02 ) ( 1.479263E-02 ,-2.707624E-02 ) (-3.379653E-02 ,-3.512894E-02 )*
'* (-5.015412E-03 ,-1.617311E-02 )                                                                *
'* (-4.225129E-03 ,-5.070878E-03 ) (-7.673959E-03 ,-1.820856E-03 ) ( 3.482056E-03 , 4.805251E-03 )*
'* ( 6.924988E-04 , 3.649074E-03 )                                                                *
'* (-1.972743E-02 ,-.0165881 ) (-1.48955E-03 , 1.880551E-02 ) ( 3.373942E-02 , 1.536814E-02 )     *
'* ( 5.363459E-03 , 1.723915E-02 )                                                             ****
'*                                                                                             *
'*  SOLUTION MATRIX X:                                                                         *
'*                                                                                             *
'* (-.9999988 , 3 )                                                                            *
'* ( 1.999996 , 10 )                                                                           *
'* ( 7 ,-5 )                                                                                   *
'* ( 17 , 6 )                                                                                  *
'*                                                                                             *
'*  DETERMINANT= ( 1.741656E+07 ,-1.059832E+07 )                                               *
'*                                                                                             *
'*  VERIFICATION AX = B:                                                                       *
'*                                                                                             *
'* ( 629 , 988.0001 )                                                                          *
'* (-180 , 825 )                                                                               *
'* ( 877 , 441 )                                                                               *
'* ( 734.0001 ,-88.00003 )                                                                     *
'*                                                                                             *
'* End of results file.                                                                        *
'*-------------------------------------------------------                                      *
'*                                                                               Feb. 2010     *
'***********************************************************************************************
DEFINT I-N       'Integers begin with I..N

OPTION BASE 0    'Table index start from zero

'variables shared with subroutines 

dim shared NMAX
      NMAX = 20
DIM SHARED MMAX
      MMAX = 4

DIM SHARED A(NMAX, NMAX, 2), A1(NMAX, NMAX, 2), D(NMAX, NMAX, 2)
DIM SHARED B(NMAX, MMAX, 2), C(NMAX, MMAX, 2)


DIM DETER(2)   'complex determinant of A
					    
      CLS
      PRINT
      INPUT " Input data file name (without .dat): ", nom$
									       
' READ MATRIX A(N,N) TO BE INVERTED                                            
' and COPY to OUTPUT FILE                                                          
      
      s$ = nom$ + ".dat"
      OPEN s$ FOR INPUT AS #1

      s1$ = nom$ + ".lst"
      OPEN s1$ FOR OUTPUT AS #2

      PRINT #2, "-------------------------------------------------------"
      PRINT #2, "  Solving a complex linear matrix system AX = B"
      PRINT #2, "      (By GAUSS-JORDAN method)\n\n"
      PRINT #2, "-------------------------------------------------------"
	    
      INPUT #1, N
      PRINT #2, "  N = "; N
      
      MATREAD N
      MATPRINT "  MATRIX A:", N
									 
									       
' READ RIGHT SIDE MATRIX B(N,M), IF M<>0                                            

      INPUT #1, M
      print #2,
      PRINT #2, "  M = "; M
      
      IF M <> 0 THEN
	MATREAD1 N, M
	MATPRINT1 "  MATRIX B :", N, M
      END IF
' END DATA SECTION    
      CLOSE #1
									       
' STORING A IN A1 (Original A is lost when calling CMATINV)                                                   

      FOR I = 0 TO N - 1
	FOR J = 0 TO N - 1
	  A1(I, J, 0) = A(I, J, 0)
	  A1(I, J, 1) = A(I, J, 1)
	NEXT J
      NEXT I

' CALL INVERSION ROUTINE                                             
									       
      CMATINV N, M, DETER()
									       
' PRINT RESULTS (INVERSE MATRIX ONLY IF M=0)                                                        
									       
      MATPRINT "  INVERSE MATRIX:", N
										
      IF M <> 0 THEN
	MATPRINT1 "  SOLUTION MATRIX X:", N, M
      END IF

      PRINT #2,
      PRINT #2, "  DETERMINANT= ("; DETER(0); ","; DETER(1); ")"

' VERIFICATION THAT NEW PRODUCT A*B EQUALS ORIGINAL B MATRIX             
' (OPTIONAL)                                                              
									       
      IF M <> 0 THEN
' CALL MULTIPLICATION ROUTINE                                       
	MATMUL1 N, M
	 
	MATPRINT2 "  VERIFICATION AX = B:", N, M
      ELSE
' CALL SQUARE MULTIPLICATION ROUTINE
	MATMUL N
	  
	MATPRINT3 "  VERIFICATION A1 * A = I :", N
      END IF

      PRINT #2,
      PRINT #2, " End of results file."
      PRINT #2, "-------------------------------------------------------"
      CLOSE #2

      PRINT
      PRINT " Results in file "; s1$

END 'of main program


  FUNCTION CABS (Z())
    CABS = SQR(Z(0) ^ 2 + Z(1) ^ 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(0) ^ 2 + Z2(1) ^ 2
    IF D < 2E-12 THEN
      PRINT #2, " Complex Divide by zero!"
    ELSE
      C1(0) = Z2(0): C1(1) = -Z2(1)
      CMUL Z1(), C1(), Z()
      Z(0) = Z(0) / D: Z(1) = Z(1) / D
    END IF
  END SUB

'************************************************                                   
'* SOLVING A COMPLEX MATRIX LINEAR SYSTEM AX=B  *
'* with Gauss-Jordan method using full pivoting *
'* at each step. During the process, original   *
'* AA and BB matrices are destroyed to spare    *
'* storage location.                            *
'* -------------------------------------------- *
'* INPUTS:    A    COMPLEX MATRIX N*N           *
'*            B    COMPLEX MATRIX N*M           *                                   
'* -------------------------------------------- *                                   
'* OUTPUTS:   A    INVERSE OF AA N*N            *                                   
'*            DET  COMPLEX DETERMINANT OF A     *                                   
'*            B    SOLUTION COMPLEX MATRIX N*M  *
'* -------------------------------------------- *                                   
'* NOTA - If M=0 inversion of A matrix only.    *                                   
'************************************************
'Note: A and B are shared with main program
SUB CMATINV (N, M, DET())
DIM IPC(NMAX), IPL(NMAX)
DIM CS(NMAX, 2)
DIM PV(2), temp(2), temp1(2), temp2(2)

      EPSMACH = 2E-12

' Initializations:                                                           
      DET(0) = 1!: DET(1) = 0!
      FOR I = 0 TO N - 1
	IPC(I) = 0
	IPL(I) = 0
	CS(I, 0) = 0!
	CS(I, 1) = 0!
      NEXT I
' Main loop                                                                       
      FOR K = 0 TO N - 1
' Searching greatest pivot:                                             
	PV(0) = A(K, K, 0)
	PV(1) = A(K, K, 1)
	IK = K
	JK = K
	PAV = CABS(PV())
	FOR I = K TO N - 1
	  FOR J = K TO N - 1
	    temp(0) = A(I, J, 0)
	    temp(1) = A(I, J, 1)
	    TMP = CABS(temp())
	    IF TMP > PAV THEN
	      PV(0) = A(I, J, 0)
	      PV(1) = A(I, J, 1)
	      PAV = CABS(PV())
	      IK = I
	      JK = J
	    END IF
	  NEXT J
	NEXT I
									       
' Search terminated, the pivot is in location I=IK, J=JK.                   
' Memorizing pivot location:

	IPC(K) = JK
	IPL(K) = IK
	
' DETERMINANT DET is actualised                                            
' If DET=0, ERROR MESSAGE and STOP                                         
' Machine dependant EPSMACH equals here 2e-12
									       
	IF IK <> K THEN
	  DET(0) = -DET(0): DET(1) = -DET(1)
	END IF

	IF JK <> K THEN
	  DET(0) = -DET(0): DET(1) = -DET(1)
	END IF

	CMUL DET(), PV(), temp()
	DET(0) = temp(0): DET(1) = temp(1)
	TMP = CABS(DET())
	IF TMP < EPSMACH THEN
	  PRINT #2,
	  PRINT #2, "  The complex determinant is too small!"
	  PRINT #2,
	  EXIT SUB
	END IF

' POSITIONNING PIVOT IN K,K:          
	IF IK <> K THEN
' EXCHANGE LINES IK and K of matrices AA and BB               
	  FOR I = 0 TO N - 1
	    'CSwap(A(IK,I),A(K,I))
	    temp(0) = A(K, I, 0): temp(1) = A(K, I, 1):
	    A(K, I, 0) = A(IK, I, 0): A(K, I, 1) = A(IK, I, 1)
	    A(IK, I, 0) = temp(0): A(IK, I, 1) = temp(1)
	  NEXT I
	  IF M <> 0 THEN
	    FOR I = 0 TO M - 1
	      'CSwap(B(IK,I),B(K,I))
	      temp(0) = B(K, I, 0): temp(1) = B(K, I, 1):
	      B(K, I, 0) = B(IK, I, 0): B(K, I, 1) = B(IK, I, 1)
	      B(IK, I, 0) = temp(0): B(IK, I, 1) = temp(1)
	    NEXT I
	  END IF
	END IF
									   
' PIVOT is at correct line              
	IF JK <> K THEN
' EXCHANGE COLUMNS JK and K of matrix AA:                     
	  FOR I = 0 TO N - 1
	    'CSwap(AA[I][JK],AA[I][K])
	    temp(0) = A(I, K, 0): temp(1) = A(I, K, 1):
	    A(I, K, 0) = A(I, JK, 0): A(I, K, 1) = A(I, JK, 1)
	    A(I, JK, 0) = temp(0): A(I, JK, 1) = temp(1)
	  NEXT I
	END IF
' The PIVOT is at correct column and is located in K,K.                                                 
									       
' Column K of matrix AA is stored in CS complex vector                           
' then column K is set to zero.                 
	FOR I = 0 TO N - 1
	  CS(I, 0) = A(I, K, 0)
	  CS(I, 1) = A(I, K, 1)
	  A(I, K, 0) = 0!
	  A(I, K, 1) = 0!
	NEXT I
									       
	CS(K, 0) = 0!: CS(K, 1) = 0!
	A(K, K, 0) = 1!: A(K, K, 1) = 0!

' Line K of matrix AA is modified:            
	TMP = CABS(PV())
	IF TMP < EPSMACH THEN
	  PRINT #2,
	  PRINT #2, " COMPLEX PIVOT TOO SMALL - STOP"
	  PRINT #2,
	  EXIT SUB
	END IF
	FOR I = 0 TO N - 1
'         A(K,I)=A(K,I)/PV
	  temp1(0) = A(K, I, 0): temp1(1) = A(K, I, 1)
	  CDIV temp1(), PV(), temp()
	  A(K, I, 0) = temp(0)
	  A(K, I, 1) = temp(1)
	NEXT I
									   
	IF M <> 0 THEN
	  FOR I = 0 TO M - 1
'           B(K,I)=B(K,I)/PV
	    temp1(0) = B(K, I, 0): temp1(1) = B(K, I, 1)
	    CDIV temp1(), PV(), temp()
	    B(K, I, 0) = temp(0)
	    B(K, I, 1) = temp(1)
	  NEXT I
	END IF

' Other lines of matrix AA are modified:                
	FOR J = 0 TO N - 1
	  IF J = K THEN J = J + 1
	  FOR I = 0 TO N - 1
' Line J of matrix AA is modified:
'           A(J,I)=A(J,I)-CS(J)*A(K,I)
	    temp1(0) = CS(J, 0): temp1(1) = CS(J, 1)
	    temp2(0) = A(K, I, 0): temp2(1) = A(K, I, 1)
	    CMUL temp1(), temp2(), temp()
	    temp2(0) = A(J, I, 0): temp2(1) = A(J, I, 1)
	    CDIF temp2(), temp(), temp1()
	    A(J, I, 0) = temp1(0)
	    A(J, I, 1) = temp1(1)
	  NEXT I
	  IF M <> 0 THEN
	    FOR I = 0 TO M - 1
'             B(J,I)=B(J,I)-CS(J)*B(K,I)
	      temp1(0) = CS(J, 0): temp1(1) = CS(J, 1)
	      temp2(0) = B(K, I, 0): temp2(1) = B(K, I, 1)
	      CMUL temp1(), temp2(), temp()
	      temp2(0) = B(J, I, 0): temp2(1) = B(J, I, 1)
	      CDIF temp2(), temp(), temp1()
	      B(J, I, 0) = temp1(0)
	      B(J, I, 1) = temp1(1)
	    NEXT I
	  END IF
	NEXT J
' Line K is ready
      NEXT K       'end of main loop

' COMPLEX MATRIX AA INVERSION IS DONE - REARRANGEMENT OF MATRIX A
									       
' EXCHANGE LINES                                                                      
      FOR I = N - 1 TO 0 STEP -1
	IK = IPC(I)
	IF IK = I THEN GOTO 10
' EXCHANGE LINES I and PC(I) of matrix AA:
	FOR J = 0 TO N - 1
	  'CSwap(A(I,J),A(IK,J))
	  temp(0) = A(IK, J, 0): temp(1) = A(IK, J, 1):
	  A(IK, J, 0) = A(I, J, 0): A(IK, J, 1) = A(I, J, 1)
	  A(I, J, 0) = temp(0): A(I, J, 1) = temp(1)
	NEXT J
	IF M <> 0 THEN
	  FOR J = 0 TO M - 1
	    'CSwap(B(I,J),B(IK,J))
	    temp(0) = B(IK, J, 0): temp(1) = B(IK, J, 1):
	    B(IK, J, 0) = B(I, J, 0): B(IK, J, 1) = B(I, J, 1)
	    B(I, J, 0) = temp(0): B(I, J, 1) = temp(1)

	  NEXT J
	END IF
' NO MORE EXCHANGE is NECESSARY                                                    
' Go to next line                                               
10    NEXT I
									       
' EXCHANGE COLUMNS:     
      FOR J = N - 1 TO 0 STEP -1
	JK = IPL(J)
	IF JK = J THEN GOTO 20
' EXCHANGE COLUMNS J ET PL(J) of matrix AA:
	FOR I = 0 TO N - 1
	  'CSwap(A(I,J),A(I,JK))
	  temp(0) = A(I, JK, 0): temp(1) = A(I, JK, 1):
	  A(I, JK, 0) = A(I, J, 0): A(I, JK, 1) = A(I, J, 1)
	  A(I, J, 0) = temp(0): A(I, J, 1) = temp(1)
	NEXT I
' NO MORE EXCHANGE is NECESSARY                                                    
' Go to next column.
20    NEXT J
' REARRANGEMENT TERMINATED.                                                      
  
END SUB 'CMATINV                                                                     

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

'*******************************************
'*   MULTIPLICATION OF TWO SQUARE COMPLEX  *
'*   MATRICES                              *
'* --------------------------------------- *
'* INPUTS:    A  MATRIX N*N                *
'*            A1 MATRIX N*N                *
'*            N  INTEGER                   *
'* --------------------------------------- *
'* OUTPUTS:   D  MATRIX N*N PRODUCT A*A1   *
'*                                         *
'*******************************************
SUB MATMUL (N)
DIM SUM(2), TMP(2), TMP1(2), TMP2(2)
  FOR I = 0 TO N - 1
    FOR J = 0 TO N - 1
      SUM(0) = 0!: SUM(1) = 0!
      FOR K = 0 TO N - 1
'       SUM=SUM+A(I,K)*A1(K,J)
	TMP1(0) = A(I, K, 0): TMP1(1) = A(I, K, 1)
	TMP2(0) = A1(K, J, 0): TMP2(1) = A1(K, J, 1)
	CMUL TMP1(), TMP2(), TMP()
	CADD SUM(), TMP(), TMP1()
	SUM(0) = TMP1(0): SUM(1) = TMP1(1)
      NEXT K
      D(I, J, 0) = SUM(0)
      D(I, J, 1) = SUM(1)
    NEXT J
  NEXT I
END SUB

'*******************************************
'*   MULTIPLICATION OF TWO REAL MATRICES   *
'* --------------------------------------- *
'* INPUTS:    A1 MATRIX N*N                *
'*            B  MATRIX N*M                *
'*            N  INTEGER                   *
'*            M  INTEGER                   *
'* --------------------------------------- *
'* OUTPUTS:   C  MATRIX N*M PRODUCT A*B    *
'*                                         *
'*******************************************
SUB MATMUL1 (N, M)
DIM SUM(2), TMP(2), TMP1(2), TMP2(2)
  FOR I = 0 TO N - 1
    FOR J = 0 TO M - 1
      SUM(0) = 0!: SUM(1) = 0!
      FOR K = 0 TO N - 1
'       SUM=SUM+A(I,K)*B(K,J)
	TMP1(0) = A1(I, K, 0): TMP1(1) = A1(I, K, 1)
	TMP2(0) = B(K, J, 0): TMP2(1) = B(K, J, 1)
	CMUL TMP1(), TMP2(), TMP()
	CADD SUM(), TMP(), TMP1()
	SUM(0) = TMP1(0): SUM(1) = TMP1(1)
      NEXT K
      C(I, J, 0) = SUM(0)
      C(I, J, 1) = SUM(1)
    NEXT J
  NEXT I
END SUB

' Print matrix A
  SUB MATPRINT (titre$, N)
     PRINT #2,
     PRINT #2, titre$
     PRINT #2,
     FOR I = 0 TO N - 1
	FOR J = 0 TO N - 1
	  PRINT #2, " ("; A(I, J, 0); ","; A(I, J, 1); ")";
	  IF J = N - 1 THEN PRINT #2,
	NEXT J
     NEXT I
  END SUB

' Print matrix B
  SUB MATPRINT1 (titre$, N, M)
     PRINT #2,
     PRINT #2, titre$
     PRINT #2,
     FOR I = 0 TO N - 1
	FOR J = 0 TO M - 1
	  PRINT #2, " ("; B(I, J, 0); ","; B(I, J, 1); ")";
	  IF J = M - 1 THEN PRINT #2,
	NEXT J
     NEXT I
  END SUB

' Print matrix C
  SUB MATPRINT2 (titre$, N, M)
     PRINT #2,
     PRINT #2, titre$
     PRINT #2,
     FOR I = 0 TO N - 1
	FOR J = 0 TO M - 1
	  PRINT #2, " ("; C(I, J, 0); ","; C(I, J, 1); ")";
	  IF J = M - 1 THEN PRINT #2,
	NEXT J
     NEXT I
  END SUB

' Print matrix D
  SUB MATPRINT3 (titre$, N)
     PRINT #2,
     PRINT #2, titre$
     PRINT #2,
     FOR I = 0 TO N - 1
	FOR J = 0 TO N - 1
	  PRINT #2, " ("; D(I, J, 0); ","; D(I, J, 1); ")";
	  IF J = N - 1 THEN PRINT #2,
	NEXT J
     NEXT I
  END SUB

'Read matrix A from input file
SUB MATREAD (N)
  FOR I = 0 TO N - 1
    FOR J = 0 TO N - 1
      INPUT #1, x, y
      A(I, J, 0) = x
      A(I, J, 1) = y
    NEXT J
  NEXT I
END SUB

'Read matrix B from input file
SUB MATREAD1 (N, M)
  FOR I = 0 TO N - 1
    FOR J = 0 TO M - 1
      INPUT #1, x, y
      B(I, J, 0) = x
      B(I, J, 1) = y
    NEXT J
  NEXT I
END SUB