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

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 "Individual Differences and Multidimensional Scaling" by J Douglas
C Carroll in "Multidimensional Scaling: Theory and Applications in the
C Behavioral Sciences", vol. 1, Theory, Seminar Press, New York and
C London, 1972
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

C     MDPREF - MULTIDIMENSIONAL ANALYSIS OF PREFERENCE DATA             MDPF  10
C     BY J D CARROLL AND J-J CHANG                                      MDPF  20
C     MAY 17,1968                                                       MDPF  30
      REAL*8 SCR(3600)                                                  MDPF  40
      DIMENSION W(60,60)                                                MDPF  50
      DIMENSION TIT(18) , BLK(6)                                        MDPF  60
      DIMENSION D(60,60)                                                MDPF  70
      DIMENSION WT(60,60)                                               MDPF  80
      DIMENSION VI(60,60)                                               MDPF  90
      REAL*4 A(120),B(120)
      INTEGER*4 SP
      DIMENSION X(60,60)                                                MDPF 100
       DIMENSION SPT(3600)                                              MDPF 110
       EQUIVALENCE  (VI(1,1),SPT(1))                                    MDPF 120
      COMMON W                                                          MDPF 130
      DATA BBK/'    '/                                                  MDPF 140
      DO 5 I=1,6                                                        MDPF 150
    5 BLK(I) = BBK                                                      MDPF 151
   60 READ303,(TIT(I),I=1,18)                                           MDPF 160
      DO 10 I=1,6                                                       MDPF 170
      IF(TIT(I)-BLK(I))11,10,11                                         MDPF 180
 10   CONTINUE                                                          MDPF 190
      GO TO 50                                                          MDPF 200
 11   READ 300,NP,NS,NF,NFP,IREAD,MDATA,NS1,NPUNCH,IPUNF                MDPF 210
C     PRINT TITLE AND OPTIONS                                           MDPF 220
 26   PRINT 301,(TIT(I),I=1,18)                                         MDPF 230
      PRINT 302,NP,NS,NF,NFP,IREAD,MDATA,NS1,NPUNCH  ,IPUNF             MDPF 240
      IF(IREAD.GT.0) GO TO 42                                           MDPF 250
C     COMPUTE FIRST SCORE MATRIX                                        MDPF 260
 41   CALL COMPW(NP,NS,D,W,WT,MDATA,IPUNF,NS1,SCR,VI)                   MDPF 270
      GO TO 51                                                          MDPF 280
C     READ IN FIRST SCORE MATRIX                                        MDPF 290
 42   CALL RFSM(W,NP,NS,IREAD,IPUNF)                                    MDPF 300
C     COMPUTE SECOND SCORE MATRIX                                       MDPF 310
 51   CALL COMPR(NP,NS,NF,W,D,X,WT,NPUNCH)                              MDPF 320
