C cancor.f
C The authors of this software are Jih-Jie Chang and J D Carroll.

C Copyright (c) 1993 by AT&T.
C Permission to use, copy, modify, and distribute this software for any
C purpose without fee is hereby granted, provided that this entire
C notice is included in all copies of any software which is or includes
C a copy or modification of this software and in all copies of the
C supporting documentation for such software.
C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP-
C LIED WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

C This software comes from the FIRST MDS Package of AT&T Bell Laboratories

C For explanation of the method this software implements, see
C Carroll, J.D. (1968), "Generalization of canonical correlation
C analysis to three or more sets of variables" in Proceedings of the
C 76th Annual Convention of the American Psychological Association, 3,
C 227-228.
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

$      FORTRAN
$      INCODE  IBMF
C     CANCOR - CANONICAL CORRELATION ANALYSIS
C     OF THREE OR MORE SETS OF VARIABLES
C     CHANG FOR CARROLL,  FEBRUARY 1969
C     NSET - NUMBER OF SETS (MAX=100)
C     NS - NUMBER OF STIMULI (MAX=100).
C     NCS - NUMBER OF COMPONENTS WITHIN EACH SET.
C     NC - NUMBER OF COMPONENTS TO BE COMPUTED
C     WF - WEIGHTING FACTOR, ONE FOR EACH SET.
C     IRX - =1  THE SET IS PUNCHED NCS BY NS
C           =0  THE SET IS PUNCHED NS BY NCS
C     ICOR - =1  COMPUTE CORRELATIONS BETWEEN SETS
C            =0  DO NOT COMPUTE CORRELATIONS BETWEEN SETS
C     X - DATA, STORES ONE SET OF VARIABLE. SCRATCH DISK IS USED TO
C     STORE ALL SETS.
C     TIT - TITLE
C     SNAME - SUB-TITLE, ONE FOR EACH SET.
C     FMT - VARIABLE FORMAT FOR READING DATA
C     CORL - CORRELATIONS BETWEEN EACH SET AND THE CANONICAL COMPONENTS
      DIMENSION X(100,10),Q(100,100),Z(100,100),Y(100,10),C(100,10)
      DIMENSION NCS(100),SNAME(12,100), FMT(12),TIT(12)
      DIMENSION CORL(10,100),SUMZ(10),SSQZ(10)
      EQUIVALENCE(C,Z)
      ND1=100
      ND2=10
      NDD=1000
      NDH=11
 500  READ100,NSET,NS,NC,ICOR
      IF(NSET.EQ.0) GOTO999
      XNS=NS
      READ101,(TIT(I),I=1,12)
      PRINT102,(TIT(I),I=1,12)
      DO5I=1,100
      DO5J=1,100
 5    Q(I,J)=0.
      PRINT113
 113  FORMAT(13H0I.  DATA SET)
      DO10K=1,NSET
      READ101,(SNAME(I,K),I=1,12)
      READ112,NCS(K),IRX,WF
      M=NCS(K)
      CALL READX(X,NS,M,IRX,ND1)
C     X SUBTRACT XMEAN
      DO15I=1,M
      SUM=0.
      DO16J=1,NS
 16   SUM=SUM&X(J,I)
      XMEAN=SUM/XNS
      DO17J=1,NS
 17   X(J,I)=X(J,I)-XMEAN
 15   CONTINUE
C     PUT X ON DISK 19
      IF(K.EQ.1) GOTO11
      CALL DOUC(19,X,NDD)
      GOTO13
 11   CALL DOUTAT(19,X,NDD,1)
 13   PRINT103,K,(SNAME(I,K),I=1,12)
      PRINT104
C     PRINT X
      DO18I=1,M
 18   PRINT105,I,(X(J,I),J=1,NS)
C     COMPUTE (XT)X AND STORE IN Y
      DO19I=1,M
      DO19L=1,M
      Y(I,L)=0.
      DO19J=1,NS
 19   Y(I,L)=Y(I,L)&X(J,I)*X(J,L)
