C	WESTERN MICHIGAN UNIVERSITY
C	CMINV.F4 (FILE NAME ON LIBRARY DECTAPE)
C	CMINV, 2.3.1 (CALLING NAME, SUBLST. #)
C	SOLUTION OF LINEAR EQS., COMPLEX COEFFS. ALLOWED
C	THIS PROGRAM IS A COMBINATION OF ONE GIVEN BY WAYNE STATE
C	 UNIVERSITY WITH REVISIONAL AND ADDITIONAL PROGRAMMING
C	 BY MRS. EVA GAINES AND B. GRANET.
C	FORWMU PROGRAMS USED:  TTYPTY, DEVCHG, EXISTS, PRINTS
C	LIBRARY DECTAPE PROGRAMS USED:  USAGE.MAC
C	APLIB PROGRAMS USED:  IOB, GETFOR
C	INTERNAL SUBR. USED:  MULT, DIV, COMPNV
C	ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
C---------------NO VARIABLES ARE DECLARED AS COMPLEX IN THIS PROG.
C---------------THIS PROG. PERFORMS COMPLEX ARITHMETIC BY ASSIGNING
C--------------- A TO REAL PART AND B TO IMAG. PART OF COMPLEX NUMBERS.
C--------------- MULT. AND DIV. OF COMPLEX NOS. ARE DONE BY SUBR. MULT.
C--------------- DIV. RESP. USING REAL ARITHMETIC.  COMPLEX ADDITION
C---------------IS DONE IN SUBR. COMPNV IN ST. #1013+5 AND +6 AND
C--------------- AND ST. #33 AND ST. #33-1.
      DIMENSION A(30,40),B(30,40),AID(30,30),BID(30,30),ICOL(30),
     1IFMT(96)
      DIMENSION ID(16)
	COMMON/IOBLK/IDLG,INT,INP,IOUT,IDEV,IDEVA,IB,ICODE,NAMI(2)
      IOUT=2
      IDLG=-1
      INT=5
C--------------- TTYPTY RETURNS ZERO - TTY JOB, MINUS ONE - BATCH JOB
      CALL TTYPTY(ICODE)
      INP=1
	WRITE(IDLG,1000)
1000	FORMAT(1X,'WMU COMPLEX MATRIX INVERSION',/)
C	CALL USAGE('CMINV')
      CALL IOB(1)
8	CALL IOB(0)
      JX=30
3     CONTINUE
      WRITE(IDLG,14041)
14041 FORMAT(' IF YOU DESIRE OUTPUT IDENTIFICATION,ENTER UP TO 80'/
     1' CHARACTERS. OTHERWISE, JUST RETURN CARRIAGE'/)
      READ(INT,99)(ID(I),I=1,16)
99    FORMAT(16A5)
      WRITE(IDLG,1001)
1001  FORMAT(' ENTER THE NUMBER OF UNKNOWNS'/)
      READ(INT,5)N
      WRITE(IDLG,11101)
11101 FORMAT(' ENTER THE NUMBER OF SYSTEMS'/)
      READ(INT,5)M
      WRITE(IDLG,20120)
20120 FORMAT(' ENTER TOLERANCE'/)
      READ(INT,10)E
5     FORMAT(2I)
      CALL GETFOR(IDLG,INT,IFMT,ISTD,96,4)
44401 NP1=N+1
      LIM=M+N
      K=0
      DETR=0.
      DETI=0.
      IF(IDEV.NE.'TTY')GO TO 1
10    FORMAT(8G)
      WRITE(IDLG,2001)
2001  FORMAT(' ENTER DATA'/,
     1	' NOTE: USE COMMA(,) INSTEAD OF SEMICOLON(;) TO SEPARATE PAIRS'/
     1	' OF COMPLEX NUMBERS',/)
4     IF(ISTD-1)86,6,86
C---------------A(I,J) CONTAINS REAL PART OF ALL THE COEFFS. OF LEFT
C--------------- AND RT. HAND SIDES.  B(I,J) CONTAINS COMPLEX PART OF
C--------------- ALL COEFFS. OF LEFT AND RT. HAND SIDES.  BOTH ARE
C--------------- IN ROW ORDER I.E. FIRST ROW FIRST ETC.
6     READ(INP,10)((A(I,J),B(I,J),J=1,LIM),I=1,N)
      GO TO 144
86    READ(INP,IFMT)((A(I,J),B(I,J),J=1,LIM),I=1,N)
7001  FORMAT(1H0)
144   WRITE(IOUT,25052)(ID(I),I=1,16)
25052 FORMAT(1H-,16A5/)
      WRITE(IOUT,11)
11    FORMAT (1H1,23X,13H INPUT MATRIX /)
      DO 100 I=1,N
100   WRITE(IOUT,7)(A(I,J),B(I,J),J=1,LIM)
7     FORMAT(1H ,2E12.4,2X,2E12.4,2X,2E12.4)
      CALL COMPNV(A,B,AID,BID,M,N,E,DETR,DETI,K,30,ICOL)
      IF(DETR) 12,50,12
50    IF(DETI) 12,40,12
12    WRITE(IOUT,15)DETR,DETI
15    FORMAT (//,' REAL PART OF DETERMINANT=',E15.8,20X/ ,' IMAGINARY',
     1 ' PART OF DETERMINANT=',E15.8//)
      WRITE(IOUT,16)
16    FORMAT(23X,23HINVERSE OF INPUT MATRIX/)
      DO 20 I=1,N
C---------------THE FIRST N ROWS AND N COLMNS OF IAD AND BID
C--------------- CONSTITUTES THE INVERSE OF THE FIRST N ROWS AND N COLS.
C--------------- OF A AND B.   AID AND BID HAVE REAL AND IMAG. PARTS
C--------------- OF THE ELEMENTS OF THE INVERSE.
20    WRITE(IOUT,7)(AID(I,J),BID(I,J),J=1,N)
      WRITE(IOUT,30)
30    FORMAT ( //23X,17H SOLUTION VECTORS /)
      DO 35 J=NP1,LIM
35    WRITE(IOUT,7)(A(I,J),B(I,J),I=1,N)
      GO TO 8
C     HERE IF DET=0. PRINT OUT SOLUTION TO HOMOGENEOUS EQUATIONS.
40    KM1=K-1
      WRITE(IOUT,45)KM1
45    FORMAT (//37H INPUT MATRIX IS SINGULAR AND OF RANK I4 /
     1  31H SOLUTION TO HOMOGENEOUS SYSTEM /)
      DO 48 J=K,N
48    WRITE(IOUT,7)(A(I,J),B(I,J),I=1,N)
      GO TO 8
1     WRITE(IDLG,2)
2     FORMAT(1X,'DATA BEING PROCESSED'/)
      GO TO 4
      END
C--------------- A, B, M, N, E, JX ARE PASSED TO COMPNV.
C--------------- AID, BID, DETR, DETI, AND K ARE RETURNED BY COMPNV.
      SUBROUTINE COMPNV(A,B,AID,BID,M,N,E,DETR,DETI,K,JX,ICOL)
C     IF DET=0, THEN (K-1)= THE NUMBER OF HOMOGENEOUS EQUATIONS AND THE
C     SOLUTION IS LEFT IN COLUMNS K THRU N OF THE ORIGINAL A-MATRIX.
      DIMENSION A(1),B(1),AID(1),BID(1),ICOL(1)
      EQUIVALENCE(XLOC,XLOC1,XLOCI,TEMP,TEMPI,TEMP1),(ISUBI1,ISUB1,ISUB2
     1 ),(TERM,TERM1,CON1,TERM3,FACTR),(TERM2,CON2,TERM4,FACTI)
      E=ABS(E)
      LIM=N+M
      NP1=N+1
      NM1=N-1
      SIGN=1.
      DETR=1.
      DETI=0.
C     PRODUCE AN IDENTITY MATRIX IN ARRAY AID.
      DO 3 I=1,N
      DO 3 J=1,N
      ISUB=I+JX*(J-1)
      IF(I-J) 1,2,1
1     AID(ISUB)=0.
      GO TO 3
2     AID(ISUB)=1.
3     CONTINUE
      DO 400 I=1,N
      DO 400 J=1,N
      ISUB=I+JX*(J-1)
400   BID(ISUB)=0.
      DO 35 K=1,N
C---------------E=TOLERANCE ENTERED BY USER.  THIS IS CRITERION
C--------------- FOR DECIDING THAT A SMALL NO. IS REALLY ZERO.
      BIG=E
      DO 5 I=K,N
      DO 5 J=K,N
      ISUB=I+JX*(J-1)
      TERM=A(ISUB)*A(ISUB)+B(ISUB)*B(ISUB)
      IF(TERM-BIG) 5,5,4
4     IR=I
      IC=J
      BIG=TERM
5     CONTINUE
      IF(BIG-E) 6,6,11
C---------------IF BIG < E, DET. =0.
C---------------SEE FOLDER FOR MATH. DOCUM.
C     HERE IF DETERMINANT=0. PRODUCE SOLUTION TO HOMOGENEOUS SYSTEM.
6     DETR=0.
      DETI=0.
      KM1=K-1
      DO 7 I=1,KM1
      DO 7 J=K,N
      ISUB=I+JX*(J-1)
      A(ISUB)=-A(ISUB)
7     B(ISUB)=-B(ISUB)
      DO 10 I=K,N
      DO 10 J=K,N
      ISUB=I+JX*(J-1)
      IF(I-J) 8,9,8
8     A(ISUB)=0.
      GO TO 10
9     A(ISUB)=1.
10    CONTINUE
C     NOW RE-ORDER THE SOLUTIONS TO THE HOMOGENEOUS SYSTEM.
      I=K-1
301   IF(ICOL(I)-I) 302,305,302
302   I1=ICOL(I)
      DO 303 J=K,N
      ISUB=I+JX*(J-1)
      ISUBI1=I1+JX*(J-1)
      XLOC=A(ISUB)
      A(ISUB)=A(ISUBI1)
      A(ISUBI1)=XLOC
      XLOC1=B(ISUB)
      B(ISUB)=B(ISUBI1)
303   B(ISUBI1)=XLOC1
305   I=I-1
      IF(I-1) 306,301,301
306   RETURN
C---------------THE PRODUCT OF THE N PIVOTAL ELEMENTS GIVES
C--------------- THE DET.
C     MULT DETERMINANT BY PIVOTAL ELEMENT.
11    CON1=DETR
      CON2=DETI
      ISUB=IR+JX*(IC-1)
      CALL MULT(CON1,CON2,A(ISUB),B(ISUB),DETR,DETI)
C---------------IF LARGEST NO. IS IN THE KTH ROW, NO INTERCHANGE OF 
C--------------- ROWS IS NECESSARY.  INTERCHANGE IS DONE TO INCREASE
C--------------- ACCURACY.
      IF(IR-K) 12,20,12
C     HERE IF A ROW INTERCHG IS NECESSARY.
C---------------INTERCHANGE OF ROWS CHANGES SIGN OF DETERMINAMT.
12    SIGN=SIGN*(-1.)
      DO 15 J=K,LIM
      ISUB=IR+JX*(J-1)
      ISUB2=K+JX*(J-1)
      XLOC=A(ISUB)
      A(ISUB)=A(ISUB2)
      A(ISUB2)=XLOC
      XLOCI=B(ISUB)
      B(ISUB)=B(ISUB2)
15    B(ISUB2)=XLOCI
      DO 16 J=1,N
      ISUB=IR+JX*(J-1)
      ISUB2=K+JX*(J-1)
      XLOC1=AID(ISUB)
      AID(ISUB)=AID(ISUB2)
      AID(ISUB2)=XLOC1
      XLOCI=BID(ISUB)
      BID(ISUB)=BID(ISUB2)
16    BID(ISUB2)=XLOCI
20    ICOL(K)=IC
C---------------IF LARGEST NO. IS IN KTH COL., INTERCHANGE OF COLUMNS
C--------------- IS NECESSARY.
C---------------INTERCHANGE IS DONE TO INCREASE ACCURACY.
C     SEE IF COLUMN INTERCHG IS NECESSARY.
      IF(IC-K) 21,24,21
C---------------INTERCHANGE OF COLS. CHANGES SIGN OF DET.
21    SIGN=SIGN*(-1.)
      DO 23 I=1,N
      ISUB=I+JX*(K-1)
      ISUB2=I+JX*(IC-1)
      XLOC=A(ISUB)
      A(ISUB)=A(ISUB2)
      A(ISUB2)=XLOC
      XLOCI=B(ISUB)
      B(ISUB)=B(ISUB2)
23    B(ISUB2)=XLOCI
C--------------- BETWEEN HERE AND ST. #35 GAUSSIAN ELIM. IS BEING USED.
C     NOW REDUCE PIVOTAL ELEMENT TO 1.
24    ISUB=K+JX*(K-1)
      IF(A(ISUB)-1.) 5002,5001,5002
5001  IF(B(ISUB)) 5002,28,5002
5002  PIVR=A(ISUB)
      PIVI=B(ISUB)
      DO 25 J=K,LIM
      ISUB=K+JX*(J-1)
      TERM1=A(ISUB)
      TERM2=B(ISUB)
25    CALL DIV(TERM1,TERM2,PIVR,PIVI,A(ISUB),B(ISUB))
      DO 26 J=1,N
      ISUB=K+JX*(J-1)
      TERM3=AID(ISUB)
      TERM4=BID(ISUB)
26    CALL DIV(TERM3,TERM4,PIVR,PIVI,AID(ISUB),BID(ISUB))
C     NOW PRODUCE ZEROS IN COLUMN IC.
28    DO 34 I=1,N
      IF(I-K) 30,34,30
30    ISUB=I+JX*(K-1)
      IF(A(ISUB)) 1013,1012,1013
1012  IF(B(ISUB)) 1013,34,1013
1013  FACTR=A(ISUB)
      FACTI=B(ISUB)
      DO 32 J=K,LIM
      ISUB=K+JX*(J-1)
      ISUB2=I+JX*(J-1)
      CALL MULT(FACTR,FACTI,A(ISUB),B(ISUB),TERMR,TERMI)
      A(ISUB2)=A(ISUB2)-TERMR
32    B(ISUB2)=B(ISUB2)-TERMI
      DO 33 J=1,N
      ISUB=K+JX*(J-1)
      ISUB2=I+JX*(J-1)
      CALL MULT(FACTR,FACTI,AID(ISUB),BID(ISUB),TERMR,TERMI)
      AID(ISUB2)=AID(ISUB2)-TERMR
33    BID(ISUB2)=BID(ISUB2)-TERMI
34    CONTINUE
35    CONTINUE
C     NOW EVALUATE DETERMINANT.
C---------------SIGN REFLECTS THE RESULT OF INTERCHANGE OF ROWS AND
C--------------- INTERCHANGE OF COLS.
      DETR=DETR*SIGN
      DETI=DETI*SIGN
C     NOW RE-ORDER SOLUTION MATRIX
C---------------THIS RE-ORDERING GIVES SOLUTION IN THE ORDER WHICH
C--------------- THE USER SPECIFIED.  THEY MAY HAVE BEEN REARRANGED BY
C--------------- THE INTERCHANGE OF ROWS AND COLS. TO ACHIEVE HIGHER
C--------------- ACCURACY.
59    IF(ICOL(K)-K) 60,64,60
60    K1=ICOL(K)
      DO 63 J=NP1,LIM
      ISUB=K+JX*(J-1)
      ISUB1=K1+JX*(J-1)
      TEMP=A(ISUB)
      A(ISUB)=A(ISUB1)
      A(ISUB1)=TEMP
      TEMPI=B(ISUB)
      B(ISUB)=B(ISUB1)
63    B(ISUB1)=TEMPI
64    K=K-1
      IF(K-1) 70,59,59
C---------------THIS RE-ORDERING GIVES AN INVERSE WHICH WOULD HAVE BEEN
C--------------- OBTAINED WITHOUT AN INTERCHANGE OF ROWS AND COLUMNS.
C     NOW RE-ORDER ROWS OF INVERSE MATRIX.
70    K=N
71    IF(ICOL(K)-K) 73,83,73
73    K1=ICOL(K)
      DO 80 J=1,N
      ISUB=K+JX*(J-1)
      ISUB2=K1+JX*(J-1)
      TEMP1=AID(ISUB)
      AID(ISUB)=AID(ISUB2)
      AID(ISUB2)=TEMP1
      TEMPI=BID(ISUB)
      BID(ISUB)=BID(ISUB2)
80    BID(ISUB2)=TEMPI
83    K=K-1
      IF(K-1) 85,71,71
85    RETURN
      END
C--------------- A, B, C, D ARE INPUT.  E, F RETURNED.
      SUBROUTINE MULT(A,B,C,D,E,F)
      INP=5
      IOUT=30
      E=A*C-B*D
      F=A*D+B*C
      RETURN
      END
C--------------- A, B, C, D ARE INPUT.  E, F ARE RETURNED.
      SUBROUTINE DIV(A,B,C,D,E,F)
      DENOM=C*C+D*D
      E=(A*C+B*D)/DENOM
      F=(B*C-A*D)/DENOM
      RETURN
      END