C                                                                       MDPF 330
C     HERE IS WHERE PLOTTING ROUTINE WILL GO                            MDPF 340
C     PEOPLE PLOT                                                       MDPF 350
      WRITE(6,400)                                                      MDPF 360
      PTWO=+1.5                                                         MDPF 380
      AMTWO=-1.5                                                        MDPF 390
      PONE = +1                                                         MDPF 400
       YMAX = +1.                                                       MDPF 410
       YMIN = -1.                                                       MDPF 420
      DO 410 I=1,NP
      A(I) = X(I,1)
  410 B(I) = X(I,2)
      CALL PLOTX(A,B,NP)
  400 FORMAT('1',10X,'PLOT OF SUBJECTS IN FIRST TWO DIMENSIONS ')
      WRITE(6,401)
  401 FORMAT('1',10X,'PLOT OF STIMULUS POINTS IN FIRST TWO DIMENSIONS')
      DO 403 I=1,NS
      A(I) = WT(I,1)
  403 B(I) = WT(I,2)
      CALL PLOTX(A,B,NS)
      WRITE(6,402)
  402 FORMAT('1',10X,'PLOT OF STIMULI AND SUBJECTS IN FIRST TWO DIMENSIO
     XNS ')
      SP = NP+NS
      DO 417 I=1,NS
      A(I) = WT(I,1)
  417 B(I) = WT(I,2)
      DO 404 I=1,NP
      J=I+NS
      A(J) = X(I,1)
  404 B(J) = X(I,2)
      CALL PLOTX(A,B,SP)
      WRITE(6,405)
      WRITE(6,405) NS,NP
  405 FORMAT(//3X,'FIRST',1X,I3,1X,'POINTS ARE STIMULI AND NEXT',1X,I3,1
     XX,'POINTS ARE SUBJECTS '/)
      GO TO 60                                                          MDPF 590
 300  FORMAT(18I4)                                                      MDPF 600
  301 FORMAT(1H1,18A4)                                                  MDPF 610
 302  FORMAT(66H0 NP  NS  NF  NFP  IREAD  MDATA  NS1  NPUNCH  IPUNF  MIDMDPF 620
     1EN  LP  KV/4I4,I6,I7,I6,3I7,I6,I4)                                MDPF 630
  303 FORMAT(18A4)                                                      MDPF 640
 50   CALL EXIT                                                         MDPF 650
      STOP                                                              MDPF 660
      END                                                               MDPF 670
      SUBROUTINE COMPW(NP,NS,D,W,WT,MDATA,IPUNF,NS1,SCR,VI)             MDPF 680
      REAL*8 SCR(3600)                                                  MDPF 690
      DIMENSION D(60,60)                                                MDPF 700
      DIMENSION WT(60,60)                                               MDPF 710
      DIMENSION W(60,60)                                                MDPF 720
      DIMENSION FMT(18),CODE(4)                                         MDPF 730
      DIMENSION Y(3600)                                                 MDPF 740
      DIMENSION VI(60,60)                                               MDPF 750
      READ 100,(FMT(I),I=1,18)                                          MDPF 760
      XM=NS                                                             MDPF 770
      READFMT,(CODE(I),I=1,4)                                           MDPF 780
      DO 18 L=1,NP                                                      MDPF 790
C     READ IN MATRIX D AND CONVERT                                      MDPF 800
      IF(MDATA.GT.0) GO TO 11                                           MDPF 810
      CALL READS(D,FMT,CODE,NS)                                         MDPF 820
      GO TO 406                                                         MDPF 830
C     READ IN DATA AND WEIGHT MATRIX                                    MDPF 840
 11   GO TO (12,13,12),MDATA                                            MDPF 850
 13   CALL READW(D,WT,FMT,CODE,NS)                                      MDPF 860
      GO TO 17                                                          MDPF 870
C     READ IN DATA AND GENERATE WEIGHT MATRIX FOR MISSING DATA          MDPF 880
 12   CALL READMW(D,WT,FMT,CODE,NS)                                     MDPF 890
      IF(MDATA.EQ.3) GO TO 19                                           MDPF 900
 17   CALL LSQMD(D,WT,Y,NS,NP,W,L,NS1,SCR,VI)                           MDPF 910
      GO TO 18                                                          MDPF 920
C     MISSING DATA SPECIAL BLOCK CONSTRUCTION                           MDPF 930
 19   CALL MDN(Y,WT,NS)                                                 MDPF 940
C     COMPUTE FIRST SCORE MATRIX                                        MDPF 950
 406  SWS=0.                                                            MDPF 960
      DO 10 I=1,NS                                                      MDPF 970
      SR=0.                                                             MDPF 980
      SC=0.                                                             MDPF 990
 9    DO 7 J=1,NS                                                       MDPF1000
      SR=SR+D(I,J)                                                      MDPF1010
 7    SC=SC+D(J,I)                                                      MDPF1020
 5    IF(MDATA.EQ.3) GO TO 8                                            MDPF1030
      W(L,I)=SR-SC                                                      MDPF1040
      GO TO 10                                                          MDPF1050
 8    W(L,I)=(SR-SC)/Y(I)                                               MDPF1060
 10   SWS=SWS+W(L,I)                                                    MDPF1070
      WMEAN=SWS/XM                                                      MDPF1080
      SUM=0.                                                            MDPF1090
      SUMN=0.                                                           MDPF1100
      DO 15 I=1,NS                                                      MDPF1110
      W(L,I)=W(L,I)-WMEAN                                               MDPF1120
      IF(MDATA.EQ.0) GO TO 15                                           MDPF1130
      WSQ=W(L,I)**2                                                     MDPF1140
      SUMN=SUMN+WSQ*Y(I)                                                MDPF1150
      SUM=SUM+WSQ                                                       MDPF1160
 15   CONTINUE                                                          MDPF1170
      IF(MDATA.EQ.0) GO TO 18                                           MDPF1180
      C=SQRT(SUMN/SUM)                                                  MDPF1190
      DO 16 I=1,NS                                                      MDPF1200
 16   W(L,I)=C*W(L,I)                                                   MDPF1210
 18   CONTINUE                                                          MDPF1220
C     PRINT W                                                           MDPF1230
      PRINT 101                                                         MDPF1240
      DO 20 I=1,NP                                                      MDPF1250
 20   PRINT 102,I,(W(I,J),J=1,NS)                                       MDPF1260
      IF(IPUNF)22,25,22                                                 MDPF1270
C     PUNCH FIRST SCORE MATRIX                                          MDPF1280
 22   PUNCH 106                                                         MDPF1290
      DO 21 I=1,NP                                                      MDPF1300
 21   PUNCH 107,I,(W(I,J),J=1,NS)                                       MDPF1310
  100 FORMAT(18A4)                                                      MDPF1320
 101  FORMAT(19H0FIRST SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS)            MDPF1330
 102  FORMAT(1H I4,10F12.4/(5X,10F12.4))                                MDPF1340
 103  FORMAT(3H C=3F13.4)                                               MDPF1350
 106  FORMAT(18HFIRST SCORE MATRIX/11H(2X,5F13.5))                      MDPF1360
 107  FORMAT(I2,5F13.5/(2X,5F13.5))                                     MDPF1370
 25   RETURN                                                            MDPF1380
      END                                                               MDPF1390
      SUBROUTINE COMPR(NP,NS,NF,W,R,X,Y,NPUNCH)                         MDPF1400
      DIMENSION Y(60,60)                                                MDPF1410
      DIMENSION R(60,60)                                                MDPF1420
      DIMENSION W(60,60)                                                MDPF1430
      DIMENSION X(60,60)                                                MDPF1440
C     COMPUTE CROSS PRODUCT MATRIX OF SUBJECTS Y=W(WT)                  MDPF1450
      DO 41 I=1,NP                                                      MDPF1460
      DO 41 J=I,NP                                                      MDPF1470
      Y(I,J)=0.                                                         MDPF1480
      DO 43 L=1,NS                                                      MDPF1490
 43   Y(I,J)=Y(I,J)+W(I,L)*W(J,L)                                       MDPF1500
 41   Y(J,I)=Y(I,J)                                                     MDPF1510
      PRINT 109                                                         MDPF1520
      DO 42 I=1,NP                                                      MDPF1530
 42   PRINT 104,I,(Y(I,J),J=1,NP)                                       MDPF1540
C     COMPUTE CORRELATION MATRIX OF SUBJECTS AND PRINT OUT              MDPF1550
      DO 44 I=2,NP                                                      MDPF1560
      II=I-1                                                            MDPF1570
      DO 44 J=1,II                                                      MDPF1580
      R(I,J)=Y(I,J)/SQRT(Y(I,I)*Y(J,J))                                 MDPF1590
      R(J,I)=R(I,J)                                                     MDPF1600
 44   CONTINUE                                                          MDPF1610
      DO 45 I=1,NP                                                      MDPF1620
 45   R(I,I)=1.                                                         MDPF1630
      PRINT 111                                                         MDPF1640
      DO 46 I=1,NP                                                      MDPF1650
 46   PRINT 104,I,(R(I,J),J=1,NP)                                       MDPF1660
C     COMPUTE CROSS PRODUCT MATRIX OF STIMULI R=WT(W)                   MDPF1670
      DO 20 I=1,NS                                                      MDPF1680
      DO 20 J=I,NS                                                      MDPF1690
      R(I,J)=0.                                                         MDPF1700
      DO 21 K=1,NP                                                      MDPF1710
 21   R(I,J)=R(I,J)+W(K,I)*W(K,J)                                       MDPF1720
 20   R(J,I)=R(I,J)                                                     MDPF1730
      PRINT 110                                                         MDPF1740
      DO 40 I=1,NS                                                      MDPF1750
 40   PRINT 104,I,(R(I,J),J=1,NS)                                       MDPF1760
C     FIND ROOTS AND VECTORS                                            MDPF1770
      M=1                                                               MDPF1780
      PRINT 311                                                         MDPF1790
      IF(NS.LE.NP) GO TO 61                                             MDPF1800
C     THE PURPOSE OF FINDING ROOTS AND VECTORS OF THE LESSER ORDER IS TOMDPF1810
C     SAVE COMPUTOR TIME                                                MDPF1820
C     FIND ROOTS AND VECTORS OF Y, IF NP.LT.NS.                         MDPF1830
C     CALL MTEVV(Y,100,100,R,100,100,NP,M)                              MDPF1840
      CALL MTEVV(Y,60,60,R,60,60,NP,M)                                  MDPF1850
C     CALL MTEVV(Y,20,20,R,20,20,NP,M)                                  MDPF1860
      PRINT 312,(Y(I,I),I=1,NP)                                         MDPF1870
      DO 47 I=1,NF                                                      MDPF1880
 47   Y(I,I)=SQRT(Y(I,I))                                               MDPF1890
C     COMPUTE POPULATION MATRIX X                                       MDPF1900
      DO 62 J=1,NP                                                      MDPF1910
      SUMQ=0.                                                           MDPF1920
      DO 67 I=1,NF                                                      MDPF1930
      X(J,I)=R(J,I)*Y(I,I)                                              MDPF1940
      SUMQ=SUMQ+X(J,I)**2                                               MDPF1950
 67   R(J,I)=R(J,I)/Y(I,I)                                              MDPF1960
      SUMQ=SQRT(SUMQ)                                                   MDPF1970
      DO 66 I=1,NF                                                      MDPF1980
 66   X(J,I)=X(J,I)/SUMQ                                                MDPF1990
 62   CONTINUE                                                          MDPF2000
C     COMPUTE STIMULUS MATRIX Y                                         MDPF2010
      DO 63 J=1,NF                                                      MDPF2020
      SUM=0.                                                            MDPF2030
      DO 68 I=1,NS                                                      MDPF2040
      Y(I,J)=0.                                                         MDPF2050
      DO 64 M=1,NP                                                      MDPF2060
 64   Y(I,J)=Y(I,J)+W(M,I)*R(M,J)                                       MDPF2070
 68   SUM=SUM+Y(I,J)**2                                                 MDPF2080
      SUM=SQRT(SUM)                                                     MDPF2090
C     NORMALIZE Y                                                       MDPF2100
      DO 65 L=1,NS                                                      MDPF2110
 65   Y(L,J)=Y(L,J)/SUM                                                 MDPF2120
 63   CONTINUE                                                          MDPF2130
      GO TO 28                                                          MDPF2140
C     FIND ROOTS AND VECTORS OF R, IF NS.LE.NP.                         MDPF2150
C61   CALL MTEVV(R,100,100,Y,100,100,NS,M)                              MDPF2160
   61 CALL MTEVV(R,60,60,Y,60,60,NS,M)                                  MDPF2170
C     COMPUTE POPULATION MATRIX X                                       MDPF2180
      PRINT 312,(R(I,I),I=1,NS)                                         MDPF2190
      DO 30 I=1,NP                                                      MDPF2200
      SUMQ=0.                                                           MDPF2210
      DO 25 K=1,NF                                                      MDPF2220
      X(I,K)=0.                                                         MDPF2230
      DO 22 J=1,NS                                                      MDPF2240
 22   X(I,K)=X(I,K)+W(I,J)*Y(J,K)                                       MDPF2250
 25   SUMQ=SUMQ+X(I,K)**2                                               MDPF2260
      SUMQ=SQRT(SUMQ)                                                   MDPF2270
      DO 29 K=1,NF                                                      MDPF2280
 29   X(I,K)=X(I,K)/SUMQ                                                MDPF2290
 30   CONTINUE                                                          MDPF2300
C     COMPUTE SECOND SCORE MATRIX AND STORE IN W                        MDPF2310
 28   DO 329 I=1,NP                                                     MDPF2320
      DO 329 J=1,NS                                                     MDPF2330
      W(I,J)=0.                                                         MDPF2340
      DO 328 K=1,NF                                                     MDPF2350
 328  W(I,J)=W(I,J)+X(I,K)*Y(J,K)                                       MDPF2360
 329  CONTINUE                                                          MDPF2370
      PRINT 300                                                         MDPF2380
      DO 35 I=1,NP                                                      MDPF2390
 35   PRINT 104,I,(W(I,J),J=1,NS)                                       MDPF2400
C     PRINT AND PUNCH X AND Y                                           MDPF2410
      PRINT 103                                                         MDPF2420
      DO 53 I=1,NF                                                      MDPF2430
 53   PRINT 104,I,(X(J,I),J=1,NP)                                       MDPF2440
      PRINT 105                                                         MDPF2450
      DO 54 I=1,NF                                                      MDPF2460
 54   PRINT 104,I,(Y(J,I),J=1,NS)                                       MDPF2470
      IF(NPUNCH)51,50,51                                                MDPF2480
 51   PUNCH 106                                                         MDPF2490
      DO 57 I=1,NF                                                      MDPF2500
 57   PUNCH 107,I,(Y(J,I),J=1,NS)                                       MDPF2510
      PUNCH 108                                                         MDPF2520
      DO 58 I=1,NF                                                      MDPF2530
 58   PUNCH 107,I,(X(J,I),J=1,NP)                                       MDPF2540
 103  FORMAT(18H0POPULATION MATRIX/7H0FACTOR)                           MDPF2550
 104  FORMAT(1H I4,10F12.4/(5X,10F12.4))                                MDPF2560
 105  FORMAT(16H0STIMULUS MATRIX/7H0FACTOR)                             MDPF2570
 106  FORMAT(15HSTIMULUS MATRIX/11H(2X,5F13.5))                         MDPF2580
 107  FORMAT(I2,5F13.5/(2X,5F13.5))                                     MDPF2590
 108  FORMAT(17HPOPULATION MATRIX/11H(2X,5F13.5))                       MDPF2600
 109  FORMAT(33H0CROSS PRODUCT MATRIX OF SUBJECTS)                      MDPF2610
 110  FORMAT(32H0CROSS PRODUCT MATRIX OF STIMULI)                       MDPF2620
 111  FORMAT(31H0CORRELATION MATRIX OF SUBJECTS)                        MDPF2630
 300  FORMAT(20H0SECOND SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS)           MDPF2640
 311  FORMAT(32H0ROOTS OF THE FIRST SCORE MATRIX)                       MDPF2650
 312  FORMAT(3X,10F12.4)                                                MDPF2660
 50   RETURN                                                            MDPF2670
      END                                                               MDPF2680
      SUBROUTINE LSQMD(D,WT,Y,NS,NP,W,L,NS1,SCR,VI)                     MDPF2690
      REAL*8 SCR(3600)                                                  MDPF2700
      DIMENSION Y(3600)                                                 MDPF2710
      DIMENSION D(60,60)                                                MDPF2720
      DIMENSION WT(60,60)                                               MDPF2730
      DIMENSION VI(60,60)                                               MDPF2740
      DIMENSION W(60,60)                                                MDPF2750
C     BEFORE MAKING ANY CHANGES IN THIS ROUTINE, CHECK THE EQUIVALENCE  MDPF2760
C     STATEMENT IN MAIN PROGRAM INVOLVING SCR,D AND VI.                 MDPF2770
C     DATA IN D AND WEIGHT IN WT                                        MDPF2780
      PRINT 100,L                                                       MDPF2790
      BIG=10.**8                                                        MDPF2800
C     COMPUTE TOTAL VARIANCE                                            MDPF2810
      TVAR=0.                                                           MDPF2820
      DO 2 I=1,NS                                                       MDPF2830
      Y(I)=0.                                                           MDPF2840
      DO 3 J=1,NS                                                       MDPF2850
      TVAR=TVAR+WT(I,J)*D(I,J)**2                                       MDPF2860
 3    Y(I)=Y(I)+WT(I,J)*D(I,J)-WT(J,I)*D(J,I)                           MDPF2870
 2    CONTINUE                                                          MDPF2880
      IF(L.EQ.1.OR.L.GT.NS1) GO TO 50                                   MDPF2890
      GO TO 52                                                          MDPF2900
C     COMPUTE U AND PUT IN D                                            MDPF2910
 50   DO 5 I=1,NS                                                       MDPF2920
      D(I,I)=0.                                                         MDPF2930
      DO 5 J=1,NS                                                       MDPF2940
      IF(I.EQ.J) GO TO 5                                                MDPF2950
      D(I,J)=-(WT(I,J)+WT(J,I))                                         MDPF2960
      D(I,I)=D(I,I)-D(I,J)                                              MDPF2970
 5    CONTINUE                                                          MDPF2980
      DO 8 I=1,NS                                                       MDPF2990
      DO 8 J=1,I                                                        MDPF3000
      WT(I,J)=ABS(D(I,J))                                               MDPF3010
      IF(D(I,J).EQ.0.) WT(I,J)=BIG                                      MDPF3020
 8    WT(J,I)=WT(I,J)                                                   MDPF3030
C     PARTITION STIMULI IN BLOCKS                                       MDPF3040
      CALL BLOCK(WT,NS)                                                 MDPF3050
 20   DO 21 I=1,NS                                                      MDPF3060
      DO 21 J=1,NS                                                      MDPF3070
      IF(I.EQ.J) GO TO 22                                               MDPF3080
      IF(WT(I,J).GE.BIG) GO TO 24                                       MDPF3090
 22   WT(I,J)=1.                                                        MDPF3100
      GO TO 21                                                          MDPF3110
 24   WT(I,J)=0.                                                        MDPF3120
 21   CONTINUE                                                          MDPF3130
C     COMPUTE V                                                         MDPF3140
      DO 25 I=1,NS                                                      MDPF3150
      DO 25 J=1,NS                                                      MDPF3160
 25   WT(I,J)=WT(I,J)+D(I,J)                                            MDPF3170
      MM=1                                                              MDPF3180
C     ON RETURN THE INVERSE IS STORED IN WT (ORDER NS).  SCR IS USED AS MDPF3190
C     SCRATCH                                                           MDPF3200
      CALL MTINV(WT,60,60,SCR,60,60,NS,NN)                              MDPF3210
C     CALL MTINV(WT,20,20,SCR,20,20,NS,NM)                              MDPF3220
C     CALL MTINV(WT,100,100,SCR,100,100,NS,MM)                          MDPF3230
      IF(MM.EQ.1) GO TO 53                                              MDPF3240
C     ON RETURN IF MM IS NOT 1 INDICATES WT IS ILL CONDITIONED.         MDPF3250
C     PROGRAM PRINTS OUT MESSAGE AND STOPS.                             MDPF3260
      PRINT 111,MM                                                      MDPF3270
      STOP                                                              MDPF3280
 53   DO 51 I=1,NS                                                      MDPF3290
      DO 51 J=1,NS                                                      MDPF3300
 51   VI(I,J)=WT(I,J)                                                   MDPF3310
 52   DO 26 I=1,NS                                                      MDPF3320
      W(L,I)=0.                                                         MDPF3330
      DO 26 J=1,NS                                                      MDPF3340
 26   W(L,I)=W(L,I)+Y(J)*VI(J,I)                                        MDPF3350
      PRINT 109                                                         MDPF3360
      PRINT 102,I,(Y(I),I=1,NS)                                         MDPF3370
C     COMPUTE VARIANCE ACCOUNTED FOR                                    MDPF3380
      VAF=0.                                                            MDPF3390
      DEN=0.                                                            MDPF3400
      DO 27 I=1,NS                                                      MDPF3410
      VAF=VAF+W(L,I)*Y(I)                                               MDPF3420
 27   DEN=DEN+W(L,I)**2                                                 MDPF3430
      PVAF=VAF/TVAR                                                     MDPF3440
C     NORMALIZE                                                         MDPF3450
      C=SQRT(VAF/DEN)                                                   MDPF3460
      DO 28 I=1,NS                                                      MDPF3470
 28   W(L,I)=W(L,I)*C                                                   MDPF3480
      PRINT 107,TVAR,VAF,PVAF,C                                         MDPF3490
 100  FORMAT(8H0SUBJECTI4)                                              MDPF3500
 102  FORMAT(I4,10F12.4/(4X,10F12.4))                                   MDPF3510
 107  FORMAT(6H0TVAR=F12.4,3X,4HVAF=F12.4,3X,5HPVAF=F12.4,3X,2HC=F12.4) MDPF3520
 109  FORMAT(9H0Y VECTOR)                                               MDPF3530
 111  FORMAT(4H0MM=I1,21HWT CANNOT BE INVERTED)                         MDPF3540
      RETURN                                                            MDPF3550
      END                                                               MDPF3560
      SUBROUTINE RFSM(W,NP,NS,IREAD,IPUNF)                              MDPF3570
      DIMENSION W(60,60)                                                MDPF3580
      DIMENSION WMEAN(100),SD(100),FMT(18)                              MDPF3590
      XM=NS                                                             MDPF3600
      DEN=1./XM                                                         MDPF3610
      READ 100,(FMT(I),I=1,18)                                          MDPF3620
      DO 5 I=1,NP                                                       MDPF3630
 5    READFMT,(W(I,J),J=1,NS)                                           MDPF3640
C     COMPUTE MEAN FOR EACH SUBJECT                                     MDPF3650
      DO 6 I=1,NP                                                       MDPF3660
      SUM=0.                                                            MDPF3670
      DO 4 J=1,NS                                                       MDPF3680
 4    SUM=SUM+W(I,J)                                                    MDPF3690
      WMEAN(I)=SUM/XM                                                   MDPF3700
C     NORMALIZE W - SUBTRACT THE ROW MEAN                               MDPF3710
      SUM=0.                                                            MDPF3720
      DO 7 J=1,NS                                                       MDPF3730
      W(I,J)=W(I,J)-WMEAN(I)                                            MDPF3740
      IF(IREAD.EQ.2)SUM=SUM+W(I,J)**2                                   MDPF3750
 7    CONTINUE                                                          MDPF3760
      IF(IREAD.EQ.1) GO TO 6                                            MDPF3770
C     IREAD=2, NORMALIZED W - DIVIDE BY THE ROW SD                      MDPF3780
      SD(I)=SQRT(SUM*DEN)                                               MDPF3790
      DO 8 J=1,NS                                                       MDPF3800
 8    W(I,J)=W(I,J)/SD(I)                                               MDPF3810
 6    CONTINUE                                                          MDPF3820
C     PRINT MEAN                                                        MDPF3830
      PRINT 103                                                         MDPF3840
      PRINT 104,(WMEAN(I),I=1,NP)                                       MDPF3850
      IF(IREAD.EQ.1) GO TO 10                                           MDPF3860
      PRINT 105                                                         MDPF3870
      PRINT 104,(SD(I),I=1,NP)                                          MDPF3880
C     PRINT W                                                           MDPF3890
 10   PRINT 101                                                         MDPF3900
      DO 20 I=1,NP                                                      MDPF3910
 20   PRINT 102,I,(W(I,J),J=1,NS)                                       MDPF3920
      IF(IPUNF)22,25,22                                                 MDPF3930
C     PUNCH FIRST SCORE MATRIX                                          MDPF3940
   22 PUNCH 106                                                         MDPF3950
      DO 21 I=1,NP                                                      MDPF3960
 21   PUNCH 107,I,(W(I,J),J=1,NS)                                       MDPF3970
  100 FORMAT(18A4)                                                      MDPF3980
 101  FORMAT(19H0FIRST SCORE MATRIX/8H0SUBJECT5X,8HSTIMULUS)            MDPF3990
 102  FORMAT(1H I4,10F12.4/(5X,10F12.4))                                MDPF4000
 103  FORMAT(37H0MEAN OF THE RAW SCORES  (BY SUBJECT))                  MDPF4010
 104  FORMAT(5X,10F12.4)                                                MDPF4020
 105  FORMAT(35H0SD OF THE RAW SCORES  (BY SUBJECT))                    MDPF4030
 106  FORMAT(18HFIRST SCORE MATRIX/11H(2X,5F13.5))                      MDPF4040
 107  FORMAT(I2,5F13.5/(2X,5F13.5))                                     MDPF4050
 25   RETURN                                                            MDPF4060
      END                                                               MDPF4070
      SUBROUTINE MDN(Y,WT,NS)                                           MDPF4080
      DIMENSION Y(3600)                                                 MDPF4090
      DIMENSION WT(60,60)                                               MDPF4100
      DIMENSION W(60,60)                                                MDPF4110
      COMMON W                                                          MDPF4120
      DO 6 I=1,NS                                                       MDPF4130
      SUM=0.                                                            MDPF4140
      SSQ=0.                                                            MDPF4150
      SDSQ=0.                                                           MDPF4160
      DO 5 J=1,NS                                                       MDPF4170
      IF(I.EQ.J) GO TO 5                                                MDPF4180
      SUM=SUM+WT(I,J)+WT(J,I)                                           MDPF4190
      SDSQ=SDSQ+(WT(J,I)-WT(I,J))**2                                    MDPF4200
      SSQ=SSQ+WT(J,I)**2+WT(I,J)**2                                     MDPF4210
 5    CONTINUE                                                          MDPF4220
      IF(SSQ)4,7,4                                                      MDPF4230
 7    Y(I)=2.                                                           MDPF4240
      GO TO 6                                                           MDPF4250
 4    Y(I)=2.+SUM-SQRT(SDSQ/SSQ)                                        MDPF4260
 6    CONTINUE                                                          MDPF4270
      PRINT 100                                                         MDPF4280
      PRINT 101,(Y(I),I=1,NS)                                           MDPF4290
 100  FORMAT(9H0VECTOR N)                                               MDPF4300
 101  FORMAT(2X,10F12.4)                                                MDPF4310
      RETURN                                                            MDPF4320
      END                                                               MDPF4330
      SUBROUTINE READS(D,FMT,CODE,NS,WT)                                MDPF4340
C     READ IN DATA MATRIX WITHOUT MISSING DATA                          MDPF4350
      DIMENSION D(60,60),WT(60,60)                                      MDPF4360
       DIMENSION  FMT(1),      CODE(1)                                  MDPF4370
      N=1                                                               MDPF4380
C     READ IN DATA MATRIX                                               MDPF4390
 70   DO 15 I=1,NS                                                      MDPF4400
 15   READFMT,(D(I,J),J=1,NS)                                           MDPF4410
      GO TO(1,2,3),N                                                    MDPF4420
C     NO MISSING DATA                                                   MDPF4430
 1    DO 50 I=1,NS                                                      MDPF4440
      DO 50 J=1,NS                                                      MDPF4450
      IF(I.EQ.J) GO TO 45                                               MDPF4460
      DO 32 K=1,4                                                       MDPF4470
      IF(D(I,J)-CODE(K))32,41,32                                        MDPF4480
 41   GO TO(10,20,45,45),K                                              MDPF4490
 10   D(I,J)=1.                                                         MDPF4500
      GO TO 50                                                          MDPF4510
 20   D(I,J)=-1.                                                        MDPF4520
      GO TO 50                                                          MDPF4530
 32   CONTINUE                                                          MDPF4540
 45   D(I,J)=0.                                                         MDPF4550
 50   CONTINUE                                                          MDPF4560
      RETURN                                                            MDPF4570
      ENTRY READMW(D,WT,FMT,CODE,NS)                                    MDPF4580
C     MISSING DATA, WEIGHT IS EITHER 0 OR 1                             MDPF4590
      N=2                                                               MDPF4600
      GO TO 70                                                          MDPF4610
 2    DO 55 I=1,NS                                                      MDPF4620
      DO 55 J=1,NS                                                      MDPF4630
      IF(I.EQ.J) GO TO 55                                               MDPF4640
      IF(D(I,J).EQ.CODE(4)) GO TO 53                                    MDPF4650
      WT(I,J)=1.                                                        MDPF4660
      GO TO 55                                                          MDPF4670
 53   WT(I,J)=0.                                                        MDPF4680
 55   CONTINUE                                                          MDPF4690
      GO TO 1                                                           MDPF4700
      ENTRY READW(D,WT,FMT,CODE,NS)                                     MDPF4710
C     READ IN WEIGHT MATRIX                                             MDPF4720
      N=3                                                               MDPF4730
      GO TO 70                                                          MDPF4740
 3    DO 62 I=1,NS                                                      MDPF4750
 62   READ 100,(WT(I,J),J=1,NS)                                         MDPF4760
 100  FORMAT(10F7.2)                                                    MDPF4770
      DO 65 I=1,NS                                                      MDPF4780
      DO 65 J=1,NS                                                      MDPF4790
      IF(I.EQ.J) GO TO 65                                               MDPF4800
      IF(D(I,J).EQ.CODE(4)) WT(I,J)=0.                                  MDPF4810
 65   CONTINUE                                                          MDPF4820
      GO TO 1                                                           MDPF4830
      END                                                               MDPF4840
      SUBROUTINE BLOCK(D,N)                                             MDPF4850
      DIMENSION D(60,60)                                                MDPF4860
C     COMPUTE MINIMUM PATH                                              MDPF4870
      DO 5 I=1,N                                                        MDPF4880
      DO 5 J=2,N                                                        MDPF4890
      IF(I.EQ.J) GO TO 5                                                MDPF4900
      JJ=J-1                                                            MDPF4910
      DO 4 K=1,JJ                                                       MDPF4920
      IF(I.EQ.J) GO TO 4                                                MDPF4930
      DSUM=D(I,J)+D(I,K)                                                MDPF4940
      IF(D(J,K).GT.DSUM) D(J,K)=DSUM                                    MDPF4950
      D(K,J)=D(J,K)                                                     MDPF4960
 4    CONTINUE                                                          MDPF4970
 5    CONTINUE                                                          MDPF4980
      RETURN                                                            MDPF4990
      END                                                               MDPF5000
       SUBROUTINE SORT (A, N, B, C, D, K, SWITCH )                      MDPF5010
CSORT                                                                   MDPF5020
C                                                                       MDPF5030
C     THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY,    MDPF5040
C     ARRAYS 'B', 'C', AND 'D', IN ORDER CORRESPONDING TO 'A'.          MDPF5050
C     N = NUMBER OF ITEMS IN 'A' (AND 'B', 'C', 'D', IF USED)           MDPF5060
C     K = 0--SORT 'A' ONLY, 1--REARRANGE 'B', 2--REARRANGE 'B' AND 'C', MDPF5070
C         3--REARRANGE 'B', 'C', AND 'D'.                               MDPF5080
C     IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER,         MDPF5090
C                 IF ZERO OR NEGATIVE, IN DESCENDING ORDER.             MDPF5100
C     ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL            MDPF5110
C                                                                       MDPF5120
      DIMENSION A(N),B(N),C(N),D(N)                                     MDPF5130
       INTEGER  SWITCH                                                  MDPF5140
 105   KP1=K+1                                                          MDPF5150
       IF(N.LE.1) GO TO 999                                             MDPF5160
      M = 1                                                             MDPF5170
  106  M = M + M                                                        MDPF5180
       IF( M .LT. N ) GO TO 106                                         MDPF5190
       M = M - 1                                                        MDPF5200
 994       M = M/2                                                      MDPF5210
       IF(M.EQ.0) GO TO 999                                             MDPF5220
       KK = N-M                                                         MDPF5230
       J = 1                                                            MDPF5240
992     I = J                                                           MDPF5250
996     IM = I + M                                                      MDPF5260
       IF(SWITCH)  810,810,800                                          MDPF5270
800     IF(A(I).GT.A(IM)) GO TO 110                                     MDPF5280
       GO TO 995                                                        MDPF5290
810    IF(A(I).LT.A(IM))  GO TO 110                                     MDPF5300
995     J = J+1                                                         MDPF5310
       IF(J.GT.KK)  GO TO 994                                           MDPF5320
       GO TO 992                                                        MDPF5330
 110   TEMP=A(I)                                                        MDPF5340
       A(I) = A(IM)                                                     MDPF5350
       A(IM) = TEMP                                                     MDPF5360
       GO TO ( 140, 130, 120, 115), KP1                                 MDPF5370
 115   TEMP = D(I)                                                      MDPF5380
       D(I) = D(IM)                                                     MDPF5390
        D(IM) = TEMP                                                    MDPF5400
 120   TEMP=C(I)                                                        MDPF5410
       C(I) = C(IM)                                                     MDPF5420
       C(IM) = TEMP                                                     MDPF5430
 130   TEMP=B(I)                                                        MDPF5440
       B(I) = B(IM)                                                     MDPF5450
       B(IM) = TEMP                                                     MDPF5460
140       I = I-M                                                       MDPF5470
       IF(I.LT.1)  GO TO 995                                            MDPF5480
        GO TO 996                                                       MDPF5490
999      RETURN                                                         MDPF5500
      END                                                               MDPF5510
      SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M)                            MDPF5520