C     INVERT (XT)X AND STORE IN Y
      MM=1
      CALL MTINV(Y,ND1,ND2,C,NDH,NDH,M,MM)
      IF(MM.EQ.1)GOTO20
      PRINT106,MM,K
      STOP
C     COMPUTE X*((XT)X) INVERSE
 20   DO21I=1,NS
      DO21J=1,M
      C(I,J)=0.
      DO21L=1,M
 21   C(I,J)=C(I,J)&X(I,L)*Y(L,J)
C     PUT C ON DISK 19
      CALL DOUC(19,C,NDD)
C     COMPUTE Q MATRIX
      DO22I=1,NS
      DO22J=1,NS
      SUM=0.
      DO23L=1,M
 23   SUM=SUM&C(I,L)*X(J,L)
 22   Q(I,J)=Q(I,J)&SUM*WF
 10   CONTINUE
      CALL DROUT
C     FACTOR Q MATRIX
      MM=1
      CALL MTEVV(Q,ND1,ND1,Z,ND1,ND1,NS,MM)
      PRINT107
      PRINT108,(Q(I,I),I=1,NS)
      PRINT109
      DO24I=1,NC
C     SET UP FOR CORRELATIONS BETWEEN Z AND EACH SET.
      SUMZ(I)=0.
      SSQZ(I)=0.
      DO25J=1,NS
      SUMZ(I)=SUMZ(I)&Z(J,I)
 25   SSQZ(I)=SSQZ(I)&Z(J,I)**2
 24   PRINT105,I,(Z(J,I),J=1,NS)
C     COMPUTE A(I), I FOR I(TH) SET, STORE A IN Q
      PRINT114
      DO30K=1,NSET
      M=NCS(K)
C     READ X AND C FROM DISK. PUT C IN Y
      IF(K.GT.1)GOTO31
      CALL DINAT(19,X,NDD,1)
      GOTO32
 31   CALL DINC(19,X,NDD)
 32   CALL DINC(19,Y,NDD)
C     COMPUTE A AND STORE A IN Q
      DO35I=1,NC
      SUM=0.
      DO37J=1,M
      Q(I,J)=0.
      DO36L=1,NS
 36   Q(I,J)=Q(I,J)&Z(L,I)*Y(L,J)
 37   SUM=SUM&Q(I,J)**2
      SUM=SQRT(SUM)
      DO38LL=1,M
 38   Q(I,LL)=Q(I,LL)/SUM
 35   CONTINUE
C     PRINT CANONICAL DIRECTIONS
      PRINT103,K,(SNAME(I,K),I=1,12)
      PRINT110
      DO39I=1,NC
 39   PRINT105,I,(Q(I,J),J=1,M)
C     COMPUTE ROTATION MATRIX AND STORE IN Y
      DO40I=1,NC
      SXY=0.
      SUMY=0.
      SSQY=0.
      DO47J=1,NS
      Y(J,I)=0.
      DO45L=1,M
 45   Y(J,I)=Y(J,I)&Q(I,L)*X(J,L)
C     COMPUTE CORRELATIONS BETWEEN Z AND EACH SET OF CANONICAL VARIATE
      SUMY=SUMY&Y(J,I)
      SSQY=SSQY&Y(J,I)**2
 47   SXY=SXY&Y(J,I)*Z(J,I)
      PP=XNS*SXY-SUMY*SUMZ(I)
      QQ=(XNS*SSQY-SUMY**2)*(XNS*SSQZ(I)-SUMZ(I)**2)
      QQ=SQRT(QQ)
      CORL(I,K)=PP/QQ
 40   CONTINUE
      IF(ICOR.EQ.0) GOTO43
C     PUT Y ON DISK 20
      IF(K.EQ.1) GOTO42
      CALL DOUC(20,Y,NDD)
      GOTO43
 42   CALL DOUTAT(20,Y,NDD,1)
