'*****************************************************
'*     Solving a complex linear system A Z = B,      *
'*     where A is complex matrix of size N*N,        *
'*     B is a complex matrix of size N*M,            *
'*     Z is a complex solution matrix of size N*M.   *
'* ------------------------------------------------- *
'* Ref.: "Mathematiques en Turbo-Pascal By M. Ducamp *
'* and A. Reverchon (vol 2), Eyrolles, Paris, 1988"  *
'* [BIBLI 05].                                       *
'* ------------------------------------------------- *
'* SAMPLE RUN:                                       *
'*                                                   *
'* Solve complex linear system:                      *
'*                                                   *
'*   z1 +  z2 + iz3 = 5                              *
'*   z1 + iz2 +  z3 = -2i                            *
'*  iz1 +  z2 +  z3 = 1                              *
'*                                                   *
'*  Here: N=3, M=1                                   *
'*                                                   *
'*  System solution:                                 *
'*                                                   *
'*  z1 =   1.500000  i  -0.500000                    *
'*  z2 =   1.000000  i   1.000000                    *
'*  z3 =  -0.500000  i  -2.500000                    *
'*                                                   *
'*  DET =   20                                       *
'*                                                   *
'*               Basic version by J-P Moreau, Paris. *
'*                        (www.jpmoreau.fr)          *
'*****************************************************
' Explanations:
' The original complex linear system to solve is:
'     ( a11 a12 ... a1N ) ( z1 )   ( b1 )
'     ( a21 a22 ... a2N ) ( z2 )   ( b2 )
'     ( ... ... ... ... ) (... ) = (... )
'     ( aN1 aN2 ... aNN ) ( zN )   ( bN )
'     a k,l = c k,l + i d k,l
'     z k = x k + i y k
'     b k = e k + i f k
'     where c, d, e, f, x, y are real numbers.
' The system is replaced by the following REAL system of size 2*N:
'     ( c11 -d11 c12 -d12 ... c1N -d1N ) ( x1 )   ( e1 )
'     ( d11  c11 d12  c12 ... d1N  c1N ) ( y1 )   ( f1 )
'     ( c21 -d21 c22 -d22 ... c2N -d2N ) ( x2 )   ( e2 )
'     ( d21  c21 d22  c22 ... d2N  c2N ) ( y2 ) = ( f2 )
'     ( ... .... ... .... ... ... .... ) ( ...)   ( ...)
'     ( cN1 -dN1 cN2 -dN2 ... cNN -dNN ) ( xN )   ( eN )
'     ( dN1  cN1 dN2  cN2 ... dNN  cNN ) ( yN )   ( fN )
' The real system is then solved by a classic Gauss-Jordan method.
'
'-----------------------------------------------------------------
defdbl a-h,o-z
defint i-n

NMAX=3 : MMAX=1

Dim  Xm(NMAX,NMAX), Ym(NMAX,NMAX)
Dim  Xv(NMAX), Yv(NMAX)

'used by subroutine 1000
Dim  A(2*NMAX,2*NMAX)
Dim  B(2*NMAX,MMAX)
Dim  PC(2*NMAX), PL(2*NMAX), CS(2*NMAX)


  li=NMAX: ico=MMAX
  'read Xm matrix                  real part of original A matrix
  Xm(1,1)= 1: Xm(1,2)=1: Xm(1,3)=0 '   1  1  0
  Xm(2,1)= 1: Xm(2,2)=0: Xm(2,3)=1 '   1  0  1
  Xm(3,1)= 0: Xm(3,2)=1: Xm(3,3)=1 '   0  1  1
  'read Ym matrix                  imaginary part of original A matrix
  Ym(1,1)= 0: Ym(1,2)=0: Ym(1,3)=1 '   0  0  1
  Ym(2,1)= 0: Ym(2,2)=1: Ym(2,3)=0 '   0  1  0
  Ym(3,1)= 1: Ym(3,2)=0: Ym(3,3)=0 '   1  0  0
  'read Xv vector                  real part of 1st column of original B matrix
  Xv(1)=5: Xv(2)=0: Xv(3)=1        '   5  0  1
  'read Yv vector             imaginary part of 1st column of original B matrix
  Yv(1)=0: Yv(2)=-2: Yv(3)=0       '   0 -2  0

  'The complex system AZ=B is replaced by real system MX=V of size 2*NMAX
  'create M matrix (name A for subroutine 1000)
  for i=1 to NMAX
    for j=1 to NMAX
      A(2*i-1,2*j-1)=Xm(i,j): A(2*i,2*j)=Xm(i,j)
      A(2*i,2*j-1)=Ym(i,j): A(2*i-1,2*j)=-Ym(i,j)
    next j
  next i
  'create V matrix (name B for subroutine 1000)
  for i=1 to NMAX
    B(2*i-1,1)=Xv(i): B(2*i,1)=Yv(i)
  next

  ' solve real linear system MX=V of zize 2*NMAX
  ' by Gauss-Jordan method - See program sysmat.bas
  gosub 1000  'call MATINV(li*2, co, M, V, det)

  ' write complex solutions stored in vector V (name B) in the form:
  '  x1,y1,x2,y2,...,xn,yn, with zk = xk + i yk...
  cls
  print
  print "  System solution:"
  print
  for i=1 to NMAX
    print using "  z# = ####.######  i  ####.######"; i; B(2*i-1,1); B(2*i,1)
  next
  print
  print "  DET = "; DET
  print
  print

