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