C     PRINT Y
 43   PRINT111
      DO41I=1,NC
 41   PRINT105,I,(Y(J,I),J=1,NS)
 30   CONTINUE
      PRINT115
      PRINT116
      DO46I=1,NSET
 46   PRINT117,I,(CORL(J,I),J=1,NC)
      IF(ICOR)44,500,44
 44   CALL DROUT
C     COMPUTE CORRELATIONS BETWEEN SETS FOR EACH COMPONENTS
      CALL  CCOR(Y,NC,NS,NDD,X,Q,ND1,NSET,XNS)
      GOTO500
 100  FORMAT(18I4)
 101  FORMAT(12A6)
 102  FORMAT(1H112A6)
 103  FORMAT(4H0SETI3,3X,12A6)
 104  FORMAT(5X,26HDATA (SUBTRACTED ROW MEAN))
 105  FORMAT(1H04X,I2,2X,10E12.4/(9X,10E12.4))
 106  FORMAT(4H0MM=I1,32HMATRIX CANNOT BE INVERTED IN SETI3)
 107  FORMAT(41H0II.  CANONICAL VALUES (EIGEN ROOTS OF Q))
 108  FORMAT(5X,10E12.4)
 109  FORMAT(47H0     CANONICAL COMPONENTS (EIGEN VECTORS OF Q))
 110  FORMAT(23H0  CANONICAL DIRECTIONS)
 111  FORMAT(28H0  CANONICAL ROTATION MATRIX)
 112  FORMAT(2I2,F10.5)
 114  FORMAT(49H0III.  CANONICAL DIRECTIONS AND ROTATION MATRICES)
 115  FORMAT(89H0IV.  CORRELATIONS BETWEEN EACH SET OF CANONICAL VARIATE
     1 AND THE CANONICAL COMPONENTS (Z))
 116  FORMAT(5H0 SET,14X,10HCOMPONENTS/6X,I9,9I12)
 117  FORMAT(I4,5X,10E12.4)
 999  STOP
      END
$      FORTRAN
      SUBROUTINE  CCOR(Y,NC,NS,NDD,X,Q,ND1,NSET,XNS)
      DIMENSION Y(ND1,NC),X(ND1,NC),Q(ND1,ND1)
      DIMENSION SUMY(10),SSQY(10)
      N=NSET-1
      IC=0
      NY=1
      DO10I=1,N
      CALL DINAT(20,Y,NDD,NY)
      NY=NY&NDD
      DO15J=1,NC
      SUMY(J)=0.
      SSQY(J)=0.
      DO15L=1,NS
      SUMY(J)=SUMY(J)&Y(L,J)
 15   SSQY(J)=SSQY(J)&Y(L,J)**2
      II=I&1
      DO20JJ=II,NSET
      IC=IC&1
      CALL DINC(20,X,NDD)
      DO20J=1,NC
      SUMX=0.
      SSQX=0.
      SXY=0.
      DO21L=1,NS
      SUMX=SUMX&X(L,J)
      SSQX=SSQX&X(L,J)**2
 21   SXY=SXY&X(L,J)*Y(L,J)
      PP=XNS*SXY-SUMY(J)*SUMX
      QQ=(XNS*SSQY(J)-SUMY(J)**2)*(XNS*SSQX-SUMX**2)
      QQ=SQRT(QQ)
 20   Q(J,IC)=PP/QQ
 10   CONTINUE
      PRINT100
 100  FORMAT(50H0V.   CORRELATIONS BETWEEN SETS FOR EACH COMPONENT)
      PRINT101,(I,I=1,NC)
 101  FORMAT(6H0 SETS13X,10HCOMPONENTS/6H  I  J,I9,9I12)
      IC=0
      DO30I=1,N
      II=I&1
      DO30J=II,NSET
      IC=IC&1
 30   PRINT102,I,J,(Q(L,IC),L=1,NC)
 102  FORMAT(2I3,3X,10E12.4)
      RETURN
      END