C      INCODE  IBMF                                                     MDPF5530
C      FORTRAN DECK,LSTOU,STAB                                          MDPF5540
C                  MTINV                                                MDPF5550
       DOUBLE PRECISION  DABS                                           MDPF5560
      REAL*8 B(I3,I4),BMAX,BMULT                                        MDPF5570
      DIMENSION A(I1,I2)                                                MDPF5580
C                                                                       MDPF5590
C     THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS-  MDPF5600
C     JORDAN ELIMINATION SCHEME.  ALL CALCULATIONS ARE DONE IN DOUBLE-  MDPF5610
C     PRECISION ARITHMETIC                                              MDPF5620
C                                                                       MDPF5630
      EQUIVALENCE (BMAX,BMULT)                                          MDPF5640
C                                                                       MDPF5650
C     SINGLE TO DOUBLE PRECISION TRANSFER                               MDPF5660
C                                                                       MDPF5670
  101 DO 103 I = 1,N                                                    MDPF5680
  102 DO 103 J = 1,N                                                    MDPF5690
  103 B(I,J) = A(I,J)                                                   MDPF5700
C                                                                       MDPF5710
C     FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMARTIX. MDPF5720
C                                                                       MDPF5730
  200 DO 430 K = 1,N                                                    MDPF5740
      A(K,1) = FLOAT(K)                                                 MDPF5750
      A(K,2) = A(K,1)                                                   MDPF5760
  203 BMAX = DABS(B(K,K))                                               MDPF5770
  210 DO 219 I = K,N                                                    MDPF5780
  211 DO 219 J = K,N                                                    MDPF5790
      IF (BMAX - DABS(B(I,J)) )   213,219,219                           MDPF5800
  213 BMAX = DABS ( B(I,J))                                             MDPF5810
  214 A(K,1) = FLOAT(I)                                                 MDPF5820
  215 A(K,2) = FLOAT(J)                                                 MDPF5830
  219 CONTINUE                                                          MDPF5840
  216 IF (BMAX  -  1.D-38)   801,301,301                                MDPF5850