END 'of main program


'*******************************************
'*  SOLVING A LINEAR MATRIX SYSTEM AX = B  *
'*  with Gauss-Jordan method using full    *
'*  pivoting at each step. During the pro- *
'* cess, original A and B matrices are     *
'* destroyed to spare storage location.    *
'* --------------------------------------- *
'* INPUTS:    A   MATRIX N*N               *
'*            B   MATRIX N*M               *
'* --------------------------------------- *
'* OUTPUS:    A   INVERSE OF A N*N         *
'*            DET  DETERMINANT OF A        *
'*            B   SOLUTION MATRIX N*M      *
'* --------------------------------------- *
'* NOTA - If M=0, inversion of A matrix    *
'*        only.                            *
'*******************************************
1000 'Labels:  100, 200, 300
  'Initializations
  N=2*li : M=ico
  DET = 1#
  FOR I = 1 TO N
    PC(I) = 0#
    PL(I) = 0#
    CS(I) = 0#
  NEXT I

  'Main loop
  FOR K = 1 TO N
    'Searching greatest pivot
    PV = A(K, K)
    IK = K: JK = K
    PAV = ABS(PV)
    FOR I = K TO N
      FOR J = K TO N
        temp = ABS(A(I, J))
        IF temp > PAV THEN
          PV = A(I, J)
          PAV = ABS(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:

    PC(K) = JK
    PL(K) = IK

    'DETERMINANT DET is actualised
    'If DET=0, ERROR MESSAGE and STOP
    'Machine dependant EPSMACH equals here 2E-16
    IF IK <> K THEN
      DET = -DET
    END IF
    IF JK <> K THEN
      DET = -DET
    END IF
    DET = DET * PV
    temp = ABS(DET)
    IF temp < EPSMACH THEN
      print #2,
      print #2, " The determinant equals ZERO !!!"
      END   'stop program
    END IF

    'POSITIONNING PIVOT IN K,K:}
    IF IK <> K THEN
      for I=1 to N
        'EXCHANGE LINES IK and K of matrix A:}
        TT = A(IK, I)
        A(IK, I) = A(K, I)
        A(K, I) = TT
      NEXT I
    END IF

    IF M <> 0 THEN
      FOR I = 1 TO M
        TT = B(IK, I)
        B(IK, I) = B(K, I)
        B(K, I) = TT
      NEXT I
    END IF

    'PIVOT is at correct line
    IF JK <> K THEN
      FOR I = 1 TO N
      'EXCHANGE COLUMNS JK and K of matrix A:
        TT = A(I, JK)
        A(I, JK) = A(I, K)
        A(I, K)=TT
      NEXT I
    END IF

    'The PIVOT is at correct column.
    'and is located in K,K.

    'Column K of matrix A is stored in CS vector
    'then column K is set to zero.
    for I=1 to N
      CS(I) = A(I, K)
      A(I, K) = 0#
    NEXT I

    CS(K) = 0#
    A(K, K) = 1#
    'Line K of matrix A is modified:
    temp = ABS(PV)
    IF temp < EPSMACH THEN
      print #2,
      print #2, " PIVOT TOO SMALL - STOP"
      END   'stop program
    END IF
    FOR I = 1 TO N
      A(K, I) = A(K, I) / PV
    NEXT I

    IF M <> 0 THEN
      FOR I = 1 TO M
        B(K, I) = B(K, I) / PV
      NEXT I
    END IF

    'Other lines of matrix A are modified:
    FOR J = 1 TO N
      IF J = K THEN
        GOTO 100
      END IF
      FOR I = 1 TO N
        'Line J of matrix A is modified:
        A(J, I) = A(J, I) - CS(J) * A(K, I)
      NEXT I
      IF M <> 0 THEN
        FOR I = 1 TO M
          B(J, I) = B(J, I) - CS(J) * B(K, I)
        NEXT I
      END IF
100 NEXT J
  'Line K is ready
  NEXT K 'end K loop

  'MATRIX A INVERSION IS DONE - REARRANGEMENT OF MATRIX A

  'EXCHANGE LINES
  FOR I = N TO 1 STEP -1
    IK = INT(PC(I))
    IF IK = I THEN
      GOTO 200
    END IF
    'EXCHANGE LINES I and PC(I) of matrix A:
    FOR J = 1 TO N
      TT = A(I, J)
      A(I, J) = A(IK, J)
      A(IK, J) = TT
    NEXT J
    IF M <> 0 THEN
      FOR J = 1 TO M
        TT = B(I, J)
        B(I, J) = B(IK, J)
        B(IK, J) = TT
      NEXT J
    END IF
    'NO MORE EXCHANGE is NECESSARY
    'Go to next line.
200 NEXT I 'End of loop i=N to 1 step -1

  'EXCHANGE COLUMNS
  FOR J = N TO 1 STEP -1
    JK = INT(PL(J))
    IF JK = J THEN
      GOTO 300
    END IF
    'EXCHANGE COLUMNS J ET PL(J) of matrix A:
    FOR I = 1 TO N
      TT = A(I, J)
      A(I, J) = A(I, JK)
      A(I, JK) = TT
    NEXT I
    'NO MORE EXCHANGE is NECESSARY
    'Go to next column.}
300 NEXT J
'REARRANGEMENT TERMINATED
RETURN  'End of subroutine 1000


' end of file zmatsys.bas