$      FORTRAN     READX(X,N,K,IRX,NL)
      SUBROUTINE READX(X,N,K,IRX,NL)
      DIMENSION X(NL,K),FMT(12)
      READ100,(FMT(I),I=1,12)
 100  FORMAT(12A6)
      IF(IRX.EQ.0) GOTO10
      DO6I=1,K
 6    READFMT,(X(J,I),J=1,N)
      GOTO20
 10   DO8I=1,N
 8    READFMT,(X(I,J),J=1,K)
 20   RETURN
      END
$      FORTRAN                                                          EIGJ00a0
$      INCODE  IBMF                                                     EIGJ00b0
*EIGENJ      EIGENVALUE AND EIGENVECTOR SUBROUTINE                      EIGJ0020
*                              CD600D4.009      DATE 05/05/65           EIGJ0030
      SUBROUTINE EIGENJ (A,V,NN,IV,E,I3,I4)                             EIGJ00e0
      DIMENSION A(5050),V( I3, I4),E(NN)                                EIGJ00f0
      LOGICAL IV                                                        EIGJ0060
*     EIGENJ FINDS THE EIGENVALUES AND EIGENVECTORS OF                  EIGJ0070
*     A SYMMETRIC MATRIX (A) USING A MODIFIED THRESHOLD JACOBI METHOD   EIGJ0080
*     ONLY THE UPPER TRIANGLE IS STORED                                 EIGJ0090
*     DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER         EIGJ0100
*     DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER         EIGJ0110
      N=NN                                                              EIGJ0120
      IND=1                                                             EIGJ0130
      NM2=N-2                                                           EIGJ0140
      AMAX=0.0                                                          EIGJ0150
      NM=N-1                                                            EIGJ0160
      IF(.NOT.IV)GO TO 5                                                EIGJ0170
*     SET UP VECTOR MATRIX                                              EIGJ0180
    1 DO 3 I=1,N                                                        EIGJ0190
      DO 2 J=1,N                                                        EIGJ0200
    2 V(I,J)=0.0                                                        EIGJ0210
    3 V(I,I)=1.0                                                        EIGJ0220
    5 K=2                                                               EIGJ0230
      DO 28 I=1,NM                                                      EIGJ0240
      IP=I&1                                                            EIGJ0250
      DO 27 J=IP,N                                                      EIGJ0260
      Y=A(K)                                                            EIGJ0270
      X=ABS(Y)                                                          EIGJ0280
      IF(X-AMAX)27,27,6                                                 EIGJ0290
    6 GO TO (7,8),IND                                                   EIGJ0300
    7 AMAX=X                                                            EIGJ0310
      ITEST=1                                                           EIGJ0360
      X=A(II)-A(JJ)                                                     EIGJ0370
      IF(X)10,11,10                                                     EIGJ0380
   10 X=Y/X                                                             EIGJ0390
      Y=.5*X                                                            EIGJ0400
      IF(ABS(Y)-.414213562)14,11,11                                     EIGJ0410
   11 C=.707106781                                                      EIGJ0420
      IF(Y)12,12,13                                                     EIGJ0430
   12 S=-C                                                              EIGJ0440
      GO TO 15                                                          EIGJ0450
   13 S=C                                                               EIGJ0460
      GO TO 15                                                          EIGJ0470
   14 X=Y*Y                                                             EIGJ0480
      C=1.&X                                                            EIGJ0490
      S=2.*Y/C                                                          EIGJ0500
      C=(1.-X)/C                                                        EIGJ0510
   15 X=S*S                                                             EIGJ0520
      Y=C*C                                                             EIGJ0530
      XY=S*C                                                            EIGJ0540
      AXY=2.*A(K)*XY                                                    EIGJ0550
      Z=A(II)*Y&AXY&A(JJ)*X                                             EIGJ0560
      W=A(II)*X-AXY&A(JJ)*Y                                             EIGJ0570
      A(K)=A(K)*(Y-X)&XY*(A(JJ)-A(II))                                  EIGJ0580
      A(II)=Z                                                           EIGJ0590
      A(JJ)=W                                                           EIGJ0600
      IF(NM2)24,24,16                                                   EIGJ0610
   16 IT=I                                                              EIGJ0620
      JT=J                                                              EIGJ0630