C                                                                       MDPF5860
C     EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL.              MDPF5870
C                                                                       MDPF5880
  301  I = IFIX(A(K,1))                                                 MDPF5890
  302  J = IFIX(A(K,2))                                                 MDPF5900
  313 DO 314 M = 1,N                                                    MDPF5910
      BMULT = B(I,M)                                                    MDPF5920
      B(I,M) = B(K,M)                                                   MDPF5930
  314 B(K,M) = BMULT                                                    MDPF5940
  303 DO 304 M = 1,N                                                    MDPF5950
      BMULT = B(M,J)                                                    MDPF5960
      B(M,J) = B(M,K)                                                   MDPF5970
  304 B(M,K) = BMULT                                                    MDPF5980
C                                                                       MDPF5990
C     ACTUAL INVERSION STEP.                                            MDPF6000
C     BEGIN ROW ITERATION                                               MDPF6010
C                              IGNORE ROW K                             MDPF6020
  401 DO 425 I = 1,N                                                    MDPF6030
C                                                                       MDPF6040
  402 IF (I  -  K) 405,425,405                                          MDPF6050
C                                       SET ROW MULTIPLIER.             MDPF6060
  405 BMULT = B(I,K) / B(K,K)                                           MDPF6070
C                                       MODIFY ROW ELEMENTS.            MDPF6080
  410 DO 420 J = 1,N                                                    MDPF6090
