Trailing-Edge
-
PDP-10 Archives
-
decuslib10-09
-
43,50466/cminv.f4
There are no other files named cminv.f4 in the archive.
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