*     TRANSFORM A                                                       EIGJ0640
      DO 23 M=1,NM2                                                     EIGJ0650
      NP=N-M                                                            EIGJ0660
      IF(JT-K)20,17,18                                                  EIGJ0670
   17 IT=IT&1                                                           EIGJ0680
      JT=JT&NP                                                          EIGJ0690
      ITEST=2                                                           EIGJ0700
   18 IF(IT-K)20,19,19                                                  EIGJ0710
   19 IT=IT&1                                                           EIGJ0720
      JT=JT&1                                                           EIGJ0730
      ITEST=3                                                           EIGJ0740
   20 X=A(IT)*C&A(JT)*S                                                 EIGJ0750
      A(JT)=A(JT)*C-A(IT)*S                                             EIGJ0760
      A(IT)=X                                                           EIGJ0770
      GO TO (21,22,23),ITEST                                            EIGJ0780
   21 IT=IT&NP                                                          EIGJ0790
      JT=JT&NP                                                          EIGJ0800
      GO TO 23                                                          EIGJ0810
   22 IT=IT&1                                                           EIGJ0820
      JT=JT&NP-1                                                        EIGJ0830
   23 CONTINUE                                                          EIGJ0840
   24 IF(.NOT.IV)GO TO 27                                               EIGJ0850
*     TRANSFORM V                                                       EIGJ0860
   25 DO 26 M=1,N                                                       EIGJ0870
      X=V(M,I)*C&V(M,J)*S                                               EIGJ0880
      V(M,J)=V(M,J)*C-V(M,I)*S                                          EIGJ0890
   26 V(M,I)=X                                                          EIGJ0900
   27 K=K&1                                                             EIGJ0910
   28 K=K&1                                                             EIGJ0920
*     SEQUENCE TO ADJUST THRESHOLD FOR EACH SWEEP                       EIGJ0930
      GO TO (29,31),IND                                                 EIGJ0940
   29 IF(AMAX)30,35,30                                                  EIGJ0950
   30 AT=AMAX                                                           EIGJ0960
      IND=2                                                             EIGJ0970
      F=.25                                                             EIGJ0980
      AMAX=F*AT                                                         EIGJ0990
      MAX=6                                                             EIGJ1000
      GO TO 5                                                           EIGJ1010
   31 MAX=MAX-1                                                         EIGJ1020
      IF(MAX)34,33,32                                                   EIGJ1030
   32 F=2.0*F*F                                                         EIGJ1040
      AMAX=AT*F                                                         EIGJ1050
   33 C=0.0                                                             EIGJ1060
      GO TO 5                                                           EIGJ1070
   34 IF(C)33,35,33                                                     EIGJ1080
*     COPY EIGENVALUES INTO E                                           EIGJ1090
   35 MM=1                                                              EIGJ1100
      NM=N&1                                                            EIGJ1110
      DO 36 M=1,N                                                       EIGJ1120
      E(M)=A(MM)                                                        EIGJ1130
   36 MM=MM&NM-M                                                        EIGJ1140
      RETURN                                                            EIGJ1150
      END                                                               EIGJ11g0
$      FORTRAN
$      INCODE  IBMF
      SUBROUTINE MTEVV(X,I1,I2,D,I3,I4,N,M)
      DIMENSION A(5050),X(I1,I2),D(I3,I4),E(100),NS(100),MS(100)
      K=1
      DO1J=1,N
      DO1I=J,N
      A(K)=X(I,J)
 1    K=K&1
      CALL EIGENJ(A,D,N,M,E,I3,I4)