C                                       IGNORE COLUMN K                 MDPF6100
  411 IF (J  -  K)  414,412,414                                         MDPF6110
  412 B(I,J) = -BMULT                                                   MDPF6120
  413 GO TO 420                                                         MDPF6130
  414 B(I,J) = B(I,J) - B(K,J) * BMULT                                  MDPF6140
  420 CONTINUE                                                          MDPF6150
  425 CONTINUE                                                          MDPF6160
C                                                                       MDPF6170
C     DIVIDE PIVOT ROW BY PIVOT ELEMENT                                 MDPF6180
C                                                                       MDPF6190
  426 BMULT = 1.D0/B(K,K)                                               MDPF6200
  427 B(K,K) = 1.D0                                                     MDPF6210
C                                       ACTUAL DIVISION                 MDPF6220
  428 DO 429 J = 1,N                                                    MDPF6230
  429 B(K,J) = B(K,J) * BMULT                                           MDPF6240
  430 CONTINUE                                                          MDPF6250
C                                                                       MDPF6260
C     REARRANGE AND REASSEMBLE INVERSE MATRIX.                          MDPF6270
C                                                                       MDPF6280
  501 DO 512 IJK = 1,N                                                  MDPF6290
  502 I = N-IJK+1                                                       MDPF6300
  503 L = IFIX(A(I,2))                                                  MDPF6310
  504 J = IFIX(A(I,1))                                                  MDPF6320
  505 DO 508 M = 1,N                                                    MDPF6330
  506 BMULT = B(I,M)                                                    MDPF6340
  507 B(I,M) = B(L,M)                                                   MDPF6350
  508 B(L,M) = BMULT                                                    MDPF6360
  509 DO 512 M = 1,N                                                    MDPF6370
  510 BMULT = B(M,I)                                                    MDPF6380
  511 B(M,I) = B(M,J)                                                   MDPF6390
  512 B(M,J) = BMULT                                                    MDPF6400
C                                                                       MDPF6410
C     MOVE A INVERSE BACK INTO A                                        MDPF6420
C                                                                       MDPF6430
  601 DO 603 I = 1,N                                                    MDPF6440
  602 DO 603 J = 1,N                                                    MDPF6450
  603 A(I,J) = B(I,J)                                                   MDPF6460
  701 M = 1                                                             MDPF6470
  702 RETURN                                                            MDPF6480
