C ccm.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. (1969), "Categorical conjoint measurement", unpublished
C manuscript, Bell Laboratories, Murray Hill, NJ.
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

$      FORTRAN
$      FORTRAN
C     CATEGORICAL CONJOINT MEASUREMENT
C     JULY, 1969, CHANG FOR CARROLL
      DIMENSION A(2000),V(2000),P(2000),ANCL(2000)
      DIMENSION TIT(12),NWW(4)
C     NWAY - NUMBER OF WAYS (OR FACTORS).  MAX.=4
C     NC - NUMBER OF CATEGORIES.  MAX.=20
C     NL - NUMBER OF SOLUTIONS.  MAX .LE. NC
C     NWW - NUMBER OF LEVELS WITHIN FACTORS.  THE PRODUCT OF LEVELS
C           MUST BE .LE. 2000.
C     FLGEOF IS SYSTEM ROUTINE TO TEST EOF, IF SO SET NEND=1.
      NEND=0
      CALL FLGEOF(5,NEND)
 50   READ (5,100) NWAY,NC,NL
      IF(NEND.EQ.1) GOTO999
      READ (5,101) TIT
      DO1I=1,4
 1    NWW(I)=1
      READ (5,100) (NWW(I),I=1,NWAY)
      N1=NWW(1)
      N2=NWW(2)
      N3=NWW(3)
      N4=NWW(4)
      CALL READA(A,N1,N2,N3,N4,NC,V)
C     FIND MAXIMUM OF NWW
      MNW=0
      ND=0
      DO2I=1,NWAY
      ND=ND&NWW(I)
      IF(NWW(I).LE.MNW) GOTO2
      MNW=NWW(I)
 2    CONTINUE
      PRINT102,TIT
      PRINT103,NWAY,NC
      PRINT104,(NWW(I),I=1,NWAY)
      CALL CCM(A,N1,N2,N3,N4,NWAY,NWW,NC,NL,V,MNW,ANCL,P,ND)
      GOTO50
 100  FORMAT(20I4)
 101  FORMAT(12A6)
 102  FORMAT(33H1CATEGORICAL CONJOINT MEASUREMENT/1H012A6)
 103  FORMAT(6H0NWAY=,I2/14H CATEGORIES - ,I3)
 104  FORMAT(10H LEVELS - ,4I3)
 999  STOP
      END
$      FORTRAN
      SUBROUTINE READA(A,N1,N2,N3,N4,NC,V)
      INTEGER A
      DIMENSION A(N1,N2,N3,N4),SYM(20),FMT(12),V(N1,N2,N3,N4)
      READ (5,100) FMT
 100  FORMAT(12A6)
      READ (5,FMT) (SYM(I),I=1,NC)
      DO5I4=1,N4
      DO5I3=1,N3
      DO5I2=1,N2
 5    READ (5,FMT) (V(I1,I2,I3,I4),I1=1,N1)
C     CHANGE FROM SYMBOLS TO NUMERALS
      DO6I4=1,N4
      DO6I3=1,N3
      DO6I2=1,N2
      DO6I1=1,N1
      DO7M=1,NC
      IF(V(I1,I2,I3,I4).EQ.SYM(M)) GOTO8
 7    CONTINUE
      PRINT101,I1,I2,I3,I4
 101  FORMAT(6H0CELL(,4I3,35H) SYMBOL DOES NOT MATCH, CHECK DATA)
      STOP
 8    A(I1,I2,I3,I4)=M
 6    CONTINUE
      RETURN
      END
$      FORTRAN     SUBROUTINE CCM, R=P*PT, DOES ANOVA
      SUBROUTINE CCM(A,N1,N2,N3,N4,NWAY,NWW,NC,NL,V,MNW,ANCL,P,ND)
      DIMENSION A(N1,N2,N3,N4),V(N1,N2,N3,N4),NWW(1),CCT(20)
      DIMENSION ANCL(MNW,4,NC),P(ND,NC),R(20,20),B(20,20)
      DIMENSION AM(50,4,10),XNK(4),XNWW(4)
      DIMENSION SS(4,10),SSE(10),XMS(4,10),EMS(10),F(4,10),NDF(4)