C     THE FOLLOWING SYSTEM SORT PACKAGE ROUTINES SORT THE FLOATING POINT
C     ARRAY E OF N ELEMENTS IN DESCENDING ORDER.
      CALL SORTFL(-1)
      IF(M)4,3,4
 3    CALL SORT(E(1),E(2),N)
      DO 50 I=1,N
 50   X(I,I)=E(I)
      RETURN
 4    DO 5 I=1,N
       MS(I)=I
 5    NS(I)=I
C     THE FOLLOWING ROUTINE SORTS THE FLOATING POINT ARRAY E OF N
C     ELEMENTS IN DESCENDING ORDER AND ARRAY MS IS THE PASSIVE LIST
C     WHICH WILL BE REARRANGED IN THE SAME ORDER AS E.
      CALL SORT(E(1),E(2),N,MS(1))
C     THE FOLLOWING ROUTINES SORT THE FIXED POINT ARRAY MS OF N ELEMENTS
C     IN ASCENDING ORDER WITH ARRAY NS AS THE PASSIVE LIST.
       CALL SORTFX(1)
       CALL SORT(MS(1),MS(2),N,NS(1))
      CALL ARREV(N,NS,D,I3,I4)
      DO 6 I=1,N
 6    X(I,I)=E(I)
      RETURN
      END
$      FORTRAN
$      INCODE  IBMF
      SUBROUTINE ARREV(N,NS,D,I3,I4)
      DIMENSION NS(N),D(I3,I4)
      DO 10 I=1,N
 1    IF(NS(I)-I)2,10,2
 2    CALL XCH(N,I,NS(I),D,NS,I3,I4)
      GOTO1
 10   CONTINUE
      RETURN
      END
$      FORTRAN
$      INCODE  IBMF
      SUBROUTINE XCH(N,L,M,D,NS,I3,I4)
      DIMENSION        NS(N),D(I3,I4)
      I=L
       J=M
      DO 1 K=1,N
      T=D(K,I)
      D(K,I)=D(K,J)
 1    D(K,J)=T
      NS(I)=NS(J)
      NS(J)=J
      RETURN
      END
$      FORTRAN
$      INCODE  IBMF
CMTINV             MTINV
      SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M)                            MTPK1115
      DIMENSION A(I1,I2), B(I3,I4)                                      MTPK1116
C                                                                       MTPK1117
C     THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS-  MTPK1118
C     JORDAN ELIMINATION SCHEME.  ALL CALCULATIONS ARE DONE IN DOUBLE-  MTPK1119
C     PRECISION ARITHMETIC                                              MTPK1120
C                                                                       MTPK1121
      DOUBLE PRECISION B,BMAX,BMULT                                     MTPK1122
      EQUIVALENCE (BMAX,BMULT)                                          MTPK1123
C                                                                       MTPK1124
C     SINGLE TO DOUBLE PRECISION TRANSFER                               MTPK1125
C                                                                       MTPK1126
  101 DO 103 I = 1,N                                                    MTPK1127
  102 DO 103 J = 1,N                                                    MTPK1128
  103 B(I,J) = A(I,J)                                                   MTPK1129
C                                                                       MTPK1130
C     FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMATRIX. MTPK1133
C                                                                       MTPK1134
  200 DO 430 K = 1,N                                                    MTPK1131
      A(K,1) = FLOAT(K)                                                 MTPK1135
      A(K,2) = A(K,1)
  203 BMAX = DABS(B(K,K))                                               MTPK1137
  210 DO 219 I = K,N                                                    MTPK1138
  211 DO 219 J = K,N                                                    MTPK1139
      IF (BMAX - DABS(B(I,J)) )  213,219,219
  213 BMAX = DABS ( B(I,J))                                             MTPK1141
  214 A(K,1) = FLOAT(I)                                                 MTPK1142
  215 A(K,2) = FLOAT(J)                                                 MTPK1143
  219 CONTINUE                                                          MTPK1145
  216 IF (BMAX  -   1.D-38)   801,301,301