C                                                                       MDPF6490
C     ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38                   MDPF6500
C                                                                       MDPF6510
  801 M = 2                                                             MDPF6520
  802 RETURN                                                            MDPF6530
      END                                                               MDPF6540
      SUBROUTINE ARREV(N,NS,D,I3,I4)                                    MDPF6550
      DIMENSION NS(N),D(I3,I4)                                          MDPF6560
      DO 10 I=1,N                                                       MDPF6570
    1 CONTINUE                                                          MDPF6580
      IF(NS(I)-I) 2,10,2                                                MDPF6590
    2 CONTINUE                                                          MDPF6600
      KK=I                                                              MDPF6610
      JJ=NS(I)                                                          MDPF6620
      CALL XCH(N,KK,JJ,D,NS,I3,I4)                                      MDPF6630
      GO TO 1                                                           MDPF6640
 10   CONTINUE                                                          MDPF6650
      RETURN                                                            MDPF6660
      END                                                               MDPF6670
      SUBROUTINE XCH(N,L,M,D,NS,I3,I4)                                  MDPF6680
      DIMENSION NS(N),D(I3,I4)                                          MDPF6690
      I=L                                                               MDPF6700
      J=M                                                               MDPF6710
      DO 1 K=1,N                                                        MDPF6720
      T=D(K,I)                                                          MDPF6730
      D(K,I)=D(K,J)                                                     MDPF6740
 1    D(K,J)=T                                                          MDPF6750
      NS(I)=NS(J)                                                       MDPF6760
      NS(J)=J                                                           MDPF6770
      RETURN                                                            MDPF6780
      END                                                               MDPF6790
       SUBROUTINE SORT1(A, N, B, C, D, K, SWITCH )                      MDPF6800
      DIMENSION A(N),B(N),C(N),D(N)                                     MDPF6810
      INTEGER SWITCH, SW, A, B, C, D, TEMP                              MDPF6820
      IF(SWITCH) 100, 105 , 100                                         MDPF6830
  100 SW = SWITCH                                                       MDPF6840
 105   KP1=K+1                                                          MDPF6850
      NM1 = N-1                                                         MDPF6860
      IF(NM1) 300,300,107                                               MDPF6870
  107 DO 150 I=1,NM1                                                    MDPF6880
      IP1 = I+1                                                         MDPF6890
      DO 140 J=IP1,N                                                    MDPF6900
      IF(A(I)-A(J)) 140 , 140 , 110                                     MDPF6910
  110 TEMP=A(I)                                                         MDPF6920
      A(I) = TEMP                                                       MDPF6930
      A(J) = TEMP                                                       MDPF6940
      GO TO(140,130, 120,115), KP1                                      MDPF6950
  115 TEMP=D(I)                                                         MDPF6960
      D(I) = D(J)                                                       MDPF6970
      D(J) = TEMP                                                       MDPF6980
  120 TEMP= C(I)                                                        MDPF6990
      C(I) = C(J)                                                       MDPF7000
      C(J) = TEMP                                                       MDPF7010
  130 TEMP= B(I)                                                        MDPF7020
      B(I) = B(J)                                                       MDPF7030
      B(J) = TEMP                                                       MDPF7040
  140 CONTINUE                                                          MDPF7050
  150 CONTINUE                                                          MDPF7060
      IF(SW ) 200,300,300                                               MDPF7070
  200 NHALF=N/2                                                         MDPF7080
      DO 230 I=1,NHALF                                                  MDPF7090
      IC=N+1-I                                                          MDPF7100
      TEMP= A(I)                                                        MDPF7110
      A(I) = A(IC)                                                      MDPF7120
      A(IC) = TEMP                                                      MDPF7130
      GO TO ( 230,220,210,205 ) ,KP1                                    MDPF7140
  205 TEMP= D(I)                                                        MDPF7150
      D(I) = D(J)                                                       MDPF7160
      D(IC) = TEMP                                                      MDPF7170
  210 TEMP= C(I)                                                        MDPF7180
      C(I) = C(IC)                                                      MDPF7190
      C(IC) = TEMP                                                      MDPF7200
  220 TEMP = B(I)                                                       MDPF7210
      B(I) = B(IC)                                                      MDPF7220
      B(IC) = TEMP                                                      MDPF7230
  230 CONTINUE                                                          MDPF7240
  300 CONTINUE                                                          MDPF7250
      RETURN                                                            MDPF7260
      END                                                               MDPF7270
      SUBROUTINE  MTEVV(X,I1,I2,D,I3,I4,N,M)                            MDPF7280
      DIMENSION A(1830),X(60,60),D(60,60),E(1830),S(220),F(60,60)       MDPF7290
      K=1                                                               MDPF7300
      DO 1 J=1,N                                                        MDPF7310
      DO 1 I=1,J                                                        MDPF7320
      A(K)=X(I,J)                                                       MDPF7330
 1    K=K+1                                                             MDPF7340
       CALL EIGEN(A,E,N,0)                                              MDPF7350
       II = 1                                                           MDPF7360
       DO 974 I=1,N                                                     MDPF7370
       DO 974 J=1,N                                                     MDPF7380
       D(J,I) = E(II)                                                   MDPF7390
       II = II + 1                                                      MDPF7400
  974  CONTINUE                                                         MDPF7410
       II = 0                                                           MDPF7420
       DO 973 J=1,N                                                     MDPF7430
       II = II + J                                                      MDPF7440
  973  X(J,J) = A(II)                                                   MDPF7450
      RETURN                                                            MDPF7460
      END                                                               MDPF7470
      SUBROUTINE EIGEN(A,R,N,MV)                                        MDPF7480
C                                                                       MDPF7490
C        SUBROUTINE EIGEN                                               MDPF7500
C                                                                       MDPF7510
C        PURPOSE                                                        MDPF7520
C           COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL SYMMETRIC    MDPF7530
C           MATRIX                                                      MDPF7540
C                                                                       MDPF7550
C        USAGE                                                          MDPF7560
C           CALL EIGEN(A,R,N,MV)                                        MDPF7570
C                                                                       MDPF7580
C        DESCRIPTION OF PARAMETERS                                      MDPF7590
C           A - ORIGINAL MATRIX (SYMMETRIC), DESTROYED IN COMPUTATION.  MDPF7600
C               RESULTANT EIGENVALUES ARE DEVELOPED IN DIAGONAL OF      MDPF7610
C               MATRIX A IN DECENDING ORDER.                            MDPF7620
C           R - RESULTANT MATRIX OF EIGENVECTORS (STORED COLUMNWISE,    MDPF7630
C               IN SAME SEQUENCE AS EIGENVALUES)                        MDPF7640
C           N - ORDER OF MATRICES A AND R                               MDPF7650
C           MV- INPUT CODE                                              MDPF7660
C                   0   COMPUTE EIGENVALUES AND EIGENVECTORS            MDPF7670
C                   1   COMPUTE EIGENVALUES ONLY (R NEED NOT BE         MDPF7680
C                       DIMENSIONED BUT MUST STILL APPEAR IN CALLING    MDPF7690
C                       SEQUENCE)                                       MDPF7700
C                                                                       MDPF7710
C        REMARKS                                                        MDPF7720
C           ORIGINAL MATRIX A MUST BE REAL SYMMETRIC (STORAGE MODE=1)   MDPF7730
C           MATRIX A CANNOT BE IN THE SAME LOCATION AS MATRIX R         MDPF7740
C                                                                       MDPF7750
C        SUBROUTINES USED BY THIS SUBROUTINE                            MDPF7760
C           NONE                                                        MDPF7770
C                                                                       MDPF7780
C        METHOD                                                         MDPF7790
C           DIAGONALIZATION METHOD ORIGINATED BY JACOBI AND ADAPTED     MDPF7800
C           BY VON NEUMANN FOR LARGE COMPUTERS AS FOUND IN MATHEMATICAL MDPF7810
C           METHODS FOR DIGITAL COMPUTERS, EDITED BY A. RALSTON AND     MDPF7820
C           H.S. WILF, JOHN WILEY AND SONS, NEW YORK, 1962, CHAPTER 7   MDPF7830
C                                                                       MDPF7840
C     ..................................................................MDPF7850
C                                                                       MDPF7860
C                                                                       MDPF7870
      DIMENSION A(1830), R(1830)                                        MDPF7880
C        GENERATE IDENTITY MATRIX                                       MDPF7890
C                                                                       MDPF7900
      IF(MV-1) 22,14,22                                                 MDPF7910
   22 IQ=-N                                                             MDPF7920
      DO 1 J=1,N                                                        MDPF7930
      IQ=IQ+N                                                           MDPF7940
      DO 1 I=1,N                                                        MDPF7950
      IJ=IQ+I                                                           MDPF7960
      R(IJ)=0.0                                                         MDPF7970
      IF(I-J) 1,31,1                                                    MDPF7980
   31 R(IJ)=1.0                                                         MDPF7990
    1 CONTINUE                                                          MDPF8000
C                                                                       MDPF8010
C        COMPUTE INITIAL AND FINAL NORMS (ANORM AND ANORMX)             MDPF8020
C                                                                       MDPF8030
   14 ANORM=0.0                                                         MDPF8040
      DO 2 I=1,N                                                        MDPF8050
      DO 2 J=I,N                                                        MDPF8060
      IF(I-J) 15,2,15                                                   MDPF8070
   15 IA=I+(J*J-J)/2                                                    MDPF8080
      ANORM=ANORM+A(IA)*A(IA)                                           MDPF8090
    2 CONTINUE                                                          MDPF8100
      IF(ANORM) 39, 39, 50                                              MDPF8110
   50 ANORM=2.0*SQRT(ANORM)                                             MDPF8120
      ANRMX =ANORM*1.0E-8                                               MDPF8130