C     NCT(I) - NUMBER OF OCCURENCES OF CATEGORY I.
C     ANCL(I,J,K) - NUMBER OF OCCURENCES OF CAT. K FOR FACTOR J LEV. I.
      INTEGER A
      NT=1
      DO16I=1,NWAY
      NDF(I)=NWW(I)-1
      XNWW(I)=NWW(I)
 16   NT=NT*NWW(I)
      NDFE=NT-1-ND&NWAY
      EDF=NDFE
      XNT=NT
      DO17I=1,NWAY
 17   XNK(I)=NT/NWW(I)
      PRINT113
      DO30I4=1,N4
      DO30I3=1,N3
      DO30I1=1,N1
 30   PRINT114,I1,(A(I1,I2,I3,I4),I2=1,N2)
C     COMPUTE NUMBER OF OCCURENCES
      CALL ZEROT(ANCL,MNW,4,NC)
      DO3I=1,NC
 3    CCT(I)=0.
      DO5I4=1,N4
      DO5I3=1,N3
      DO5I2=1,N2
      DO5I1=1,N1
      M=A(I1,I2,I3,I4)
      CCT(M)=CCT(M)&1.
      GOTO(6,7,8,9),NWAY
 9    ANCL(I4,4,M)=ANCL(I4,4,M)&1.
 8    ANCL(I3,3,M)=ANCL(I3,3,M)&1.
 7    ANCL(I2,2,M)=ANCL(I2,2,M)&1.
 6    ANCL(I1,1,M)=ANCL(I1,1,M)&1.
 5    CONTINUE
      PRINT101,(I,I=1,NC)
      PRINT102,( CCT(I),I=1,NC)
      PRINT103
      DO1I=1,NWAY
      L=NWW(I)
      DO1J=1,L
 1    PRINT104,I,J,(ANCL(J,I,M),M=1,NC)
      DO10M=1,NC
      L=0
      DO10I=1,NWAY
      TEMP=SQRT(1./(XNK(I)*CCT(M)))
      TEM=CCT(M)/XNWW(I)
      K=NWW(I)
      DO10J=1,K
      L=L&1
 10   P(L,M)=TEMP*(ANCL(J,I,M)-TEM)
      L=0
      DO20I=1,NWAY
 20   XNWW(I)=XNWW(I)-1.
C     COMPUTE R
      DO11I=1,NC
      DO11J=1,NC
      R(I,J)=0.
      DO11M=1,ND
 11   R(I,J)=R(I,J)&P(M,I)*P(M,J)
      MM=1
      CALL MTEVV(R,20,20,B,20,20,NC,MM)
      PRINT108,(R(I,I),I=1,NC)
      PRINT109
      DO14I=1,NC
      CCT(I)=SQRT(CCT(I))
 14   PRINT107,I,(B(J,I),J=1,NC)
      DO42I=1,NC
      DO42J=1,NC
 42   B(J,I)=B(J,I)/CCT(J)
      PRINT122
      DO45I=1,NC
 45   PRINT107,I,(B(J,I),J=1,NC)
      CALL ZEROT(AM,50,4,10)
      DO13L=1,NL
      DO12I4=1,N4
      DO12I3=1,N3
      DO12I2=1,N2
      DO12I1=1,N1
      M=A(I1,I2,I3,I4)
 12   V(I1,I2,I3,I4)=B(M,L)
C     PUNCH REAL VALUED FACTORS
      PUNCH115,L
      PUNCH118
      DO15I4=1,N4
      DO15I3=1,N3
      DO15I1=1,N1
 15   PUNCH117,I1,(V(I1,I2,I3,I4),I2=1,N2)