C                                                                       MTPK1146
C     EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL.              MTPK1147
C                                                                       MTPK1148
  301 I = IFIX(A(K,1))                                                  MTPK1149
  302 J = IFIX(A(K,2))                                                  MTPK1150
  313 DO 314 M = 1,N                                                    MTPK1151
      BMULT = B(I,M)                                                    MTPK1152
      B(I,M) = B(K,M)                                                   MTPK1153
  314 B(K,M) = BMULT                                                    MTPK1155
  303 DO 304 M = 1,N                                                    MTPK1156
      BMULT = B(M,J)                                                    MTPK1157
      B(M,J) = B(M,K)                                                   MTPK1158
  304 B(M,K) = BMULT                                                    MTPK1159
C                                                                       MTPK1160
C     ACTUAL INVERSION STEP.                                            MTPK1161
C     BEGIN ROW ITERATION                                               MTPK1162
C                                                                       MTPK1163
  401 DO 425 I = 1,N                                                    MTPK1164
C                                      IGNORE ROW K.                    MTPK1165
  402 IF (I  -   K) 405,425,405
C                                      SET ROW MULTIPLIER.              MTPK1167
 405  BMULT = B(I,K) / B(K,K)
C                                      MODIFY ROW ELEMENTS.             MTPK1169
  410 DO 420 J = 1,N                                                    MTPK1170
C                                      IGNORE COLUMN K                  MTPK1171
  411 IF (J  -   K)  414,412,414
  412 B(I,J) = -BMULT                                                   MTPK1173
  413 GO TO 420                                                         MTPK1174
  414 B(I,J) = B(I,J) - B(K,J) * BMULT                                  MTPK1175
  420 CONTINUE                                                          MTPK1176
  425 CONTINUE                                                          MTPK1177
C                                                                       MTPK1178
C     DIVIDE PIVOT ROW BY PIVOT ELEMENT                                 MTPK1179
C                                                                       MTPK1180
  426 BMULT = 1.D0/ B(K,K)
  427 B(K,K) = 1.D0
C                                      ACTUAL DIVISION                  MTPK1183
  428 DO 429 J = 1,N                                                    MTPK1184
  429 B(K,J) = B(K,J) * BMULT                                           MTPK1185
  430 CONTINUE                                                          MTPK1186
C                                                                       MTPK1187
C     REARRANGE AND REASSEMBLE INVERSE MATRIX.                          MTPK1188
C                                                                       MTPK1189
  501 DO 512 IJK = 1,N                                                  MTPK1190
  502 I = N-IJK&1                                                       MTPK1191
  503 L = IFIX(A(I,2))                                                  MTPK1192
  504 J = IFIX(A(I,1))                                                  MTPK1193
  505 DO 508 M = 1,N                                                    MTPK1194
  506 BMULT = B(I,M)                                                    MTPK1195
  507 B(I,M) = B(L,M)                                                   MTPK1196
  508 B(L,M) = BMULT                                                    MTPK1197
  509 DO 512 M = 1,N                                                    MTPK1198
  510 BMULT = B(M,I)                                                    MTPK1199
  511 B(M,I) = B(M,J)                                                   MTPK1200
  512 B(M,J) = BMULT                                                    MTPK1201
C                                                                       MTPK1202
C     MOVE A INVERSE BACK INTO A                                        MTPK1203
C                                                                       MTPK1204
  601 DO 603 I = 1,N                                                    MTPK1205
  602 DO 603 J = 1,N                                                    MTPK1206
  603 A(I,J) = B(I,J)
  701 M = 1                                                             MTPK1208
  702 RETURN                                                            MTPK1209
C                                                                       MTPK1210
C     ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38
C                                                                       MTPK1212
  801 M = 2                                                             MTPK1213
  802 RETURN                                                            MTPK1214
  999 END                                                               MTPK1215