C                                                                       MDPF8140
C        INITIALIZE INDICATORS AND COMPUTE THRESHOLD, THR               MDPF8150
C                                                                       MDPF8160
      IND=0                                                             MDPF8170
      THR=ANORM                                                         MDPF8180
    3 THR=THR/FLOAT(N)                                                  MDPF8190
   13 L=1                                                               MDPF8200
    4 M=L+1                                                             MDPF8210
C                                                                       MDPF8220
C        COMPUTE SIN AND COS                                            MDPF8230
C                                                                       MDPF8240
    5 MQ=(M*M-M)/2                                                      MDPF8250
      LQ=(L*L-L)/2                                                      MDPF8260
      LM=L+MQ                                                           MDPF8270
      IF(ABS(A(LM))-THR) 7,32,32                                        MDPF8280
   32 IND=1                                                             MDPF8290
      LL=L+LQ                                                           MDPF8300
      MM=M+MQ                                                           MDPF8310
      X=0.5*(A(LL)-A(MM))                                               MDPF8320
      Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X)                                    MDPF8330
      IF(X) 33,34,34                                                    MDPF8340
   33 Y=-Y                                                              MDPF8350
   34 SINX=Y/SQRT(2.0*(1.0+(SQRT(1.0-Y*Y))))                            MDPF8360
      SINX2=SINX*SINX                                                   MDPF8370
      COSX=SQRT(1.0-SINX2)                                              MDPF8380
      COSX2=COSX*COSX                                                   MDPF8390
      SINCS =SINX*COSX                                                  MDPF8400
C                                                                       MDPF8410
C        ROTATE L AND M COLUMNS                                         MDPF8420
C                                                                       MDPF8430
      ILQ=N*(L-1)                                                       MDPF8440
      IMQ=N*(M-1)                                                       MDPF8450
      DO 6 I=1,N                                                        MDPF8460
      IQ=(I*I-I)/2                                                      MDPF8470
      IF(I-L) 16,11,16                                                  MDPF8480
   16 IF(I-M) 17,11,18                                                  MDPF8490
   17 IM=I+MQ                                                           MDPF8500
      GO TO 19                                                          MDPF8510
   18 IM=M+IQ                                                           MDPF8520
   19 IF(I-L) 37,23,23                                                  MDPF8530
   37 IL=I+LQ                                                           MDPF8540
      GO TO 24                                                          MDPF8550
   23 IL=L+IQ                                                           MDPF8560
   24 X=A(IL)*COSX-A(IM)*SINX                                           MDPF8570
      A(IM)=A(IL)*SINX+A(IM)*COSX                                       MDPF8580
      A(IL)=X                                                           MDPF8590
   11 IF(MV-1) 26,6,26                                                  MDPF8600
   26 ILR=ILQ+I                                                         MDPF8610
      IMR=IMQ+I                                                         MDPF8620
      X=R(ILR)*COSX-R(IMR)*SINX                                         MDPF8630
      R(IMR)=R(ILR)*SINX+R(IMR)*COSX                                    MDPF8640
      R(ILR)=X                                                          MDPF8650
    6 CONTINUE                                                          MDPF8660
      X=2.0*A(LM)*SINCS                                                 MDPF8670
      Y=A(LL)*COSX2+A(MM)*SINX2-X                                       MDPF8680
      X=A(LL)*SINX2+A(MM)*COSX2+X                                       MDPF8690
      A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2)                     MDPF8700
      A(LL)=Y                                                           MDPF8710
      A(MM)=X                                                           MDPF8720
C                                                                       MDPF8730
C        TESTS FOR COMPLETION                                           MDPF8740
C                                                                       MDPF8750
C        TEST FOR M = LAST COLUMN                                       MDPF8760
    7 IF(M-N) 35,8,35                                                   MDPF8770
   35 M=M+1                                                             MDPF8780
      GO TO 5                                                           MDPF8790
C        TEST FOR L = SECOND FROM LAST COLUMN                           MDPF8800
    8 IF(L-(N-1)) 36,9,36                                               MDPF8810
   36 L=L+1                                                             MDPF8820
      GO TO 4                                                           MDPF8830
    9 IF(IND-1) 10,38,10                                                MDPF8840
   38 IND=0                                                             MDPF8850
      GO TO 13                                                          MDPF8860
C        COMPARE THRESHOLD WITH FINAL NORM                              MDPF8870
   10 IF(THR-ANRMX ) 39,39,3                                            MDPF8880