C     COMPUTE MARGINAL MEANS
      DO21I4=1,N4
      DO21I3=1,N3
      DO21I2=1,N2
      DO21I1=1,N1
      Z=V(I1,I2,I3,I4)
      GOTO(22,23,24,25),NWAY
 25   AM(I4,4,L)=AM(I4,4,L)&Z
 24   AM(I3,3,L)=AM(I3,3,L)&Z
 23   AM(I2,2,L)=AM(I2,2,L)&Z
 22   AM(I1,1,L)=AM(I1,1,L)&Z
 21   CONTINUE
      SSQ=0.
      DO18I=1,NWAY
      M=NWW(I)
      SS(I,L)=0.
      DO36J=1,M
      AM(J,I,L)=AM(J,I,L)/XNK(I)
 36   SS(I,L)=SS(I,L)&AM(J,I,L)**2
      SS(I,L)=SS(I,L)*XNK(I)
      XMS(I,L)=SS(I,L)/XNWW(I)
 18   SSQ=SSQ&SS(I,L)
      SSE(L)=1.-SSQ
      EMS(L)=SSE(L)/EDF
      DO37I=1,NWAY
 37   F(I,L)=XMS(I,L)/EMS(L)
      DO32I4=1,N4
      DO32I3=1,N3
      DO32I2=1,N2
      DO32I1=1,N1
 32   V(I1,I2,I3,I4)=AM(I1,1,L)&AM(I2,2,L)&AM(I3,3,L)&AM(I4,4,L)
      PRINT125,L
      DO33I4=1,N4
      DO33I3=1,N3
      DO33I1=1,N1
 33   PRINT107,I1,(V(I1,I2,I3,I4),I2=1,N2)
 13   CONTINUE
      PRINT110
      DO19I=1,NWAY
      PRINT111,I
      PRINT112,(L,L=1,NL)
      M=NWW(I)
      DO19J=1,M
 19   PRINT116,J,(AM(J,I,L),L=1,NL)
      PRINT119
      DO38L=1,NL
      PRINT120,L
      PRINT121
      DO40I=1,NWAY
 40   PRINT123,I,SS(I,L),NDF(I),XMS(I,L),F(I,L)
 38   PRINT124,SSE(L),NDFE,EMS(L)
 101  FORMAT(38H0NUMBER OF OCCURENCES OF EACH CATEGORY/(10I6)/)
 102  FORMAT(10F6.1)
 103  FORMAT(43H0OCCURRENCES OF EACH CATEGORY ON EACH LEVEL/)
 104  FORMAT(2I2,3X,10E12.4/(7X,10E12.4))
 107  FORMAT(I3,3X,10E12.4/(6X,10E12.4))
 108  FORMAT(12H0EIGEN ROOTS/(2X,10E12.4))
 109  FORMAT(14H0EIGEN VECTORS)
 110  FORMAT(15H1MARGINAL MEANS)
 111  FORMAT(7H0FACTOR,I3)
 112  FORMAT(2X,5HLEVEL,5X,8HSOLUTION,I2,I10,6I15)
 113  FORMAT(5H0DATA)
 114  FORMAT(I4,2X,40I3/(6X,40I3))
 115  FORMAT(36H REAL VALUED FACTORS - FROM SOLUTION,I2)
 116  FORMAT(I5,1X,8E15.5)
 117  FORMAT(I3,2X,5E13.5/(5X,5E13.5))
 118  FORMAT(11H(5X,5E13.5))
 119  FORMAT(21H1ANALYSIS OF VARIANCE)
 120  FORMAT(9H0SOLUTIONI3)
 121  FORMAT(4H0WAY,8X,3HSSQ,10X,2HDF,11X,2HMS,15X,1HF)
 122  FORMAT(16H0CATEGORY VALUES)
 123  FORMAT(I3,3X,E15.5,I6,2E18.5)
 124  FORMAT(4H RES,2X,E15.5,I6,3X,E15.5)
 125  FORMAT(26H0ADDITIVE TABLE - SOLUTION,I2)
      RETURN
      END

$      FORTRAN
      SUBROUTINE ZEROS(A,N)
      DIMENSION A(1)
      K=N
 10   DO5II=1,K
 5    A(II)=0.
      GOTO50
      ENTRY ZEROD(A,I,J)
      K=I*J
      GOTO10
      ENTRY ZEROT(A,I,J,M)
      K=I*J*M
      GOTO10
 50   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
$      FORTRAN
$      INCODE  IBMF
*EIGENJ      EIGENVALUE AND EIGENVECTOR SUBROUTINE                      EIGJ0020
*                              CD600D4.009      DATE 05/05/65           EIGJ0030
      SUBROUTINE EIGENJ (A,V,NN,IV,E,I3,I4)
      DIMENSION A(5050),V( I3, I4),E(NN)
      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
      GO TO 27                                                          EIGJ0320
*     TRANSFORMATION SETUP FOLLOWS                                      EIGJ0330
    8 II=K&I-J                                                          EIGJ0340
      JJ=(J*(N&N-J&3))/2-N                                              EIGJ0350
      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
$      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)
      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
      CALL SORT(E(1),E(2),N,MS(1))
       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