Trailing-Edge
-
PDP-10 Archives
-
decuslib20-05
-
decus/20-0137/matrix/matrix.for
There is 1 other file named matrix.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C MATRIX.F4 (FILENAME ON LIBRARY DECTAPE)
C MATRIX, 2.8.1. (CALLING NAME, SUBLST#)
C MATRIX OPERATIONS
C MATRIX WAS PROGRAMMED BY SAM ANEMA
C LIBRARY DECTAPE PROGRAMS USED: USAGE.MAC
C APLIB PROGS. USED: IO, GETFOR,
C INTERNAL SUBR: MINVSQ
C FORWMU PROGS. USED: DEVCHG, EXISTS, PRINTS
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
DIMENSION SUM(20),IFMT(32)
DATA BLANK/' '/
C---------------A(20,20) IS PLACE WHERE USER'S MATRIX IS STORED
C--------------- PRIOR TO CALL MINSQV
DIMENSION VAR(10),A(20,20)
DIMENSION MC(20),MR(20)
C---------------X(10,20,20) DETERMINES THE LIMITATION OF TEN
C--------------- 20X20 MATRICES.
IPT=0
DIMENSION X(10,20,20),NDIM(2,10),RED(70)
TYPE 2
2 FORMAT('1','---WMU MATRIX ALGEBRA PROGRAM'/12X,'(20X20)'//' FOR HE
1LP TYPE HELP')
C CALL USAGE('MATRIX')
1 TYPE 10
10 FORMAT(' *',$)
READ(5,11)(RED(I),I=1,70)
11 FORMAT(70A1)
K=1
DO 13 I=1,70
IF(RED(I).EQ.BLANK)GO TO 13
RED(K)=RED(I)
K=K+1
13 CONTINUE
DO 16 I=K,70
16 RED(I)=' '
IF(RED(1).EQ.'H'.AND.RED(2).EQ.'E'.AND.RED(3).EQ.'L'.AND.
1RED(4).EQ.'P')GO TO 1500
IF(RED(1).EQ.'R'.AND.RED(2).EQ.'E'.AND.RED(3).EQ.'A'.AND.
1RED(4).EQ.'D'.AND.RED(6).EQ.BLANK)GO TO 101
IF(RED(1).EQ.'T'.AND.RED(2).EQ.'Y'.AND.RED(3).EQ.'P'.AND.
1RED(4).EQ.'E'.AND.RED(6).EQ.BLANK)GO TO 201
IF(RED(1).EQ.'D'.AND.RED(2).EQ.'E'.AND.RED(3).EQ.'T'.AND.
1RED(5).EQ.BLANK)GO TO 301
IF(RED(2).EQ.'=')GO TO 401
IF(RED(1).EQ.'D'.AND.RED(2).EQ.'A'.AND.RED(3).EQ.'T'.AND.RED(4)
1.EQ.'A'.AND.RED(5).EQ.BLANK)GO TO 991
IF(RED(1).EQ.'T'.AND.RED(2).EQ.'R'.AND.RED(3).EQ.'A'.AND.RED(4).EQ
1.'C'.AND.RED(5).EQ.'E'.AND.RED(7).EQ.' ')GO TO 1011
TYPE 1002
1002 FORMAT(' ILLEGAL STATEMENT')
GO TO 1
C
C READ IN A MATRIX
C
101 IF(IPT.EQ.0)GO TO 106
DO 109 I2=1,IPT
IF(VAR(I2).EQ.RED(5))GO TO 108
109 CONTINUE
106 IPT=IPT+1
IF(IPT.GT.10)GO TO 994
VAR(IPT)=RED(5)
I2=IPT
108 TYPE 102
102 FORMAT(' DIMENSIONS--',$)
READ(5,103)(NDIM(I,I2),I=1,2)
103 FORMAT(2I)
IF(NDIM(1,I2).GT.20.OR.NDIM(2,I2).GT.20)GO TO 150
TYPE 104
104 FORMAT(' MATRIX?'/)
DO 107 J=1,NDIM(1,I2)
107 READ(5,105)(X(I2,J,K),K=1,NDIM(2,I2))
105 FORMAT(20F)
GO TO 1
150 TYPE 251
251 FORMAT(' DIMENSIONS MUST BE LESS THAN 20')
GO TO 1
C
C TYPE OUT A MATRIX
C
201 DO 202 I=1,IPT
IF(VAR(I).EQ.RED(5))GO TO 203
202 CONTINUE
TYPE 204,RED(5)
204 FORMAT(' MATRIX ',A1,' NOT FOUND')
GO TO 1
203 DO 207 J=1,NDIM(1,I)
TYPE 206,(X(I,J,K),K=1,NDIM(2,I))
206 FORMAT(' ',5F14.6)
207 TYPE 208
208 FORMAT(' ')
GO TO 1
C
C FIND THE DETERMINANT OF A MATRIX
C
301 DO 302 I=1,IPT
IF(VAR(I).EQ.RED(4))GO TO 303
302 CONTINUE
TYPE 204,RED(4)
GO TO 1
303 IF(NDIM(1,I).NE.NDIM(2,I))GO TO 350
DO 304 J=1,NDIM(1,I)
DO 305 K=1,NDIM(2,I)
A(J,K)=X(I,J,K)
305 CONTINUE
304 CONTINUE
CALL MINVSQ(A,NDIM(1,I),.1,MC,MR,20,5,0,DET,IEXP)
DET=DET*10.**IEXP
TYPE 310,RED(4),DET
310 FORMAT(' DETERMINANT OF ',A1,' IS'G)
GO TO 1
350 TYPE 351,RED(4)
351 FORMAT(' THE MATRIX IS NOT SQUARE')
GO TO 1
C
C SEARCH FOR OPERATOR
C
401 IF(RED(3).EQ.'T'.AND.RED(4).EQ.'R'.AND.RED(5).EQ.'A'.AND.
1RED(6).EQ.'N'.AND.RED(7).EQ.'S'.AND.RED(9).EQ.BLANK)GO TO 501
IF(RED(3).EQ.'I'.AND.RED(4).EQ.'N'.AND.RED(5).EQ.'V'.AND.
1RED(7).EQ.BLANK)GO TO 601
IF(RED(4).EQ.'+'.AND.RED(6).EQ.BLANK)GO TO 701
IF(RED(4).EQ.'*'.AND.RED(6).EQ.BLANK)GO TO 801
IF(RED(4).EQ.'-'.AND.RED(6).EQ.BLANK)GO TO 901
TYPE 1002
GO TO 1
C
C TRANSPOSE A MATRIX
C
501 DO 504 I2=1,IPT
IF(VAR(I2).EQ.RED(8))GO TO 510
504 CONTINUE
TYPE 204,RED(8)
GO TO 1
510 DO 502 I1=1,IPT
IF(VAR(I1).EQ.RED(1))GO TO 503
502 CONTINUE
520 IPT=IPT+1
IF(IPT.GT.10)GO TO 994
VAR(IPT)=RED(1)
I1=IPT
503 DO 530 I=1,NDIM(2,I2)
DO 531 J=1,NDIM(1,I2)
A(I,J)=X(I2,J,I)
531 CONTINUE
530 CONTINUE
DO 532 I=1,NDIM(2,I2)
DO 533 J=1,NDIM(1,I2)
X(I1,I,J)=A(I,J)
533 CONTINUE
532 CONTINUE
NDS=NDIM(2,I2)
NDIM(2,I1)=NDIM(1,I2)
NDIM(1,I1)=NDS
GO TO 1
IND=0
C
C INVERT A MATRIX
C
601 DO 602 I2=1,IPT
IF(VAR(I2).EQ.RED(6))GO TO 603
602 CONTINUE
TYPE 204,RED(6)
GO TO 1
603 DO 604 I1=1,IPT
IF(VAR(I1).EQ.RED(1))GO TO 605
604 CONTINUE
IPT=IPT+1
IF(IPT.GT.10)GO TO 994
IND=1
VAR(IPT)=RED(1)
I1=IPT
605 IF(NDIM(1,I2).NE.NDIM(2,I2))GO TO 608
DO 606 I=1,NDIM(1,I2)
DO 606 J=1,NDIM(2,I2)
606 A(I,J)=X(I2,I,J)
CALL MINVSQ(A,NDIM(1,I2),.1,MC,MR,20,5,0,DET,IEXP)
DO 607 I=1,NDIM(1,I2)
DO 607 J=1,NDIM(2,I2)
607 X(I1,I,J)=A(I,J)
NDIM(1,I1)=NDIM(1,I2)
NDIM(2,I1)=NDIM(2,I2)
GO TO 1
608 TYPE 351
IF(IND.NE.0)IPT=IPT-1
GO TO 1
C
C ADD OR SUBTRACT TWO MATRICES
C
C---------------COME HERE FROM ST. 401+4
701 IND=0
FM=1
720 DO 702 I2=1,IPT
IF(VAR(I2).EQ.RED(3))GO TO 703
702 CONTINUE
TYPE 204,RED(3)
GO TO 1
703 DO 704 I3=1,IPT
IF(VAR(I3).EQ.RED(5))GO TO 705
704 CONTINUE
TYPE 204,RED(5)
GO TO 1
705 IF(NDIM(1,I2).NE.NDIM(1,I3).OR.NDIM(2,I2).NE.NDIM(2,I3))GO TO 707
DO 710I1=1,IPT
IF(VAR(I1).EQ.RED(1))GO TO 712
710 CONTINUE
IPT=IPT+1
IF(IPT.GT.10)GO TO 994
IND=1
VAR(IPT)=RED(1)
I1=IPT
712 DO 713 I=1,NDIM(1,I2)
DO 713 J=1,NDIM(2,I2)
713 X(I1,I,J)=X(I2,I,J)+FM*X(I3,I,J)
NDIM(1,I1)=NDIM(1,I2)
NDIM(2,I1)=NDIM(2,I2)
GO TO 1
707 TYPE 708,RED(3),RED(5)
708 FORMAT(' DIMENSIONS OF ',A1,' AND ',A1,' ARE NOT THE SAME')
GO TO 1
C
C SUBTRACT
C
C---------------COME HERE FROM ST. 401+6
901 FM=-1
IND=0
GO TO 720
C
C FIND THE PRODUCT OF TWO MATRICES
C
801 IND=0
DO 802 I2=1,IPT
IF(VAR(I2).EQ.RED(3))GO TO 803
802 CONTINUE
TYPE 204,RED(3)
GO TO 1
803 DO 804 I3=1,IPT
IF(VAR(I3).EQ.RED(5))GO TO 805
804 CONTINUE
TYPE 204 ,RED(5)
GO TO 1
805 IF(NDIM(2,I2).NE.NDIM(1,I3))GO TO 850
DO 806 I1=1,IPT
IF(VAR(I1).EQ.RED(1))GO TO 807
806 CONTINUE
IPT=IPT+1
IF(IPT.GT.10)GO TO 994
IND=1
VAR(IPT)=RED(1)
I1=IPT
807 ND1=NDIM(1,I2)
ND2=NDIM(2,I3)
DO 820 I=1,ND1
DO 820 J=1,ND2
A(I,J)=0.0
DO 820 K=1,NDIM(2,I2)
820 A(I,J)=A(I,J)+X(I2,I,K)*X(I3,K,J)
DO 821 I=1,ND1
DO 822 J=1,ND2
X(I1,I,J)=A(I,J)
822 CONTINUE
821 CONTINUE
NDIM(1,I1)=ND1
NDIM(2,I1)=ND2
GO TO 1
850 TYPE 851,RED(3),RED(5)
851 FORMAT(' ILLEGAL DIMENSIONS ON ',A1,' OR ',A1)
GO TO 1
C
C DATA ENTRY OPTION
C
991 CALL IO(2,IDV,-1,-4,0,0)
CALL GETFOR(-1,-4,IFMT,ISTD,32,2)
IF(ISTD.EQ.1)IFMT(1)='(10F)'
TYPE 902
902 FORMAT(' ENTER TYPE OF MATRIX'/3X,
1'1 - RAW SUMS OF SQUARES AND CROSS PRODUCTS'/3X,
2'2 - STANDARDIZED SUMS OF SQUARES AND CROSS PRODUCTS'/3X,
3'3 - COVARIANCES'/3X,
4'4 - CORRELATIONS'/)
READ(5,903)NTO
903 FORMAT(I)
998 TYPE 904
904 FORMAT(' ENTER NUMBER OF VARIABLES'/)
READ(5,903)NVAR
IF(NVAR.GT.20.OR.NVAR.LT.1) GO TO 998
TYPE 905
905 FORMAT(' ENTER MATRIX NAME'/)
READ(5,906)ANA
906 FORMAT(A1)
DO 907 I2=1,IPT
IF(VAR(I2).EQ.ANA)GO TO 908
907 CONTINUE
IPT=IPT+1
IF(IPT.GT.10)GO TO 994
I2=IPT
908 VAR(I2)=ANA
NDIM(1,I2)=NVAR
NDIM(2,I2)=NVAR
K=0
DO 912 I=1,NVAR
SUM(I)=0.0
DO 912 J=1,NVAR
912 X(I2,I,J)=0.0
IF(IDV.EQ.'TTY')TYPE 1017
1017 FORMAT(' ENTER DATA'/)
IF(IDV.NE.'TTY') TYPE 1019
1019 FORMAT(' YOUR DATA IS BEING READ'/)
IF(ISTD.EQ.1)TYPE 985
985 FORMAT(' (10 PER LINE)'/)
911 READ(2,IFMT,END=950)(RED(I),I=1,NVAR)
910 FORMAT(10F)
K=K+1
DO 914 I=1,NVAR
SUM(I)=SUM(I)+RED(I)
DO 914 J=1,NVAR
914 X(I2,I,J)=X(I2,I,J)+RED(I)*RED(J)
GO TO 911
950 GO TO (1,952,952,952),NTO
GO TO 992
952 DO 956 I=1,NVAR
DO 956 J=1,NVAR
956 X(I2,I,J)=X(I2,I,J)-(SUM(I)*SUM(J)/FLOAT(K))
GO TO (1,1,953,953),NTO
953 DO 957 I=1,NVAR
DO 957 J=1,NVAR
957 X(I2,I,J)=X(I2,I,J)/FLOAT(K-1)
GO TO (1,1,1,954),NTO
954 DO 958 I=1,NVAR
958 SUM(I)=SQRT(X(I2,I,I))
DO 959 I=1,NVAR
DO 959 J=1,NVAR
959 X(I2,I,J)=X(I2,I,J)/(SUM(I)*SUM(J))
GO TO 1
1011 DO 1012 I=1,IPT
IF(VAR(I).EQ.RED(6))GO TO 1013
1012 CONTINUE
TYPE 204,RED(6)
GO TO 1
1013 IF(NDIM(1,I).NE.NDIM(2,I))GO TO 350
SUMM=0.0
DO 1014J=1,NDIM(1,I)
1014 SUMM=SUMM+X(I,J,J)
TYPE 1015,VAR(I),SUMM
1015 FORMAT(' TRACE OF ',A1,' IS ',G)
GO TO 1
992 TYPE 993
993 FORMAT(' NO SUCH OPTION')
GO TO 1
994 TYPE 995
995 FORMAT(' MORE THAN 10 MATRICES IN USE!'/)
IPT=10
GO TO 1
1500 TYPE 1501
1501 FORMAT(//7X,'COMMANDS:'//' READ A -- READ IN A MATRIX AND CALL I
1T A'/' TYPE A -- TYPE OUT MATRIX A'/' DET A --- CALCULATE THE
1DETERMINANT OF A'/' TRACE A - TYPE OUT THE TRACE OF SQUARE MATRI
1X A'/' B=TRANS A LET MATRIX B BE THE TRANSPOSE OF M
1ATRIX A'/' B=INV A - LET MATRIX B BE THE INVERSE OF MATRIX A'/'
1C=A+B --- LET MATRIX C BE THE SUM OF MATRICES A AND B'/' C=A-B
1--- LET MATRIX C BE THE DIFFERENCE OF MATRICES A AND B'/' C=A*B
1--- LET MATRIX C BE THE PRODUCT OF MATRICES A AND B'/' DATA ----
1 ENTER DATA FROM WHICH A MATRIX WILL BE GENERATED'/' HELP ----
1TYPE OUT THIS COMMAND LIST'//)
GO TO 1
END
C---------------A, N, NDIM, METHOD, IOUT, TOL ARE INPUT.
C--------------- MR, MC, DET, IEXP ARE OUTPUT. A IS MODIFIED.
SUBROUTINE MINVSQ(A,N,TOL,MC,MR,NDIM,IOUT,METHOD,DET,IEXP)
C THIS SUBROUTINE INVERTS A SQUARE MATRIX WITHIN ITSELF.
C IT IS A MODIFICATION OF MINV IN THE IBM SSP MANUAL.
C A...IS MATRIX TO BE INVERTED.
C DIMENSION FOR A IN MAINLINE MUST BE NDIM BY NDIM
C MC,MR...ARE VECTORS REQUIRED BY THE SUBROUTINE FOR BOOKEEPING.
C DIMENSION FOR MC,MR IN MAINLINE MUST BE NDIM (OR LARGER)
C N...IS NUMBER OF ROWS IN MATRIX TO BE INVERTED. (N.LE.NDIM)
C TOL...SHOULD BE SIMILAR IN MAGNITUDE TO SMALLEST NONZERO ELEMENT
C IN MATRIX
C DET...IS DETERMINANT OF MATRIX (IF DET RETURNED EQUALS ZERO, THE
C INVERSE MATRIX DID NOT EXIST).
C (DET IS STORED IN SCIENTIFIC NOTATION, WITH THE POWER OF 10
C STORED IN IEXP. THIS WAS NECESSITATED BY FREQUENT OVERFLOWS
C ON THE EXPONENT)
C METHOD...IS A SWITCH TO SELECT PIVOTING TECHNIQUE.
C METHOD=0 FASTEST TECHNIQUE, LEAST ACCURATE
C METHOD=1 SLOWER THAN 0, BUT MORE ACCURATE
C METHOD=2 SLOWEST BUT BEST ACCURACY
C IOUT...IS OUTPUT DEVICE SPECIFICATION
C
DIMENSION A(1),MC(1),MR(1)
ZTOL=TOL*1.E-6
DET=1.
IEXP=0
DO 1 I=1,N
MR(I)=I
1 MC(I)=I
KSGN=1
KK=-NDIM
7 DO 100 K=1,N
KK=KK+NDIM
NL=N
IF (METHOD-1) 4,3,2
3 NL=K
2 AMAX=0.
JJ=KK-NDIM
DO 6 J=K,N
JJ=JJ+NDIM
DO 5 I=K,NL
IJ=I+JJ
IF (ABS(A(IJ)).LE.AMAX)GO TO 5
AMAX=ABS(A(IJ))
NR=I
NC=J
5 CONTINUE
6 CONTINUE
GO TO 10
4 II=KK-NDIM
DO 62 I=K,N
II=II+NDIM
KI=II+K
IF (ABS(A(KI)).GT.ZTOL) GO TO 11
62 CONTINUE
13 WRITE(IOUT,51)
51 FORMAT(1X,'INVERSE DOES NOT EXIST')
DET=0.
IEXP=0
RETURN
11 AMAX=A(KI)
NR=K
NC=I
10 IF (ABS(AMAX).LE.ZTOL) GO TO 13
IF (NR.EQ.K) GO TO 9
KSGN=-KSGN
MR(K)=NR
JJ=-NDIM
DO 12 J=1,N
JJ=JJ+NDIM
KJ=K+JJ
NRJ=NR+JJ
B=A(KJ)
A(KJ)=A(NRJ)
12 A(NRJ)=B
9 IF (NC.EQ.K) GO TO 22
KSGN=-KSGN
MC(K)=NC
NCNC=(NC-1)*NDIM
DO 20 J=1,N
JK=J+KK
JNC=J+NCNC
B=A(JK)
A(JK)=A(JNC)
20 A(JNC)=B
22 KKK=K+KK
D=A(KKK)
DET=DET*D
205 IF (ABS(DET).LT.10.) GO TO 200
DET=DET/10.
IEXP=IEXP+1
GO TO 205
200 IF (ABS(DET).GE.1.) GO TO 210
DET=DET*10.
IEXP=IEXP-1
IF (IEXP-40) 200,13,13
210 IF (DET.EQ.0.) GO TO 13
DO 30 I=1,N
IK=I+KK
30 A(IK)=A(IK)/D
A(KKK)=1./D
II=-NDIM
DO 40 I=1,N
II=II+NDIM
KI=K+II
C=A(KI)
IF (I-K) 42,40,42
42 IF (C) 41,40,41
41 DO 50 J=1,N
JI=J+II
JK=J+KK
A(JI)=A(JI)-C*A(JK)
50 CONTINUE
A(KI)=-C/D
40 CONTINUE
100 CONTINUE
DET=DET*FLOAT(KSGN)
DO 155 K=N,1,-1
L=MC(K)
IF (K-L) 150,155,150
150 DO 170 I=N,1,-1
II=(I-1)*NDIM
LI=L+II
KI=K+II
B=A(LI)
A(LI)=A(KI)
170 A(KI)=B
155 CONTINUE
DO 175 K=N,1,-1
L=MR(K)
IF (K-L) 180,175,180
180 LL=(L-1)*NDIM
KK=(K-1)*NDIM
DO 185 I=1,N
IL=I+LL
IK=I+KK
B=A(IL)
A(IL)=A(IK)
185 A(IK)=B
175 CONTINUE
RETURN
END