C                                                                       MDPF8890
C        SORT EIGENVALUES AND EIGENVECTORS                              MDPF8900
C                                                                       MDPF8910
   39 IQ=-N                                                             MDPF8920
      DO 30 I=1,N                                                       MDPF8930
      IQ=IQ+N                                                           MDPF8940
      LL=I+(I*I-I)/2                                                    MDPF8950
      JQ=N*(I-2)                                                        MDPF8960
      DO 30 J=I,N                                                       MDPF8970
      JQ=JQ+N                                                           MDPF8980
      MM=J+(J*J-J)/2                                                    MDPF8990
      IF(A(LL)-A(MM)) 40,30,30                                          MDPF9000
   40 X=A(LL)                                                           MDPF9010
      A(LL)=A(MM)                                                       MDPF9020
      A(MM)=X                                                           MDPF9030
      IF(MV-1) 28,30,28                                                 MDPF9040
   28 DO 20 K=1,N                                                       MDPF9050
      ILR=IQ+K                                                          MDPF9060
      IMR=JQ+K                                                          MDPF9070
      X=R(ILR)                                                          MDPF9080
      R(ILR)=R(IMR)                                                     MDPF9090
   20 R(IMR)=X                                                          MDPF9100
   30 CONTINUE                                                          MDPF9110
      RETURN                                                            MDPF9120
      END                                                               MDPF9130
      SUBROUTINE PLOTX(A,B,N)
      DIMENSION A(1),B(1)
      CALL GRAPH(+1.2,-1.2,1.0,-1.0)
      CALL AXES
      CALL LOC2(A,B,N,1)
      CALL PLOT
      CALL IDENT
      RETURN
      END
      SUBROUTINE GRAPH(XMX,XMN,YMX,YMN)                                 GRPH  10
      DIMENSION Z(1),X(1),Y(1),SMALL(13),CHAR(50),HOLL(11)              GRPH  20
      REAL ITEM(55,121)                                                 GRPH  30
      DATA CHAR     /'1','2','3','4','5','6','7','8','9','A','B','C',   GRPH  50
     1 'D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S', GRPH  60
     2 'T','U','V','W','X','Y','Z','+','/','=','*','&','$','@','%','?', GRPH  70
     3 '<','(',')','"',';',' '/                                         GRPH  80
      DATA HOLL     /' ','X','2','3','4','5','6','7','8','9','M'/       GRPH  90
      DATA BLANK,DD,PLUS,AIE,AMINUS,ORIG,AMULT,GT50/' ','.','+','\',    GRPH 100
     1 '-','0','#','>'/                                                 GRPH 110
      DO 115 I = 1,55                                                   GRPH 120
      DO 115 J = 1,121                                                  GRPH 130
  115 ITEM(I,J) = BLANK                                                 GRPH 140
      XMAX = XMX                                                        GRPH 150
      XMIN = XMN                                                        GRPH 160
      YMAX = YMX                                                        GRPH 170
      YMIN = YMN                                                        GRPH 180
      RETURN                                                            GRPH 190
      ENTRY LIMIT(Z,NP,IS)                                              GRPH 200
      ZMAX = -1.0E9                                                     GRPH 210
      ZMIN = +1.0E9                                                     GRPH 220
      DO 10 I = 1,NP                                                    GRPH 230
      IF(Z(I) .GT. ZMAX) ZMAX = Z(I)                                    GRPH 240
   10 IF(Z(I) .LT. ZMIN) ZMIN = Z(I)                                    GRPH 250
      IF(IS .EQ. 2) GO TO 11                                            GRPH 260
      IF(ZMAX .GT. XMAX) XMAX = ZMAX                                    GRPH 270
      IF(ZMIN .LT. XMIN) XMIN = ZMIN                                    GRPH 280
      GO TO 12                                                          GRPH 290
   11 IF(ZMAX .GT. YMAX) YMAX = ZMAX                                    GRPH 300
      IF(ZMIN .LT. YMIN) YMIN = ZMIN                                    GRPH 310
   12 RETURN                                                            GRPH 320
      ENTRY SCALE(IP)                                                   GRPH 330
      XEXP = 0.05 * (XMAX - XMIN)                                       GRPH 340
      XMAX = XMAX + XEXP                                                GRPH 350
      XMIN = XMIN - XEXP                                                GRPH 360
      YEXP = 0.05 * (YMAX - YMIN)                                       GRPH 370
      YMAX = YMAX + YEXP                                                GRPH 380
      YMIN = YMIN - YEXP                                                GRPH 390
      IF(IP .NE. 1) RETURN                                              GRPH 400
      IF(YMAX - XMAX) 20,21,22                                          GRPH 410
   20 YMAX = XMAX                                                       GRPH 420
      GO TO 21                                                          GRPH 430
   22 XMAX = YMAX                                                       GRPH 440
   21 CONTINUE                                                          GRPH 450
      IF(YMIN - XMIN) 23,24,25                                          GRPH 460
   23 XMIN = YMIN                                                       GRPH 470
      GO TO 24                                                          GRPH 480
   25 YMIN = XMIN                                                       GRPH 490
   24 EXP = 0.1667 * (XMAX - XMIN)                                      GRPH 500
      XMAX = XMAX + EXP                                                 GRPH 510
      XMIN = XMIN - EXP                                                 GRPH 520
      RETURN                                                            GRPH 530
      ENTRY AXES                                                        GRPH 540
      DELX = (XMAX - XMIN)/120.0                                        GRPH 550
      DELY = (YMAX - YMIN)/54.0                                         GRPH 560
      K = 0                                                             GRPH 570
      M = 0                                                             GRPH 580
      IF(XMIN .LT. 0.0 .AND. XMAX .GT. 0.0) GO TO 30                    GRPH 590
      GO TO 31                                                          GRPH 600
   30 K = (-XMIN/DELX) + 1.5                                            GRPH 610
      DO 32 I = 1,55                                                    GRPH 620
   32 ITEM(I,K) = AIE                                                   GRPH 630
   31 CONTINUE                                                          GRPH 640
      IF(YMIN .LT. 0.0 .AND. YMAX .GT. 0.0) GO TO 33                    GRPH 650
      GO TO 34                                                          GRPH 660
   33 M = (YMAX/DELY) + 1.5                                             GRPH 670
      DO 35 J = 1,121                                                   GRPH 680
   35 ITEM(M,J) = AMINUS                                                GRPH 690
      IF(K .EQ. 0 .OR. M .EQ. 0) GO TO 34                               GRPH 700
      ITEM(M,K) = ORIG                                                  GRPH 710
   34 RETURN                                                            GRPH 720
      ENTRY LOC1(X,Y,NPOI)                                              GRPH 730
      DELX = (XMAX - XMIN)/120.0                                        GRPH 740
      DELY = (YMAX - YMIN)/54.0                                         GRPH 750
      DO 50 II = 1,NPOI                                                 GRPH 760
      I = (YMAX - Y(II))/DELY + 1.5                                     GRPH 770
      J = (X(II) - XMIN)/DELX + 1.5                                     GRPH 780
      IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 50              GRPH 790
      DO 52 JJ = 1,10                                                   GRPH 800
      IF(ITEM(I,J) .EQ. HOLL(JJ)) GO TO 51                              GRPH 810
   52 CONTINUE                                                          GRPH 820
      IF(ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS.OR.ITEM(I,J).EQ.ORIG)  GRPH 830
     1  ITEM(I,J) = HOLL(2)                                             GRPH 840
      GO TO 50                                                          GRPH 850
   51 ITEM(I,J) = HOLL(JJ+1)                                            GRPH 860
   50 CONTINUE                                                          GRPH 870
      RETURN                                                            GRPH 880
      ENTRY LOC2(X,Y,NPOI,ISTART)                                       GRPH 890
      DELX = (XMAX - XMIN)/120.0                                        GRPH 900
      DELY = (YMAX - YMIN)/54.0                                         GRPH 910
      ISPN = ISTART + NPOI - 1                                          GRPH 920
      DO 53 II = ISTART,ISPN                                            GRPH 930
      I = (YMAX - Y(II))/DELY + 1.5                                     GRPH 940
      J = (X(II) - XMIN)/DELX + 1.5                                     GRPH 950
      IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 53              GRPH 960
      IF(ITEM(I,J).EQ.BLANK.OR.ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS  GRPH 970
     1 .OR.ITEM(I,J).EQ.ORIG) GO TO 54                                  GRPH 980
      ITEM(I,J) = AMULT                                                 GRPH 990
      GO TO 53                                                          GRPH1000
   54 CONTINUE                                                          GRPH1010
      IF(II .GT. 50) GO TO 55                                           GRPH1020
      ITEM(I,J) = CHAR(II)                                              GRPH1030
      GO TO 53                                                          GRPH1040
   55 ITEM(I,J) = GT50                                                  GRPH1050
   53 CONTINUE                                                          GRPH1060
      RETURN                                                            GRPH1070
      ENTRY LOC3(X,Y,NPOI,POINT)                                        GRPH1080
      DELX = (XMAX - XMIN)/120.0                                        GRPH1090
      DELY = (YMAX - YMIN)/54.0                                         GRPH1100
      DO 56 II = 1,NPOI                                                 GRPH1110
      I = (YMAX - Y(II))/DELY + 1.5                                     GRPH1120
      J = (X(II) - XMIN)/DELX + 1.5                                     GRPH1130
      IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1)  GO TO 56             GRPH1140
      IF(ITEM(I,J).EQ.BLANK.OR.ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS  GRPH1150
     1 .OR. ITEM(I,J) .EQ. ORIG .OR. ITEM(I,J) .EQ. POINT) GO TO 57     GRPH1160
      ITEM(I,J) = AMULT                                                 GRPH1170
      GO TO 56                                                          GRPH1180
   57 ITEM(I,J) = POINT                                                 GRPH1190
   56 CONTINUE                                                          GRPH1200
      RETURN                                                            GRPH1210
      ENTRY PLOT                                                        GRPH1220
      SMALL(1) = XMIN                                                   GRPH1230
      DO 120 I = 1,12                                                   GRPH1240
  120 SMALL(I+1) = SMALL(I) + 10.0 * DELX                               GRPH1250
      VALUE = YMAX + DELY                                               GRPH1260
      WRITE(6,301)                                                      GRPH1270
  301 FORMAT(9X,' +.........+.........+.........+.........+.........+...GRPH1280
     1......+.........+.........+.........+.........+.........+.........GRPH1290
     2+')                                                               GRPH1300
      DO 330 I = 1,55                                                   GRPH1310
      VALUE = VALUE - DELY                                              GRPH1320
      II = I - 1                                                        GRPH1330
      L = II + 6                                                        GRPH1340
      IF(L/6*6 - L) 320,310,320                                         GRPH1350
  310 A = PLUS                                                          GRPH1360
      WRITE(6,302) VALUE,A,(ITEM(I,J), J=1,121),A                       GRPH1370
  302 FORMAT(1X,F8.3,123A1)                                             GRPH1380
      GO TO 330                                                         GRPH1390
  320 A = DD                                                            GRPH1400
      WRITE(6,303)  A,(ITEM(I,J), J=1,121),A                            GRPH1410
  303 FORMAT(9X,123A1)                                                  GRPH1420
  330 CONTINUE                                                          GRPH1430
      WRITE(6,301)                                                      GRPH1440
      WRITE(6,304)  (SMALL(K), K=1,12)                                  GRPH1450
  304 FORMAT(4X,12F10.3)                                                GRPH1460
      WRITE(6,305) SMALL(13)                                            GRPH1470
  305 FORMAT(1H+,123X,F8.2/)                                            GRPH1480
      RETURN                                                            GRPH1490
      ENTRY IDENT                                                       GRPH1500
      WRITE(6,400) (J, J=1,25)                                          GRPH1510
  400 FORMAT(/35X,'*****IDENTIFICATION KEY FOR PLOTS WITH IDENTIFIED POIGRPH1520
     1NTS*****'//' PT #',25I5)                                          GRPH1530
      WRITE(6,401)  (CHAR(I), I=1,25)                                   GRPH1540
  401 FORMAT(' CHAR',25(4X,A1)/)                                        GRPH1550
      WRITE(6,402)  (J, J=26,50)                                        GRPH1560
  402 FORMAT(' PT #', 25I5)                                             GRPH1570
      WRITE(6,403) (CHAR(I), I=26,50)                                   GRPH1580
  403 FORMAT(' CHAR', 25(4X,A1)//)                                      GRPH1590
      WRITE(6,404) GT50, AMULT                                          GRPH1600
  404 FORMAT(' POINT NUMBERS ABOVE 50 IDENTIFIED AS  ',A1/' MULTIPLE POIGRPH1610
     1NTS IDENTIFIED AS  ',A1/)                                         GRPH1620
      RETURN                                                            GRPH1630
      END                                                               GRPH1640