C sindscal.f
C The author of this software is Sandra Pruzansky.

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 SECOND MDS Package of AT&T Bell Laboratories

C For explanation of the method this software implements, see
C Arabie,P., Carroll,J.D., & DeSarbo,W.S. (1987). Three-Way Scaling and
C Clustering. Newbury Park, CA: Sage, and

C Carroll,J.D. & Chang,J.J. (1970), "Analysis of individual differences
C in multidimensional scaling via an N-way generalization of
C 'Eckart-Young' decomposition" in Psychometrika, 35,283-319, and

C Carroll,J.D. (1972). Individual differences and multidimensional
C scaling, in R.N.Shepard, A.K.Romney & S.B. Nerlove (Eds.),
C Multidimensional Scaling: Theory and Applications in the Behavioral
C Sciences (Vol.1, pp.105-155). New York and London: Seminar Press.

C The latter two have been reprinted in
C P. Davies & A.P.M. Coxon (Eds.), (1984)
C "Key Texts in Multidimensional Scaling" Exeter, NH: Heinemann.
C First one: pp.  229-252;  second: pp. 267-301.
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@
C          MACHINE DEPENDENT CONSTANT, XMAXN,  AT LINE   XA00  12
C            CURRENTLY SET FOR THE VAX.
C            REPRESENTS MINIMUM OF LARGEST POSITIVE MACHINE NUMBER
C            AND MINUS LARGEST NEGATIVE MACHINE NUMBER
C          MODIFIED MAY, 1987; CHANGE IN RDAT
C            RDAT  44 COMMENTED OUT TO AVOID LOADER ERROR;
C            IF A USER-WRITTEN MYREAD SUBROUTINE IS
C            TO BE USED, REMOVE C FROM FIRST COLUMN
C            AND RECOMPILE
C          MODIFIED MAY, 1982; CHANGE IN CANDE
C             ADDED RELAXATION FACTOR TO SPEED UP CONVERGENCE;
C                CHANGED CONVERGENCE CRIT TO 1.E-6 FROM 1.E-5
C                 THIS WILL INCREASE CONSIDERABLY, THE NUMBER OF
C                 ITERATIONS FOR CONVERGENCE BUT WILL REDUCE THE
C                 POSSIBILITY OF LOCAL MINIMA
C          MODIFIED OCT., 1977; CHANGE IN SERCH                         X000   3
C          MODIFIED AUG., 1977; CHANGES IN SORT, RUNIFV, SMTINV         X000   4
C  WRITTEN BY SANDRA PRUZANSKY                                          X000   5
C             BELL TELEPHONE LABORATORIES                               X000   6
C             MURRAY HILL, NEW JERSEY                                   X000   7
C                                                                       X000   8
C   BASED ON ALGORITHM BY CARROLL AND CHANG  PUBLISHED IN               X000   9
C     PSYCHOMETRIKA,35:1970,283-319.                                    X000  10
C                                                                       X000  11
C                                                                       X000  12
C     PARAMETER DEFINITIONS                                             X000  13
C***************************************                                X000  14
C     MAXDIM- MAXIMUM NUMBER OF DIMENSIONS   (MAX. =10)                 X000  15
C     MINDIM- MINIMUM NUMBER OF DIMENSIONS                              X000  16
C     NWT(1)  - NUMBER OF MATRICES                                      X000  17
C     NWT(2)- NUMBER OF STIMULI                                         X000  18
C     IRDATA- 1, READ IN LOWERHALF OF SIMILARITY MATRIX WITHOUT DIAGS.  X000  19
C     IRDATA- 2, READ IN DISSIMILARITIES, LOWERHALF WITHOUT DIAGS.      X000  20
C     IRDATA- 3, READ IN EUCLIDEAN DISTANCES, LOWERHALF WITHOUT DIAGS.  X000  21
C     IRDATA- 4, READ IN LOWERHALF CORRELATION MATRIX WITHOUT DIAGS.    X000  22
C            -4, SAME AS ABOVE BUT PROGRAM WILL NORMALIZE SO SSQ=1.     X000  23
C     IRDATA- 5, READ IN LOWERHALF COVARIANCE MATRIX WITH DIAGS.        X000  24
C            -5, SAME AS ABOVE BUT PROGRAM WILL NORMALIZE SO SSQ=1.     X000  25
C     IRDATA- 6, READ IN FULL SYMMETRIC SIMILARITY MATRIX, DIAG. IGNOREDX000  26
C     IRDATA- 7, READ IN FULL SYMMETRIC DISSIM. MATRIX,DIAG. IGNORED    X000  27
C     ITMAX - MAXIMUM NUMBER OF ITERATIONS.                             X000  28
C     IF ITMAX=0 PROGRAM WILL SOLVE FOR WEIGHTS ONLY.                   X000  29
C  PUNCH OPTIONS - ALL PUNCH OPTIONS INCLUDE PUNCHING NORMALIZED        X000  30
C                  SOLUTION EXCEPT -1.                                  X000  31
C     IPUNSP  -1, NO PUNCHED OUTPUT                                     X000  32
C     IPUNSP - 0, PUNCH NORMALIZED  SOLUTION                            X000  33
C     IPUNSP - 1, PUNCH SCALAR PRODUCT MATRICES.                        X000  34
C     IPUNSP - 2, PUNCH UNNORMALIZED MATRICES - A1,A2,A3                X000  35
C     IPUNSP - 3, PUNCH ALL OF THE ABOVE                                X000  36
C      IPLOT - -1, NO PLOT                                              X000  37
C      IPLOT - 0, PLOT WITH POINTS                                      X000  38
C     IRN - 0, READ IN INITIAL CONFIGURATION.                           X000  39
C     IRN - AN INTEGER FOR GENERATING A RANDOM STARTING                 X000  40
C           CONFIGURATION. (MAXIMUM 4 DIGITS)                           X000  41
C     ARRAYS                                                            X000  42
C     Y - DATA.  Y IS A 2 SUBSCRIPTED VARIABLE.                         X000  43
C                Y(N1,LSD), WHERE N1 IS NO. OF MATRICES                 X000  44
C                AND LSD IS (N2(N2+1))/2;(N2=# OF STIM)                 X000  45
C     A - COORDINATE MATRIX.                                            X000  46
C         A(I,J,K), WHERE I - ELEMENTS WITHIN MATRICES                  X000  47
C                         J - DIMENSIONS.                               X000  48
C                         K - MATRIX.                                   X000  49
C      AA - SCRATCH MATRIX                                              X000  50
C          AA(I,J) - WHERE I - MAXIMUM OF # STIM AND # OF MATRICES      X000  51
C                          J - MAXIMUM # OF DIMENSIONS                  X000  52
C     TIT - TITLE, USE NO MORE THAN ONE CARD.                           X000  53
C                                                                       X000  54
C ----- DATA CARD INPUT FORMAT: -----                                   X000  55
C                                                                       X000  56
C     FIRST CARD (4I4): MAXDIM,MINDIM,NWT(1),NWT(2)                     X000  57
C     SECOND CARD (5I4): ITMAX,IRDATA,IPUNSP,IPLOT,IRN                  X000  58
C     THIRD CARD (80A1) : TITLE CARD FOR THIS DATA SET                  X000  59
C     FOURTH CARD (80A1) : (DATA FORMAT ENCLOSED IN PARENTHESES)        X000  60
C                                    OR                                 X000  61
C                          MYREAD (IN COLUMNS 1-6)                      X000  62
C     FIFTH TO NTH CARD : DATA IN FORMAT GIVEN ON FOURTH CARD           X000  63
C     OPTIONAL CONFIGURATION CARDS - (IF IRN=0)                         X000  64
C       FIRST CONFIG. CARD - (80A1)  :FORMAT FOR STARTING CONFIGURATION X000  65
C        SECOND THRU LAST CONFIG. CARD  : STARTING STIMULUS             X000  66
C                            CONFIGURATION IN FORMAT SPECIFIED BY       X000  67
C                            PRECEDING CARD.                            X000  68
C     LAST CARD:  MUST BE BLANK IN COLUMNS 1-20                         X000  69
C*************************                                              X000  70
C                                                                       X000  71
C                                                                       X000  72
C        MAIN PROGRAM FOR SINDSCAL - (FOR PORTABLE USE)                 X000  73
C              WRITTEN BY S. PRUZANSKY                                  X000  74
C            LATEST CHANGE 9-16-75                                      X000  75
      DIMENSION ARRAY(5000)                                             X000  76
      COMMON /PARAM/ MAXDIM, MINDIM, ITMAX, IRDATA, N, IPUNSP,          X000  77
     *    IPLOT, IRN, CRIT, NWT(3)                                      X000  78
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             X000  79
      DATA MAXCOR /5000/                                                X000  80
C     DESIGNATE DEVICE NUMBERS FOR INPUT, OUTPUT                        X000  81
      LREAD = 5                                                         X000  82
      LPRINT = 6                                                        X000  83
      LPUNCH = 7                                                        X000  84
C  READ INPUT PARMETERS NEEDED FOR STORAGE ALLOCATION                   X000  85
  10  READ (LREAD,9999) MAXDIM, MINDIM, NWT(1), NWT(2)                  X000  86
      IF (MAXDIM.EQ.0) STOP                                             X000  87
C  CHECK VALIDITY OF PARAMETER CARD                                     X000  88
      IF (MAXDIM.GT.10) GO TO 20                                        X000  89
      IF (MINDIM.GT.10) GO TO 20                                        X000  90
9999  FORMAT (4I4)                                                      X000  91
      READ (LREAD,9998) ITMAX, IRDATA, IPUNSP, IPLOT, IRN               X000  92
9998  FORMAT (5I4)                                                      X000  93
C  COMPUTE AMOUNT OF CORE NEEDED FOR ARRAYS                             X000  94
      LY = NWT(1)*((NWT(2)*(NWT(2)+1))/2)                               X000  95
      N = MAX0(NWT(1),NWT(2))                                           X000  96
      LA = MAX0(N*MAXDIM*3,NWT(2)**2)                                   X000  97
      LAA = N*MAXDIM                                                    X000  98
      MEMORY = LY + LA + LAA                                            X000  99
      IF (MEMORY.GT.MAXCOR) GO TO 30                                    X000 100
C  COMPUTE STARTING LOCATIONS FOR THREE STORAGE ARRAYS                  X000 101
      LOCY = 1                                                          X000 102
      LOCA = LOCY + LY                                                  X000 103
      LOCAA = LOCA + LA                                                 X000 104
      CRIT = 1.E-6                                                      X000 105
      CALL SINSCL(ARRAY(LOCY), LY, ARRAY(LOCA), LA,                     X000 106
     *    ARRAY(LOCAA), LAA)                                            X000 107
      GO TO 10                                                          X000 108
  20  WRITE (LPRINT,9997)                                               X000 109
9997  FORMAT (45H0REQUEST EXCEEDS MAXIMUM NUMBER OF DIMENSIONS/         X000 110
     *    23H ***CHECK DECK SETUP***)                                   X000 111
      STOP                                                              X000 112
  30  WRITE (LPRINT,9996)                                               X000 113
9996  FORMAT (46H0CORE REQUESTED EXCEEDS LIMITS OF DIMENSIONED ,        X000 114
     *    5HARRAY, 32H ***CHECK VARIABLES LY,LA,LAA***)                 X000 115
      STOP                                                              X000 116
      END                                                               X000 117
C                                                                       SINS   1
      SUBROUTINE SINSCL(Y, LY, A, LA, AA, LAA)                          SINS   2
C SINSCL       FOR SINDSCAL              SEPTEMBER, 1975                SINS   3
C WRITTEN BY S. PRUZANSKY                                               SINS   4
C         CONTROL PROGRAM FOR INDSCAL ANALYSIS                          SINS   5
      COMMON /PARAM/ MAXDIM, MINDIM, ITMAX, IRDATA, MNWT,               SINS   6
     *    IPUNSP, IPLOT, IRN, CRIT, NWT(3)                              SINS   7
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             SINS   8
      DIMENSION Y(LY), A(LA), AA(LAA)                                   SINS   9
      INTEGER TIT(80)                                                   SINS  10
C  READ TITLE                                                           SINS  11
      READ (LREAD,9999) (TIT(I),I=1,80)                                 SINS  12
C  SET PARAMETERS                                                       SINS  13
      NA = 2                                                            SINS  14
      NWT(3) = NWT(2)                                                   SINS  15
      N = 3                                                             SINS  16
      N1 = NWT(1)                                                       SINS  17
      N2 = NWT(2)                                                       SINS  18
      LSD = (N2*(N2+1))/2                                               SINS  19
      NDIM = 0                                                          SINS  20
C  DO INDSCAL ANALYSES IN MAXDIM TO MINDIM DIMENSIONS                   SINS  21
      MMDIM = MAXDIM + MINDIM                                           SINS  22
      DO 60 IDIM=MINDIM,MAXDIM                                          SINS  23
        NF = MMDIM - IDIM                                               SINS  24
        MAXIT = ITMAX                                                   SINS  25
        IF (ITMAX.EQ.0) GO TO 10                                        SINS  26
        ISTI = 0                                                        SINS  27
        GO TO 20                                                        SINS  28
  10    ISTI = 1                                                        SINS  29
  20    IF (IDIM.GT.MINDIM) WRITE (LPRINT,9992)                         SINS  30
C  PRINT TITLE AND PARAMETER INFO.                                      SINS  31
        WRITE (LPRINT,9997)                                             SINS  32
        WRITE (LPRINT,9998) (TIT(I),I=1,80)                             SINS  33
        WRITE (LPRINT,9995)                                             SINS  34
        WRITE (LPRINT,9994)                                             SINS  35
        WRITE (LPRINT,9996) NF, IRDATA, MAXIT, IPUNSP, IPLOT,           SINS  36
     *      IRN                                                         SINS  37
        WRITE (LPRINT,9993) (NWT(I),I=1,NA)                             SINS  38
        WRITE (LPRINT,9995)                                             SINS  39
        IF (ISTI.NE.0) MAXIT = 0                                        SINS  40
        JRDATA = IABS(IRDATA)                                           SINS  41
        IF (NDIM.EQ.0) GO TO 30                                         SINS  42
C  REDUCE DIMENSIONALITY OF STARTING CONFIG.                            SINS  43
        NFO = NF + 1                                                    SINS  44
        NQ = NFO*N                                                      SINS  45
        CALL DIMA(A, MNWT, NQ, NFO, NF)                                 SINS  46
        GO TO 40                                                        SINS  47
C  READ AND PREPROCESS INPUT DATA                                       SINS  48
  30    CALL RDATA(Y, N1, N2, LSD, JRDATA, A, IPUNSP)                   SINS  49
C  PERFORM CANONICAL DECOMPOSITION                                      SINS  50
  40    CALL CANDE(Y, N1, N2, LSD, NA, A, MNWT, NF, N, NWT,             SINS  51
     *      IRDATA, ISTI, MAXIT, CRIT, IRN, NDIM, FVAF, AA,             SINS  52
     *      IPUNSP)                                                     SINS  53
C   NORMALIZE FINAL WEIGHTS AND STIMULUS MATRICES                       SINS  54
        CALL NORMA(A, AA, MNWT, NF, N1, N2, NA, FVAF, AA,               SINS  55
     *      IPUNSP)                                                     SINS  56
C  COMPUTE CORRELATIONS BETWEEN SCALAR PROD. AND SOLUTION FOR EACH      SINS  57
C   DATA MATRIX.                                                        SINS  58
        CALL SUBR(A, MNWT, NF, N, Y, N1, N2, LSD)                       SINS  59
        IF (IPLOT.LT.0) GO TO 50                                        SINS  60
C  PLOT FINAL WEIGHTS AND STIMULUS SPACES                               SINS  61
        CALL PPLOT(A, MNWT, NF, TIT, NWT)                               SINS  62
  50    NDIM = 1                                                        SINS  63
  60  CONTINUE                                                          SINS  64
9999  FORMAT (80A1)                                                     SINS  65
9998  FORMAT (1H0, 80A1)                                                SINS  66
9997  FORMAT (1H0, 23X, 17HSYMMETRIC INDSCAL)                           SINS  67
9996  FORMAT (1H , I3, 4I7, 2X, I7)                                     SINS  68
9995  FORMAT (46H0*********************************************,        SINS  69
     *    5H*****)                                                      SINS  70
9994  FORMAT (11H PARAMETERS/20H  DIM  IRDATA  ITMAX, 7H IPUNCH,        SINS  71
     *    14H  IPLOT    IRN)                                            SINS  72
9993  FORMAT (18H NO. OF MATRICES =, I4, 16H  NO. OF STIM. =,           SINS  73
     *    I4)                                                           SINS  74
9992  FORMAT (1H1)                                                      SINS  75
      RETURN                                                            SINS  76
      END                                                               SINS  77
C                                                                       RDAT   1
      SUBROUTINE RDATA(Y, N1, N2, LSD, IRDATA, D, IPUNSP)               RDAT   2
C RDATA          FOR SINDSCAL              SEPTEMBER, 1975              RDAT   3
C WRITTEN BY S. PRUZANSKY                                               RDAT   4
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             RDAT   5
      DIMENSION Y(N1,LSD), D(N2,N2), MY(2)                              RDAT   6
      INTEGER FMT(80)                                                   RDAT   7
      DATA IC /1HI/, LP /1H)/, MY(1) /1HM/, MY(2) /1HY/                 RDAT   8
C     IRDATA- 1, READ IN LOWERHALF OF SIMILARITY MATRIX WITHOUT DIAGS.  RDAT   9
C     IRDATA- 2, READ IN DISSIMILARITIES, LOWERHALF WITHOUT DIAGS.      RDAT  10
C     IRDATA- 3, READ IN EUCLIDEAN DISTANCES, LOWERHALF WITHOUT DIAGS.  RDAT  11
C     IRDATA- 4, READ IN LOWERHALF CORRELATION MATRIX WITHOUT DIAGS.    RDAT  12
C     IRDATA- 5, READ IN LOWERHALF COVARIANCE MATRIX WITH DIAGS.        RDAT  13
C     IRDATA- 6, READ IN FULL SYMMETRIC SIMILARITY MATRIX, DIAG. IGNOREDRDAT  14
C     IRDATA- 7, READ IN FULL SYMMETRIC DISSIMILARITY MATRIX,DIAG. IGNORRDAT  15
C     IADDC=1, COMPUTE ADDITIVE CONSTANT.                               RDAT  16
C     IADDC=0, DO NOT COMPUTE ADDITIVE CONSTANT.                        RDAT  17
C     LSD - LENGTH OF SYM. DATA VECTOR                                  RDAT  18
C*********************************************                          RDAT  19
      MIP2 = MOD(IPUNSP,2)                                              RDAT  20
      IF (MIP2.LE.0) GO TO 10                                           RDAT  21
      WRITE (LPUNCH,9997)                                               RDAT  22
  10  IADDC = 1                                                         RDAT  23
      IFIRST = 2                                                        RDAT  24
      II = N2                                                           RDAT  25
      IF (IRDATA.EQ.3) IADDC = 0                                        RDAT  26
      IF (IRDATA.GE.5) IFIRST = 1                                       RDAT  27
      READ (LREAD,9999) (FMT(I),I=1,80)                                 RDAT  28
C  BYPASS FORMAT CHECKING                                               RDAT  29
9999  FORMAT (80A1)                                                     RDAT  30
      IF (FMT(1).EQ.MY(1) .AND. FMT(2).EQ.MY(2)) GO TO 30               RDAT  31
C  CHECK THAT DATA FORMAT IS PROPER                                     RDAT  32
      DO 20 I=1,80                                                      RDAT  33
        IF (FMT(I).EQ.LP) GO TO 30                                      RDAT  34
        IF (FMT(I).NE.IC) GO TO 20                                      RDAT  35
        WRITE (LPRINT,9998)                                             RDAT  36
        STOP                                                            RDAT  37
  20  CONTINUE                                                          RDAT  38
9998  FORMAT (38H***  FORMAT FOR INPUT DATA IS IMPROPER,                RDAT  39
     *    24H - MUST BE E OR F FORMAT)                                  RDAT  40
  30  DO 190 I1=1,N1                                                    RDAT  41
C  READ USERS DATA READ SUBROUTINE                                      RDAT  42
        IF (FMT(1).NE.MY(1) .OR. FMT(2).NE.MY(2)) GO TO 40              RDAT  43
C       CALL MYREAD(D, N2, I1)                                          RDAT  44
        GO TO 90                                                        RDAT  45
  40    IN = 0                                                          RDAT  46
        DO 80 J=IFIRST,N2                                               RDAT  47
          IF (IRDATA-5) 50, 60, 70                                      RDAT  48
  50      II = J - 1                                                    RDAT  49
          GO TO 70                                                      RDAT  50
  60      II = J                                                        RDAT  51
  70      READ (LREAD,FMT) (D(J,M),M=1,II)                              RDAT  52
  80    CONTINUE                                                        RDAT  53
  90    DO 110 J=2,N2                                                   RDAT  54
          JJ = J - 1                                                    RDAT  55
          DO 100 M=1,JJ                                                 RDAT  56
            IF (IRDATA.EQ.1 .OR. IRDATA.EQ.6) D(J,M) = -D(J,M)          RDAT  57
            D(M,J) = D(J,M)                                             RDAT  58
 100      CONTINUE                                                      RDAT  59
 110    CONTINUE                                                        RDAT  60
        GO TO (120, 120, 120, 140, 160, 120, 120), IRDATA               RDAT  61
C     COMPUTE ADDITIVE CONSTANT AND SCALAR PRODUCT MATRICES.            RDAT  62
 120    CALL TORGB(D, N2, IADDC)                                        RDAT  63
        IF (MIP2.LE.0) GO TO 160                                        RDAT  64
C  PUNCH LOWER-HALF SCALAR PROD. MATRIX WITH DIAG                       RDAT  65
        DO 130 I=1,N2                                                   RDAT  66
          WRITE (LPUNCH,9996) I1, I, (D(I,J),J=1,I)                     RDAT  67
 130    CONTINUE                                                        RDAT  68
        GO TO 160                                                       RDAT  69
 140    DO 150 I=1,N2                                                   RDAT  70
          D(I,I) = 1.                                                   RDAT  71
 150    CONTINUE                                                        RDAT  72
C STORE LOWER HALF SCALAR PROD. MATRIX WITH DIAG IN Y                   RDAT  73
 160    K = 0                                                           RDAT  74
        DO 180 I=1,N2                                                   RDAT  75
          DO 170 J=1,I                                                  RDAT  76
            K = K + 1                                                   RDAT  77
            Y(I1,K) = D(I,J)                                            RDAT  78
 170      CONTINUE                                                      RDAT  79
 180    CONTINUE                                                        RDAT  80
 190  CONTINUE                                                          RDAT  81
9997  FORMAT (38H(6X,8F8.4)     SCALAR PRODUCT MATRICES)                RDAT  82
9996  FORMAT (2I3, 8F8.4/(F14.4, 7F8.4))                                RDAT  83
      RETURN                                                            RDAT  84
      END                                                               RDAT  85
C                                                                       TORG   1
      SUBROUTINE TORGB(A, N, IADDC)                                     TORG   2
C TORGB          FOR SINDSCAL              SEPTEMBER, 1975              TORG   3
C WRITTEN BY J.J. CHANG FOR INDSCAL                                     TORG   4
C MODIFIED FOR SINDSCAL BY S. PRUZANSKY                                 TORG   5
C       TORGERSON METHOD TO COMPUTE SCALAR PRODUCT MATRICES             TORG   6
      DIMENSION DR(200), A(N,N)                                         TORG   7
      XN = N                                                            TORG   8
      DO 10 I=1,N                                                       TORG   9
        DR(I) = 0.                                                      TORG  10
        A(I,I) = 0.                                                     TORG  11
  10  CONTINUE                                                          TORG  12
      SUMD = 0.                                                         TORG  13
      IF (IADDC) 110, 110, 20                                           TORG  14
C  COMPUTE ADDITIVE CONSTANT                                            TORG  15
  20  ADDC = -100.                                                      TORG  16
      N1 = N - 1                                                        TORG  17
      N2 = N - 2                                                        TORG  18
      DO 100 I=1,N2                                                     TORG  19
        JJ = I + 1                                                      TORG  20
        DO 90 J=JJ,N1                                                   TORG  21
          KK = J + 1                                                    TORG  22
          DO 80 K=KK,N                                                  TORG  23
            A1 = A(I,K)                                                 TORG  24
            A2 = A(J,K)                                                 TORG  25
            A3 = A(I,J)                                                 TORG  26
            DO 70 M=1,3                                                 TORG  27
              GO TO (30, 40, 50), M                                     TORG  28
  30          CIJK = A1 - A2 - A3                                       TORG  29
              GO TO 60                                                  TORG  30
  40          CIJK = A2 - A1 - A3                                       TORG  31
              GO TO 60                                                  TORG  32
  50          CIJK = A3 - A2 - A1                                       TORG  33
  60          IF (CIJK.GT.ADDC) ADDC = CIJK                             TORG  34
  70        CONTINUE                                                    TORG  35
  80      CONTINUE                                                      TORG  36
  90    CONTINUE                                                        TORG  37
 100  CONTINUE                                                          TORG  38
 110  DO 140 I=2,N                                                      TORG  39
        II = I - 1                                                      TORG  40
        DO 130 J=1,II                                                   TORG  41
          IF (IADDC.EQ.0) GO TO 120                                     TORG  42
          A(I,J) = A(I,J) + ADDC                                        TORG  43
          A(J,I) = A(I,J)                                               TORG  44
 120      SUMD = SUMD + A(I,J)**2                                       TORG  45
 130    CONTINUE                                                        TORG  46
 140  CONTINUE                                                          TORG  47
      ADIJ = 2.*SUMD/XN**2                                              TORG  48
      DO 160 I=1,N                                                      TORG  49
        DO 150 J=1,N                                                    TORG  50
          DR(I) = DR(I) + A(I,J)**2                                     TORG  51
 150    CONTINUE                                                        TORG  52
 160  CONTINUE                                                          TORG  53
C  COMPUTE SCALAR PRODUCT MATRIX                                        TORG  54
      SUM = 0.                                                          TORG  55
      SUMD = 0.                                                         TORG  56
      DO 190 I=1,N                                                      TORG  57
        DO 180 J=1,I                                                    TORG  58
          QQ = (DR(I)+DR(J))/XN                                         TORG  59
          A(I,J) = (QQ-A(I,J)**2-ADIJ)*(.5)                             TORG  60
          IF (I.EQ.J) GO TO 170                                         TORG  61
          SUM = SUM + A(I,J)**2                                         TORG  62
          GO TO 180                                                     TORG  63
 170      SUMD = SUMD + A(I,J)**2                                       TORG  64
 180    CONTINUE                                                        TORG  65
 190  CONTINUE                                                          TORG  66
C  NORMALIZE SCALAR PRODUCT MATRIX                                      TORG  67
      Q = 1./SQRT(SUMD+2.*SUM)                                          TORG  68
      DO 210 I=1,N                                                      TORG  69
        DO 200 J=1,I                                                    TORG  70
          A(I,J) = A(I,J)*Q                                             TORG  71
 200    CONTINUE                                                        TORG  72
 210  CONTINUE                                                          TORG  73
      RETURN                                                            TORG  74
      END                                                               TORG  75
C                                                                       CAND   1
      SUBROUTINE CANDE(Y, N1, N2, LSD, NA, A, MNWT, ND, N, NWT,         CAND   2
     *    IRDATA, JSTI, NAXIT, CRIT, IRN, NDIM, FVAF, S, IPUN)          CAND   3
C CANDE          FOR SINDSCAL              SEPTEMBER, 1975              CAND   4
C WRITTEN BY S. PRUZANSKY                                               CAND   5
C        LATEST CHANGE 5-20-82 BY SP                                    CAND   6
C       RELAXATION FACTOR INTRODUCED TO SPEED UP CONVERGENCE
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             CAND   7
      DIMENSION Y(N1,LSD), A(MNWT,ND,3), S(MNWT,ND), R(55)              CAND   8
      DIMENSION NWT(3), SCR(55), SUMM(2)                                CAND   9
      INTEGER FMT(80)                                                   CAND  10
      ITCOM = 0                                                         CAND  11
C  THE FOLLOWING STATEMENT WAS INSERTED 5-20-82
      ALPHA=1.5
      MAXIT = NAXIT                                                     CAND  12
      ISTI = JSTI                                                       CAND  13
      NF = ND                                                           CAND  14
      IF (NDIM.EQ.1) GO TO 40                                           CAND  15
      SSQY = N1                                                         CAND  16
C  COMPUTE SSQ FOR COVAR. OR CORR. DATA (WITH OPTION TO NORMALIZE)      CAND  17
      IAD = IABS(IRDATA)                                                CAND  18
      IF (IAD.EQ.4 .OR. IAD.EQ.5) CALL SSQC(Y, N1, LSD, S, N2,          CAND  19
     *    SSQY, IRDATA)                                                 CAND  20
      IF (IRN) 10, 20, 10                                               CAND  21
C  GENERATE RANDOM STARTING CONFIG.                                     CAND  22
  10  CALL RCON(A, MNWT, NF, N2, IRN)                                   CAND  23
      GO TO 70                                                          CAND  24
C  READ IN STARTING STIM. CONFIG.                                       CAND  25
  20  READ (LREAD,9997) (FMT(I),I=1,80)                                 CAND  26
      DO 30 K=1,NF                                                      CAND  27
        READ (LREAD,FMT) (A(J,K,2),J=1,N2)                              CAND  28
  30  CONTINUE                                                          CAND  29
C  SET TABLE 2 EQ. TABLE 3 TO START                                     CAND  30
  40  DO 60 K=1,NF                                                      CAND  31
        DO 50 J=1,N2                                                    CAND  32
          A(J,K,3) = A(J,K,2)                                           CAND  33
  50    CONTINUE                                                        CAND  34
  60  CONTINUE                                                          CAND  35
  70  WRITE (LPRINT,9991)                                               CAND  36
      DO 80 J=1,NF                                                      CAND  37
        WRITE (LPRINT,9996) J, (A(M,J,2),M=1,N2)                        CAND  38
  80  CONTINUE                                                          CAND  39
C  BEGIN FITTING CYCLE                                                  CAND  40
  90  IT = 0                                                            CAND  41
 100  DO 350 JZ=1,N                                                     CAND  42
C   THE FOLLOWING STATEMENT MODIFIED 5-20-82
        JB = JZ                                                         CAND  43
        IF (IT.NE.0 .AND. ISTI.EQ.0) GO TO 110                          CAND  44
        IF (JB.GE.2) GO TO 350                                          CAND  45
C  COMPUTE INVERSE                                                      CAND  46
 110    KL = 0                                                          CAND  47
        DO 150 K1=1,NF                                                  CAND  48
          DO 140 K2=1,K1                                                CAND  49
            KL = KL + 1                                                 CAND  50
            SUM1 = 1.                                                   CAND  51
            DO 130 I=1,N                                                CAND  52
              IF (I.EQ.JB) GO TO 130                                    CAND  53
              MM = NWT(I)                                               CAND  54
              SUM = 0.                                                  CAND  55
              DO 120 J=1,MM                                             CAND  56
                SUM = SUM + A(J,K1,I)*A(J,K2,I)                         CAND  57
 120          CONTINUE                                                  CAND  58
              SUM1 = SUM*SUM1                                           CAND  59
 130        CONTINUE                                                    CAND  60
            SCR(KL) = SUM1                                              CAND  61
 140      CONTINUE                                                      CAND  62
 150    CONTINUE                                                        CAND  63
        CALL SMTINV(SCR, SCR, NF, R, IER)                               CAND  64
        IF (IER.EQ.0) GO TO 160                                         CAND  65
C  PRINT ERROR MESSAGE                                                  CAND  66
        WRITE (LPRINT,9995)                                             CAND  67
        STOP                                                            CAND  68
 160    IF (JB.GT.1) GO TO 230                                          CAND  69
C  COMPUTE S FOR WEIGHTS MATRIX (JB=1)                                  CAND  70
        DO 220 I1=1,N1                                                  CAND  71
          DO 170 L=1,NF                                                 CAND  72
            S(I1,L) = 0                                                 CAND  73
 170      CONTINUE                                                      CAND  74
          LC = 0                                                        CAND  75
          DO 210 J=1,N2                                                 CAND  76
            DO 200 K=1,J                                                CAND  77
              LC = LC + 1                                               CAND  78
              DO 190 L1=1,NF                                            CAND  79
                IF (J.EQ.K) GO TO 180                                   CAND  80
                S(I1,L1) = Y(I1,LC)*A(J,L1,2)*A(K,L1,3) +               CAND  81
     *              S(I1,L1)                                            CAND  82
 180            S(I1,L1) = Y(I1,LC)*A(K,L1,2)*A(J,L1,3) +               CAND  83
     *              S(I1,L1)                                            CAND  84
 190          CONTINUE                                                  CAND  85
 200        CONTINUE                                                    CAND  86
 210      CONTINUE                                                      CAND  87
 220    CONTINUE                                                        CAND  88
        GO TO 300                                                       CAND  89
C  COMPUTE S FOR STIM MATRICES (JB>1)                                   CAND  90
 230    JC = 5 - JB                                                     CAND  91
        JRB = 1                                                         CAND  92
        DO 290 J=1,N2                                                   CAND  93
          LC = JRB                                                      CAND  94
          DO 240 L=1,NF                                                 CAND  95
            S(J,L) = 0                                                  CAND  96
 240      CONTINUE                                                      CAND  97
          DO 280 K=1,N2                                                 CAND  98
            DO 260 I1=1,N1                                              CAND  99
              DO 250 L1=1,NF                                            CAND 100
                S(J,L1) = Y(I1,LC)*A(I1,L1,1)*A(K,L1,JC) +              CAND 101
     *              S(J,L1)                                             CAND 102
 250          CONTINUE                                                  CAND 103
 260        CONTINUE                                                    CAND 104
            IF (K.LT.J) GO TO 270                                       CAND 105
            LC = LC + K                                                 CAND 106
            GO TO 280                                                   CAND 107
 270        LC = LC + 1                                                 CAND 108
 280      CONTINUE                                                      CAND 109
          JRB = JRB + J                                                 CAND 110
 290    CONTINUE                                                        CAND 111
C  COMPUTE A                                                            CAND 112
 300    NN = NWT(JB)                                                    CAND 113
        DO 340 I=1,NN                                                   CAND 114
          JRI = 1                                                       CAND 115
          DO 330 J=1,NF                                                 CAND 116
            KL = JRI                                                    CAND 117
            SUM = 0.                                                    CAND 118
            DO 320 K=1,NF                                               CAND 119
              SUM = SUM + S(I,K)*R(KL)                                  CAND 120
              IF (K.LT.J) GO TO 310                                     CAND 121
              KL = KL + K                                               CAND 122
              GO TO 320                                                 CAND 123
 310          KL = KL + 1                                               CAND 124
 320        CONTINUE                                                    CAND 125
            JRI = JRI + J                                               CAND 126
C  THE FOLLOWING STATEMENT WAS INSERTED 5-20-82
            IF(IT.NE.0) A(I,J,JB)=ALPHA*SUM+A(I,J,JB)*(1.-ALPHA)
C   THE FOLLOWING STATEMENT WAS MODIFIED 5-20-82
            IF(IT.EQ.0) A(I,J,JB) = SUM                                 CAND 127
 330      CONTINUE                                                      CAND 128
 340    CONTINUE                                                        CAND 129
 350  CONTINUE                                                          CAND 130
C  COMPUTE Y HAT AND ERROR                                              CAND 131
      ERR = 0.                                                          CAND 132
      SSQYH = 0.                                                        CAND 133
      DO 410 I1=1,N1                                                    CAND 134
        LC = 0                                                          CAND 135
        DO 400 I2=1,N2                                                  CAND 136
          DO 390 I3=1,I2                                                CAND 137
            KE = 2                                                      CAND 138
            LC = LC + 1                                                 CAND 139
            SUMM(1) = 0                                                 CAND 140
            SUMM(2) = 0.                                                CAND 141
            DO 370 I=1,NF                                               CAND 142
              IF (I3.EQ.I2) GO TO 360                                   CAND 143
              SUMM(2) = SUMM(2) + A(I1,I,1)*A(I3,I,2)*A(I2,I,3)         CAND 144
 360          SUMM(1) = SUMM(1) + A(I1,I,1)*A(I2,I,2)*A(I3,I,3)         CAND 145
 370        CONTINUE                                                    CAND 146
            YS = Y(I1,LC)                                               CAND 147
            IF (I2.EQ.I3) KE = 1                                        CAND 148
            DO 380 I=1,KE                                               CAND 149
              DIF = YS - SUMM(I)                                        CAND 150
              SSQYH = SSQYH + SUMM(I)**2                                CAND 151
              ERR = ERR + DIF**2                                        CAND 152
 380        CONTINUE                                                    CAND 153
 390      CONTINUE                                                      CAND 154
 400    CONTINUE                                                        CAND 155
 410  CONTINUE                                                          CAND 156
      YY = (ERR-SSQY-SSQYH)*(-.5)                                       CAND 157
      CORYY = YY/SQRT(SSQY*SSQYH)                                       CAND 158
      RSQ = CORYY**2                                                    CAND 159
      ERR = ERR/FLOAT(N1)                                               CAND 160
      IF (ISTI.LT.2 .AND. IT.EQ.0) WRITE (LPRINT,9992)                  CAND 161
      IF (ISTI.GT.0) GO TO 420                                          CAND 162
      WRITE (LPRINT,9990) IT, CORYY, RSQ, ERR                           CAND 163
      GO TO 430                                                         CAND 164
 420  WRITE (LPRINT,9989) CORYY, RSQ, ERR                               CAND 165
      FVAF = RSQ                                                        CAND 166
      GO TO 520                                                         CAND 167
 430  CONTINUE                                                          CAND 168
      IF (IT.LE.1) GO TO 440                                            CAND 169
      ETEST = PERR - ERR                                                CAND 170
      IF (ETEST.GT.CRIT) GO TO 440                                      CAND 171
      ITCOM = 1                                                         CAND 172
      WRITE (LPRINT,9994) IT                                            CAND 173
      GO TO 450                                                         CAND 174
 440  PERR = ERR                                                        CAND 175
      IF (IT.LT.MAXIT) GO TO 460                                        CAND 176
      WRITE (LPRINT,9993)                                               CAND 177
      ITCOM = 2                                                         CAND 178
 450  FVAF = RSQ                                                        CAND 179
 460  IF ((IPUN-1).LE.0) GO TO 500                                      CAND 180
      IF (ITCOM.GT.0) GO TO 470                                         CAND 181
C  PRINT UNNORMALIZED A MATRICES EVERY 10 TH ITERATION                  CAND 182
      IF (IT.EQ.0) GO TO 510                                            CAND 183
      IF (MOD(IT,10).NE.0) GO TO 510                                    CAND 184
 470  WRITE (LPUNCH,9999) IT                                            CAND 185
      DO 490 I=1,N                                                      CAND 186
        NN = NWT(I)                                                     CAND 187
        DO 480 J=1,NF                                                   CAND 188
          WRITE (LPUNCH,9998) I, J, (A(M,J,I),M=1,NN)                   CAND 189
 480    CONTINUE                                                        CAND 190
 490  CONTINUE                                                          CAND 191
9999  FORMAT (38H UNNORMALIZED MATRICES AFTER ITERATION, I4)            CAND 192
9998  FORMAT (I3, 1X, I3, 10F7.3/(7X, 10F7.3))                          CAND 193
 500  IF (ITCOM.GT.0) GO TO 520                                         CAND 194
 510  IT = IT + 1                                                       CAND 195
      GO TO 100                                                         CAND 196
C  SET TABLE 2=TABLE 3; ITERATE AGAIN TO COMPUTE SUBJ. WEIGHTS          CAND 197
 520  IF (ISTI.NE.0) GO TO 550                                          CAND 198
      DO 540 I=1,NF                                                     CAND 199
        DO 530 J=1,N2                                                   CAND 200
          A(J,I,2) = A(J,I,3)                                           CAND 201
 530    CONTINUE                                                        CAND 202
 540  CONTINUE                                                          CAND 203
      MAXIT = 0                                                         CAND 204
      ISTI = 2                                                          CAND 205
      GO TO 90                                                          CAND 206
 550  RETURN                                                            CAND 207
9997  FORMAT (80A1)                                                     CAND 208
9996  FORMAT (1H , I3, 10F7.3/(1H 3X, 10F7.3))                          CAND 209
9995  FORMAT (21H0R CANNOT BE INVERTED)                                 CAND 210
9994  FORMAT (31H0REACHED CRITERION ON ITERATION, I3)                   CAND 211
9993  FORMAT (26HREACHED MAXIMUM ITERATIONS)                            CAND 212
9992  FORMAT (1H0, 23X, 22HHISTORY OF COMPUTATION/10H ITERATION,        CAND 213
     *    8X, 20HCORRELATIONS BETWEEN, 7X, 3HVAF, 17X, 4HLOSS/          CAND 214
     *    19X, 16HY(DATA) AND YHAT, 8X, 6H(R**2), 12X, 6H(Y-YHA,        CAND 215
     *    5HT)**2)                                                      CAND 216
9991  FORMAT (24H0INITIAL STIMULUS MATRIX)                              CAND 217
9990  FORMAT (I7, 3X, 3F20.6)                                           CAND 218
9989  FORMAT (4X, 5HFINAL, 1X, 3F20.6)                                  CAND 219
      END                                                               CAND 220
C                                                                       SSQC   1
      SUBROUTINE SSQC(Y, N1, LSD, USSQ, N2, TSSQ, IRDATA)               SSQC   2
C SSQC           FOR SINDSCAL              SEPTEMBER, 1975              SSQC   3
C WRITTEN BY S. PRUZANSKY                                               SSQC   4
C  SUBROUTINE TO COMPUTE SSQ (WITH OPTION TO NORMALIZE)                 SSQC   5
C  FOR CORR. AND COVAR. DATA                                            SSQC   6
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             SSQC   7
      DIMENSION Y(N1,LSD), USSQ(N2)                                     SSQC   8
      TSSQ = 0                                                          SSQC   9
      DO 70 I1=1,N1                                                     SSQC  10
        SSQD = 0                                                        SSQC  11
        SSQ = 0                                                         SSQC  12
        K = 0                                                           SSQC  13
        DO 30 J=1,N2                                                    SSQC  14
          DO 20 L=1,J                                                   SSQC  15
            K = K + 1                                                   SSQC  16
            IF (J.EQ.L) GO TO 10                                        SSQC  17
            SSQ = Y(I1,K)**2 + SSQ                                      SSQC  18
            GO TO 20                                                    SSQC  19
  10        SSQD = Y(I1,K)**2 + SSQD                                    SSQC  20
  20      CONTINUE                                                      SSQC  21
  30    CONTINUE                                                        SSQC  22
        USSQ(I1) = SSQ*2. + SSQD                                        SSQC  23
        IF (IRDATA.GT.0) GO TO 60                                       SSQC  24
C NORMALIZE DATA SO SSQ=1                                               SSQC  25
        TSSQ = N1                                                       SSQC  26
        Q = 1./SQRT(USSQ(I1))                                           SSQC  27
        K1 = 0                                                          SSQC  28
        DO 50 I=1,N2                                                    SSQC  29
          DO 40 J=1,I                                                   SSQC  30
            K1 = K1 + 1                                                 SSQC  31
            Y(I1,K1) = Y(I1,K1)*Q                                       SSQC  32
  40      CONTINUE                                                      SSQC  33
  50    CONTINUE                                                        SSQC  34
        GO TO 70                                                        SSQC  35
  60    TSSQ = TSSQ + USSQ(I1)                                          SSQC  36
  70  CONTINUE                                                          SSQC  37
      IF (IRDATA.LT.0) RETURN                                           SSQC  38
C PRINT SSQ FOR UNNORMALIZED DATA                                       SSQC  39
      WRITE (LPRINT,9999)                                               SSQC  40
      DO 80 I=1,N1,10                                                   SSQC  41
        IE = I + 9                                                      SSQC  42
        IF (IE.GT.N1) IE = N1                                           SSQC  43
        WRITE (LPRINT,9998) I, (USSQ(L),L=I,IE)                         SSQC  44
  80  CONTINUE                                                          SSQC  45
9999  FORMAT (46H0SUM OF SQUARES OF UNNORMALIZED CORR. OR COVAR,        SSQC  46
     *    2H. , 8HMATRICES)                                             SSQC  47
9998  FORMAT (1H I3, 10F7.2)                                            SSQC  48
      RETURN                                                            SSQC  49
      END                                                               SSQC  50
C                                                                       NORM   1
      SUBROUTINE NORMA(A, AA, MNWT, K, N1, N2, NA, FVAF, S,             NORM   2
     *    IPUN)                                                         NORM   3
C NORMA          FOR SINDSCAL              SEPTEMBER, 1975              NORM   4
C WRITTEN BY J.J. CHANG                                                 NORM   5
C MODIFIED FOR SINDSCAL BY S. PRUZANSKY                                 NORM   6
C     PROGRAM TO NORMALIZE A MATRICES FROM THE CANNONICAL DECOMPOSITION NORM   7
C     OF 3 WAY TABLES - TABLE 2 MUST EQ. TABLE 3                        NORM   8
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             NORM   9
      DIMENSION NWT(3), A(MNWT,K,3), S(K,K)                             NORM  10
      DIMENSION AA(MNWT,K), DIAG(10), LPAS(10), V(10)                   NORM  11
      NWT(1) = N1                                                       NORM  12
      NWT(2) = N2                                                       NORM  13
C  NORMALIZE STIM VECTORS SO SSQ OF EACH DIM EQ 1                       NORM  14
      XN1 = NWT(1)                                                      NORM  15
      I = 2                                                             NORM  16
      DO 30 J=1,K                                                       NORM  17
        SUM = 0.                                                        NORM  18
        DO 10 M=1,N2                                                    NORM  19
          SUM = SUM + A(M,J,I)**2                                       NORM  20
  10    CONTINUE                                                        NORM  21
        V(J) = SUM                                                      NORM  22
        SUM = SQRT(SUM)                                                 NORM  23
        DO 20 M=1,N2                                                    NORM  24
          A(M,J,I) = A(M,J,I)/SUM                                       NORM  25
  20    CONTINUE                                                        NORM  26
  30  CONTINUE                                                          NORM  27
C  MULT. ELEMENTS OF WEIGHTS VECTOR BY SSQ OF CORRES. DIM               NORM  28
      DO 50 J=1,K                                                       NORM  29
        DO 40 M=1,N1                                                    NORM  30
          A(M,J,1) = A(M,J,1)*V(J)                                      NORM  31
  40    CONTINUE                                                        NORM  32
  50  CONTINUE                                                          NORM  33
C  COMPUTE SUM SQUARES OF WEIGHTS FOR EACH DIM.                         NORM  34
      DO 70 I=1,K                                                       NORM  35
        LPAS(I) = I                                                     NORM  36
        DIAG(I) = 0.                                                    NORM  37
        DO 60 J=1,N1                                                    NORM  38
          DIAG(I) = DIAG(I) + A(J,I,1)**2                               NORM  39
  60    CONTINUE                                                        NORM  40
  70  CONTINUE                                                          NORM  41
      CALL SORT(DIAG, K, LPAS, -1)                                      NORM  42
C  PERMUTE DIMENSIONS ACCORDING TO SUM SQUARES IN DESCENDING ORDER      NORM  43
      DO 120 I=1,NA                                                     NORM  44
        NN = NWT(I)                                                     NORM  45
        DO 90 J=1,K                                                     NORM  46
          DO 80 M=1,NN                                                  NORM  47
            AA(M,J) = A(M,J,I)                                          NORM  48
  80      CONTINUE                                                      NORM  49
  90    CONTINUE                                                        NORM  50
        DO 110 J=1,K                                                    NORM  51
          JJ = LPAS(J)                                                  NORM  52
          DO 100 M=1,NN                                                 NORM  53
            A(M,J,I) = AA(M,JJ)                                         NORM  54
 100      CONTINUE                                                      NORM  55
 110    CONTINUE                                                        NORM  56
 120  CONTINUE                                                          NORM  57
C  PRINT AND PUNCH NORMALIZED SOLUTION                                  NORM  58
      WRITE (LPRINT,9999)                                               NORM  59
      IF (IPUN.GE.0) WRITE (LPUNCH,9994)                                NORM  60
      DO 190 I=1,NA                                                     NORM  61
        NN = NWT(I)                                                     NORM  62
        IF (IPUN.LT.0) GO TO 140                                        NORM  63
        DO 130 J=1,K                                                    NORM  64
          WRITE (LPUNCH,9993) J, I, (A(M,J,I),M=1,NN)                   NORM  65
 130    CONTINUE                                                        NORM  66
 140    GO TO (150, 160), I                                             NORM  67
 150    WRITE (LPRINT,9996)                                             NORM  68
        GO TO 170                                                       NORM  69
 160    WRITE (LPRINT,9992)                                             NORM  70
        GO TO 170                                                       NORM  71
 170    DO 180 J=1,K                                                    NORM  72
          WRITE (LPRINT,9997) J, (A(M,J,I),M=1,NN)                      NORM  73
 180    CONTINUE                                                        NORM  74
 190  CONTINUE                                                          NORM  75
C  COMPUTE SUMS OF PRODUCTS AND SUMS OF SQUARES OF A MATRIX             NORM  76
      DO 290 I=1,NA                                                     NORM  77
        IF (I.GT.1) WRITE (LPRINT,9998)                                 NORM  78
        NN = NWT(I)                                                     NORM  79
        DO 220 L=1,K                                                    NORM  80
          DO 210 J=1,K                                                  NORM  81
            S(L,J) = 0.                                                 NORM  82
            DO 200 M=1,NN                                               NORM  83
              S(L,J) = S(L,J) + A(M,L,I)*A(M,J,I)                       NORM  84
 200        CONTINUE                                                    NORM  85
 210      CONTINUE                                                      NORM  86
 220    CONTINUE                                                        NORM  87
        SUM = 0.                                                        NORM  88
        DO 240 L=1,K                                                    NORM  89
          IF (I.EQ.1) GO TO 230                                         NORM  90
          WRITE (LPRINT,9997) L, (S(L,J),J=1,L)                         NORM  91
          GO TO 240                                                     NORM  92
 230      V(L) = S(L,L)                                                 NORM  93
          SUM = SUM + V(L)                                              NORM  94
 240    CONTINUE                                                        NORM  95
        IF (I.GT.1) GO TO 290                                           NORM  96
        DO 260 II=1,K                                                   NORM  97
          DO 250 JJ=1,II                                                NORM  98
            S(II,JJ) = S(II,JJ)/SQRT(V(II)*V(JJ))                       NORM  99
 250      CONTINUE                                                      NORM 100
 260    CONTINUE                                                        NORM 101
        WRITE (LPRINT,9990)                                             NORM 102
        DO 270 L=1,K                                                    NORM 103
          WRITE (LPRINT,9997) L, (S(L,J),J=1,L)                         NORM 104
 270    CONTINUE                                                        NORM 105
        DO 280 II=1,K                                                   NORM 106
          V(II) = (V(II)*FVAF)/SUM                                      NORM 107
 280    CONTINUE                                                        NORM 108
 290  CONTINUE                                                          NORM 109
      DO 300 II=1,K                                                     NORM 110
        LPAS(II) = II                                                   NORM 111
 300  CONTINUE                                                          NORM 112
      WRITE (LPRINT,9991) (LPAS(II),II=1,K)                             NORM 113
      WRITE (LPRINT,9995) (V(II),II=1,K)                                NORM 114
9999  FORMAT (1H0, 23X, 19HNORMALIZED SOLUTION)                         NORM 115
9998  FORMAT (26H0SUM OF PRODUCTS (STIMULI))                            NORM 116
9997  FORMAT (I3, 3X, 10F7.3/(6X, 10F7.3))                              NORM 117
9996  FORMAT (23H0SUBJECTS WEIGHT MATRIX)                               NORM 118
9995  FORMAT (6X, 10F7.3)                                               NORM 119
9994  FORMAT (35H(2X,5E13.5)     NORMALIZED SOLUTION)                   NORM 120
9993  FORMAT (2I1, 5E13.5/(2X, 5E13.5))                                 NORM 121
9992  FORMAT (16H0STIMULUS MATRIX)                                      NORM 122
9991  FORMAT (46H0APPROXIMATE PROPORTION OF TOTAL VARIANCE ACCO,        NORM 123
     *    9HUNTED FOR, 18H BY EACH DIMENSION/6X, 10(I4, 3X))            NORM 124
9990  FORMAT (38H0NORMALIZED SUM OF PRODUCTS (SUBJECTS))                NORM 125
      RETURN                                                            NORM 126
      END                                                               NORM 127
C                                                                       SUBR   1
      SUBROUTINE SUBR(A, MNWT, NF, N, Y, N1, N2, LSD)                   SUBR   2
C SUBR           FOR SINDSCAL              SEPTEMBER, 1975              SUBR   3
C WRITTEN BY J.J. CHANG                                                 SUBR   4
C MODIFIED FOR SINDSCAL BY S. PRUZANSKY                                 SUBR   5
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             SUBR   6
      DIMENSION A(MNWT,NF,N), Y(N1,LSD)                                 SUBR   7
      XN = N2*N2                                                        SUBR   8
      WRITE (LPRINT,9999)                                               SUBR   9
      DO 50 I1=1,N1                                                     SUBR  10
        SUMX = 0.                                                       SUBR  11
        SUMY = 0.                                                       SUBR  12
        SSQX = 0.                                                       SUBR  13
        SSQY = 0.                                                       SUBR  14
        SXY = 0.                                                        SUBR  15
        LC = 0                                                          SUBR  16
        DO 40 I2=1,N2                                                   SUBR  17
          DO 30 I3=1,I2                                                 SUBR  18
            KE = 2                                                      SUBR  19
            LC = LC + 1                                                 SUBR  20
            SUM = 0                                                     SUBR  21
            YY = Y(I1,LC)                                               SUBR  22
            DO 10 I=1,NF                                                SUBR  23
              SUM = SUM + A(I1,I,1)*A(I2,I,2)*A(I3,I,2)                 SUBR  24
  10        CONTINUE                                                    SUBR  25
            IF (I2.EQ.I3) KE = 1                                        SUBR  26
            DO 20 I=1,KE                                                SUBR  27
              SUMX = SUMX + SUM                                         SUBR  28
              SUMY = SUMY + YY                                          SUBR  29
              SSQX = SSQX + SUM**2                                      SUBR  30
              SSQY = SSQY + YY**2                                       SUBR  31
              SXY = SXY + SUM*YY                                        SUBR  32
  20        CONTINUE                                                    SUBR  33
  30      CONTINUE                                                      SUBR  34
  40    CONTINUE                                                        SUBR  35
        R = (XN*SXY-SUMX*SUMY)/SQRT((XN*SSQX-SUMX**2)*(XN*              SUBR  36
     *      SSQY-SUMY**2))                                              SUBR  37
        WRITE (LPRINT,9998) I1, R                                       SUBR  38
  50  CONTINUE                                                          SUBR  39
9999  FORMAT (46H0CORRELATION BETWEEN COMPUTED SCORES AND SCALA,        SUBR  40
     *    8HR PROD. , 12HFOR SUBJECTS/1H0)                              SUBR  41
9998  FORMAT (I4, 3X, F10.6)                                            SUBR  42
      RETURN                                                            SUBR  43
      END                                                               SUBR  44
C                                                                       DIMA   1
      SUBROUTINE DIMA(A, MNWT, NQ, NFO, NF)                             DIMA   2
C DIMA           FOR SINDSCAL              SEPTEMBER, 1975              DIMA   3
C WRITTEN BY S. PRUZANSKY                                               DIMA   4
      DIMENSION A(MNWT,NQ)                                              DIMA   5
      NN = NFO                                                          DIMA   6
      L = NF                                                            DIMA   7
      DO 20 I=1,NF                                                      DIMA   8
        L = L + 1                                                       DIMA   9
        NN = NN + 1                                                     DIMA  10
        DO 10 J=1,MNWT                                                  DIMA  11
          A(J,L) = A(J,NN)                                              DIMA  12
  10    CONTINUE                                                        DIMA  13
  20  CONTINUE                                                          DIMA  14
      RETURN                                                            DIMA  15
      END                                                               DIMA  16
C                                                                       SMTI   1
      SUBROUTINE SMTINV(A, UL, N, AINV, IER)                            SMTI   2
C-SMTINV-------ADAPTED FROM IMSL LIBRARY 2---------------------------   SMTI   3
C  SMTINV MERGES THE TWO IMSL ROUTINES LINV1P AND LUDECP.               SMTI   4
C    THESE IMSL ROUTINES ARE PROPRIETARY SOFTWARE, OWNED BY IMSL.       SMTI   5
C    THEY MAY BE USED ONLY IN THE CODE IN WHICH THEY ARE EMBEDDED.      SMTI   6
C    NO OTHER USE OF THESE ROUTINES IS PERMITTED.                       SMTI   7
C                                                                       SMTI   8
C   FUNCTION            - INVERSION OF A POSITIVE DEFINITE SYMMETRIC    SMTI   9
C                           MATRIX - SYMMETRIC STORAGE MODE - SPACE     SMTI  10
C                           ECONOMIZER SOLUTION                         SMTI  11
C   USAGE               - CALL SMTINV(A,UL,N,AINV,IER)                  SMTI  12
C   PARAMETERS   A      - INPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING    SMTI  13
C                           THE N BY N POSITIVE DEFINITE SYMMETRIC      SMTI  14
C                           MATRIX TO BE INVERTED. A IS STORED IN       SMTI  15
C                           SYMMETRIC STORAGE MODE.                     SMTI  16
C               UL         ON OUTPUT, A IS REPLACED BY THE LOWER        SMTI  17
C                           TRIANGULAR MATRIX L WHERE A = L*L-TRANSPOSE.SMTI  18
C                           L IS STORED IN SYMMETRIC STORAGE MODE WITH  SMTI  19
C                           THE DIAGONAL ELEMENTS OF L STORED IN        SMTI  20
C                           RECIPROCAL FORM.                            SMTI  21
C                N      - ORDER OF A.(INPUT)                            SMTI  22
C                AINV   - OUTPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING   SMTI  23
C                           THE N BY N INVERSE OF A. AINV IS STORED IN  SMTI  24
C                           SYMMETRIC STORAGE MODE.                     SMTI  25
C                IER    - ERROR PARAMETER.                              SMTI  26
C                         TERMINAL ERROR = 128+N.                       SMTI  27
C   PRECISION           - SINGLE                                        SMTI  28
C   REQ'D IMSL ROUTINES - LUELMP                                        SMTI  29
C-----------------------------------------------------------------------SMTI  30
      DIMENSION A(1), AINV(1), UL(1)                                    SMTI  31
      DOUBLE PRECISION X                                                SMTI  32
      DATA ZERO, ONE, FOUR, SIXTN, SIXTH /0.0,1.,4.,16.,.0625/          SMTI  33
      IER = 0                                                           SMTI  34
      L = 1                                                             SMTI  35
      N1 = N - 1                                                        SMTI  36
      LN = N                                                            SMTI  37
C                                  DECOMPOSE A                          SMTI  38
      D1 = ONE                                                          SMTI  39
      D2 = ZERO                                                         SMTI  40
      IP = 1                                                            SMTI  41
      IER = 0                                                           SMTI  42
      DO 90 I=1,N                                                       SMTI  43
        IQ = IP                                                         SMTI  44
        IR = 1                                                          SMTI  45
        DO 80 J=1,I                                                     SMTI  46
          X = A(IP)                                                     SMTI  47
          IF (J.EQ.1) GO TO 20                                          SMTI  48
          DO 10 K=IQ,IP1                                                SMTI  49
            X = X - DBLE(UL(K))*DBLE(UL(IR))                            SMTI  50
            IR = IR + 1                                                 SMTI  51
  10      CONTINUE                                                      SMTI  52
  20      IF (I.NE.J) GO TO 60                                          SMTI  53
          D1 = D1*X                                                     SMTI  54
          IF (X.LE.0.D0) GO TO 100                                      SMTI  55
  30      IF (ABS(D1).LT.ONE) GO TO 40                                  SMTI  56
          D1 = D1*SIXTH                                                 SMTI  57
          D2 = D2 + FOUR                                                SMTI  58
          GO TO 30                                                      SMTI  59
  40      IF (ABS(D1).GE.SIXTH) GO TO 50                                SMTI  60
          D1 = D1*SIXTN                                                 SMTI  61
          D2 = D2 - FOUR                                                SMTI  62
          GO TO 40                                                      SMTI  63
  50      UL(IP) = 1.D0/DSQRT(X)                                        SMTI  64
          GO TO 70                                                      SMTI  65
  60      UL(IP) = X*UL(IR)                                             SMTI  66
  70      IP1 = IP                                                      SMTI  67
          IP = IP + 1                                                   SMTI  68
          IR = IR + 1                                                   SMTI  69
  80    CONTINUE                                                        SMTI  70
  90  CONTINUE                                                          SMTI  71
      GO TO 110                                                         SMTI  72
 100  IER = 129                                                         SMTI  73
 110  CONTINUE                                                          SMTI  74
      IF (IER.NE.0) GO TO 140                                           SMTI  75
      DO 130 I=1,N                                                      SMTI  76
        DO 120 J=L,LN                                                   SMTI  77
          AINV(J) = ZERO                                                SMTI  78
 120    CONTINUE                                                        SMTI  79
        LI = L + I - 1                                                  SMTI  80
        AINV(LI) = ONE                                                  SMTI  81
C                                  OBTAIN THE SOLUTION                  SMTI  82
        CALL LUELMP(UL, AINV(L), N, AINV(L))                            SMTI  83
        L = L + I                                                       SMTI  84
        LN = L + N1                                                     SMTI  85
 130  CONTINUE                                                          SMTI  86
 140  RETURN                                                            SMTI  87
      END                                                               SMTI  88
C                                                                       LUEL   1
      SUBROUTINE LUELMP(A, B, N, X)                                     LUEL   2
C-LUELMP--------FROM IMSL LIBRARY 2----------------------------------   LUEL   3
C  THE IMSL ROUTINE LUELMP IS PROPRIETARY SOFTWARE, OWNED BY IMSL.      LUEL   4
C    IT MAY BE USED ONLY IN THE CODE IN WHICH IT IS EMBEDDED.           LUEL   5
C    NO OTHER USE OF THIS ROUTINE IS PERMITTED.                         LUEL   6
C                                                                       LUEL   7
C   FUNCTION            - ELIMINATION PART OF THE SOLUTION OF AX=B -    LUEL   8
C                           SYMMETRIC STORAGE MODE                      LUEL   9
C   USAGE               - CALL LUELMP (A,B,N,X)                         LUEL  10
C   PARAMETERS   A      - INPUT VECTOR OF LENGTH N(N+1)/2 CONTAINING    LUEL  11
C                           THE N BY N MATRIX L WHERE A = L*L-TRANSPOSE.LUEL  12
C                           L IS A LOWER TRIANGULAR MATRIX STORED IN    LUEL  13
C                           SYMMETRIC STORAGE MODE. THE MAIN DIAGONAL   LUEL  14
C                           ELEMENTS OF L ARE STORED IN RECIPROCAL      LUEL  15
C                           FORM. MATRIX L MAY BE OBTAINED FROM IMSL    LUEL  16
C                           ROUTINE LUDECP.                             LUEL  17
C                B      - VECTOR OF LENGTH N CONTAINING THE RIGHT HAND  LUEL  18
C                           SIDE OF THE EQUATION AX = B. (INPUT)        LUEL  19
C                N      - ORDER OF A AND THE LENGTH OF B AND X. (INPUT) LUEL  20
C                X      - VECTOR OF LENGTH N CONTAINING THE SOLUTION TO LUEL  21
C                           THE EQUATION AX = B. (OUTPUT)               LUEL  22
C   PRECISION           - SINGLE                                        LUEL  23
C   LANGUAGE            - FORTRAN                                       LUEL  24
C-----------------------------------------------------------------------LUEL  25
C   LATEST REVISION     - FEBRUARY 12,1973                              LUEL  26
C                                                                       LUEL  27
      DIMENSION A(1), B(1), X(1)                                        LUEL  28
      DOUBLE PRECISION T                                                LUEL  29
C                                  SOLUTION OF LY = B                   LUEL  30
      IP = 1                                                            LUEL  31
      DO 30 I=1,N                                                       LUEL  32
        IQ = I - 1                                                      LUEL  33
        T = B(I)                                                        LUEL  34
        IF (IQ.LT.1) GO TO 20                                           LUEL  35
        DO 10 K=1,IQ                                                    LUEL  36
          T = T - DBLE(A(IP))*DBLE(X(K))                                LUEL  37
          IP = IP + 1                                                   LUEL  38
  10    CONTINUE                                                        LUEL  39
  20    X(I) = T*A(IP)                                                  LUEL  40
        IP = IP + 1                                                     LUEL  41
  30  CONTINUE                                                          LUEL  42
C                                  SOLUTION OF UX = Y                   LUEL  43
      N1 = N + 1                                                        LUEL  44
      DO 60 I=1,N                                                       LUEL  45
        II = N1 - I                                                     LUEL  46
        IP = IP - 1                                                     LUEL  47
        IS = IP                                                         LUEL  48
        IQ = II + 1                                                     LUEL  49
        T = X(II)                                                       LUEL  50
        IF (N.LT.IQ) GO TO 50                                           LUEL  51
        KK = N                                                          LUEL  52
        DO 40 K=IQ,N                                                    LUEL  53
          T = T - DBLE(A(IS))*DBLE(X(KK))                               LUEL  54
          KK = KK - 1                                                   LUEL  55
          IS = IS - KK                                                  LUEL  56
  40    CONTINUE                                                        LUEL  57
  50    X(II) = T*A(IS)                                                 LUEL  58
  60  CONTINUE                                                          LUEL  59
      RETURN                                                            LUEL  60
      END                                                               LUEL  61
C                                                                       SORT   1
CSORT                                                                   SORT   2
      SUBROUTINE SORT(A, N, LPAS, SWITCH)                               SORT   3
C SORT             FOR KYST                  JANUARY,1973               SORT   4
C WRITTEN BY J.KRUSKAL  MODIFIED FOR SINDSCAL 8-18-77 BY SP             SORT   5
C                                                                       SORT   6
C     THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY,    SORT   7
C     ARRAY LPAS, IN ORDER CORRESPONDING TO 'A'.                        SORT   8
C     N = NUMBER OF ITEMS IN 'A' (AND LPAS         , IF USED)           SORT   9
C     IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER,         SORT  10
C                 IF ZERO OR NEGATIVE, IN DESCENDING ORDER.             SORT  11
C     ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL            SORT  12
C                                                                       SORT  13
      DIMENSION A(1), LPAS(1)                                           SORT  14
      INTEGER SWITCH                                                    SORT  15
      IF (N.LE.1) GO TO 90                                              SORT  16
      M = 1                                                             SORT  17
  10  M = M + M                                                         SORT  18
      IF (M.LE.N) GO TO 10                                              SORT  19
      M = M - 1                                                         SORT  20
  20  M = M/2                                                           SORT  21
      IF (M.EQ.0) GO TO 90                                              SORT  22
      KK = N - M                                                        SORT  23
      J = 1                                                             SORT  24
  30  I = J                                                             SORT  25
  40  IM = I + M                                                        SORT  26
      IF (SWITCH) 60, 60, 50                                            SORT  27
  50  IF (A(I).GT.A(IM)) GO TO 80                                       SORT  28
      GO TO 70                                                          SORT  29
  60  IF (A(I).LT.A(IM)) GO TO 80                                       SORT  30
  70  J = J + 1                                                         SORT  31
      IF (J.GT.KK) GO TO 20                                             SORT  32
      GO TO 30                                                          SORT  33
  80  TEMP = A(I)                                                       SORT  34
      A(I) = A(IM)                                                      SORT  35
      A(IM) = TEMP                                                      SORT  36
      ITEMP = LPAS(I)                                                   SORT  37
      LPAS(I) = LPAS(IM)                                                SORT  38
      LPAS(IM) = ITEMP                                                  SORT  39
      I = I - M                                                         SORT  40
      IF (I.LT.1) GO TO 70                                              SORT  41
      GO TO 40                                                          SORT  42
  90  RETURN                                                            SORT  43
      END                                                               SORT  44
C                                                                       RCON   1
      SUBROUTINE RCON(A, MNWT, NF, N2, IRN)                             RCON   2
C RCON             FOR SINDSCAL              SEPTEMBER, 1975            RCON   3
C WRITTEN BY S. PRUZANSKY                                               RCON   4
C  SUBROUTINE TO GENERATE RANDOM STARTING CONFIGURATION                 RCON   5
C       USES PORTABLE RANDOM NUMBER GENERATOR                           RCON   6
      DIMENSION A(MNWT,NF,3)                                            RCON   7
      IF (IRN.LT.100) IRN = IRN*11                                      RCON   8
      DO 10 I=1,IRN                                                     RCON   9
        TEMP = RUNIFV(DUM)                                              RCON  10
  10  CONTINUE                                                          RCON  11
      DO 40 L=2,3                                                       RCON  12
        DO 30 K=1,NF                                                    RCON  13
          DO 20 J=1,N2                                                  RCON  14
            TEMP = RUNIFV(DUM)                                          RCON  15
            A(J,K,L) = TEMP - .5                                        RCON  16
  20      CONTINUE                                                      RCON  17
  30    CONTINUE                                                        RCON  18
  40  CONTINUE                                                          RCON  19
      RETURN                                                            RCON  20
      END                                                               RCON  21
C                                                                       RUNI   1
      FUNCTION RUNIFV(A)                                                RUNI   2
C RUNIFV           FOR KYST                  JANUARY,1973               RUNI   3
C WRITTEN BY J.KRUSKAL                                                  RUNI   4
C      THIS IS A SIMPLE ROUTINE FOR GENERATING RANDOM NUMBERS.          RUNI   5
C      THEY ARE UNIFORMLY DISTRIBUTED BETWEEN 0 AND 1.                  RUNI   6
C      IT IS NOT FAST, NOR DOES IT PRODUCE VERY HIGH QUALITY NUMBERS.   RUNI   7
C      ITS MAIN VIRTUE IS THAT IT IS IN FORTRAN AND WILL WORK ON        RUNI   8
C      ALMOST ANY MACHINE ON WHICH FORTRAN IS AVAILABLE.                RUNI   9
C                                                                       RUNI  10
      DATA MODULO, FLMOD, K /8192,8192.0,1/                             RUNI  11
C      MODULO   2**13                                                   RUNI  12
C                                                                       RUNI  13
      DO 10 I=1,15                                                      RUNI  14
        K = MOD(5*K,MODULO)                                             RUNI  15
  10  CONTINUE                                                          RUNI  16
      Z = FLOAT(K)/FLMOD                                                RUNI  17
      RUNIFV = Z                                                        RUNI  18
      RETURN                                                            RUNI  19
      END                                                               RUNI  20
C                                                                       PPLO   1
      SUBROUTINE PPLOT(X, MNWT, NF, TIT, N)                             PPLO   2
C PPLOT          FOR SINDSCAL              SEPTEMBER, 1975              PPLO   3
C WRITTEN BY S. PRUZANSKY                                               PPLO   4
C  SUBROUTINE TO PLOT CONFIGURATION ON PRINTER                          PPLO   5
C     REQUIRES 120 CHARACTERS PER LINE                                  PPLO   6
C        WRITTEN 2-6-75 BY SP                                           PPLO   7
      DIMENSION X(MNWT,NF,3), N(1)                                      PPLO   8
      INTEGER TIT(1)                                                    PPLO   9
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             PPLO  10
      DO 30 L1=1,2                                                      PPLO  11
        DO 20 J=2,NF                                                    PPLO  12
          J1 = J - 1                                                    PPLO  13
          DO 10 K=1,J1                                                  PPLO  14
C  DETERMINE MAX. ABSOLUTE VALUE FOR PLANE  TO BE PLOTTED               PPLO  15
            CALL MAXI(X(1,K,L1), X(1,J,L1), XMAX, N(L1))                PPLO  16
            XMIN = -XMAX                                                PPLO  17
C  PRINT TITLE INFORMATION AT TOP OF PAGE                               PPLO  18
            IF (L1.EQ.1) WRITE (LPRINT,9999) J, K,                      PPLO  19
     *          (TIT(II),II=1,80)                                       PPLO  20
            IF (L1.EQ.2) WRITE (LPRINT,9998) J, K,                      PPLO  21
     *          (TIT(II),II=1,80)                                       PPLO  22
C  PLOT A PLANE AT A TIME                                               PPLO  23
            CALL PLOT(X(1,K,L1), X(1,J,L1), XMAX, XMAX, XMIN,           PPLO  24
     *          XMIN, N(L1), 1, 2, 66, 1)                               PPLO  25
  10      CONTINUE                                                      PPLO  26
  20    CONTINUE                                                        PPLO  27
  30  CONTINUE                                                          PPLO  28
9999  FORMAT (26H1WEIGHTS SPACE   DIMENSION, I2, 11H (Y-AXIS) V,        PPLO  29
     *    12HS. DIMENSION, I3, 9H (X-AXIS), /, 30X, 80A1)               PPLO  30
9998  FORMAT (27H1STIMULUS SPACE   DIMENSION, I2, 10H (Y-AXIS) ,        PPLO  31
     *    13HVS. DIMENSION, I3, 9H (X-AXIS), /, 30X, 80A1)              PPLO  32
      RETURN                                                            PPLO  33
      END                                                               PPLO  34
C                                                                       MAXI   1
      SUBROUTINE MAXI(X, Y, Z, N)                                       MAXI   2
C MAXI           FOR SINDSCAL              SEPTEMBER, 1975              MAXI   3
C WRITTEN BY J.J. CHANG                                                 MAXI   4
C     MAXI FINDS THE MAGNITUDE OF ARRAYS X AND Y EACH OF LENGTH N.      MAXI   5
C     ON RETURN ABSOLUTE MAXIMUM IS STORED IN Z.                        MAXI   6
      DIMENSION X(1), Y(1)                                              MAXI   7
      Z = ABS(X(1))                                                     MAXI   8
      DO 10 I=1,N                                                       MAXI   9
        AXO = ABS(X(I))                                                 MAXI  10
        IF (AXO.GT.Z) Z = AXO                                           MAXI  11
        AYO = ABS(Y(I))                                                 MAXI  12
        IF (AYO.GT.Z) Z = AYO                                           MAXI  13
  10  CONTINUE                                                          MAXI  14
      RETURN                                                            MAXI  15
      END                                                               MAXI  16
C                                                                       XA00   1
      BLOCK DATA                                                        XA00   2
C BLDA             FOR KYST                  JANUARY,1973               XA00   3
C WRITTEN BY J.KRUSKAL                                                  XA00   4
C MODIFIED FOR SINDSCAL BY S PRUZANSKY JUNE,1975                        XA00   5
C                                                                       XA00   6
      REAL PTID(100), ITEM(101)                                         XA00   7
C                                                                       XA00   8
      COMMON /ACCUR/ XMAXN                                              XA00   9
      COMMON /PLTCHR/ PTID, ITEM                                        XA00  10
C                                                                       XA00  11
C  XMAXN = MINIMUM OF LARGEST POSITIVE MACHINE NUMBER AND MINUS
C   LARGEST NEGATIVE MACHINE NUMBER
C     DATA XMAXN /1.0E38/                                               XA00  12
      DATA PTID /1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,           XA00  13
     *    1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,          XA00  14
     *    1HY,1HZ,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,          XA00  15
     *    1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,          XA00  16
     *    1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2,          XA00  17
     *    1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HA,1HB,1HC,1HD,1HE,1HF,          XA00  18
     *    1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO,1HP,1HQ,1HR,1HS,          XA00  19
     *    1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H1,1H2,1H3,1H4/                  XA00  20
      DATA ITEM /101*1H /                                               XA00  21
      END                                                               XA00  22
C                                                                       PLOT   1
      SUBROUTINE PLOT(X, Y, XA, YA, XI, YI, NPOI, NY, ID, LONG,         PLOT   2
     *    MX)                                                           PLOT   3
C PLOT             FOR KYST                  JANUARY,1973               PLOT   4
C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965                           PLOT   5
C     UPDATED OCTOBER 11,1971  JAY R LEVINSOHN                          PLOT   6
C ALTERED FOR KYST BY J.KRUSKAL AND J.SEERY  JANUARY,1973               PLOT   7
C                                                                       PLOT   8
C ROUTINE TO GENERATE A ONE PAGE PLOT OF A MATRIX Y VS SOME VECTOR X    PLOT   9
C                                                                       PLOT  10
C Y MAY BE A MATRIX AND IT IS ASSUMED THAT X IS A VECTOR                PLOT  11
C                                                                       PLOT  12
C THE PARAMETERS XA,YA,XI,YI INDICATE THE UPPER                         PLOT  13
C AND LOWER BOUNDS FOR EACH AXIS.                                       PLOT  14
C THE USER MAY PASS THE MIN AND MAX FOR X AND Y OR HE MAY LET THE       PLOT  15
C PROGRAM FIND THEM.  IF XA=XI THEN THE PROGRAM WILL FIND THE MAX AND MIPLOT  16
C FOR THE X VECTOR,  IF YA=YI THE PROGRAM WILL FIND THE MAX AND MIN     PLOT  17
C FOR THE Y MATRIX.                                                     PLOT  18
C                                                                       PLOT  19
C IF  -ID-  IS POSITIVE, AXES WILL BE INCLUDED ON THE PLOT,             PLOT  20
C IF  -ID-  IS NEGATIVE, NO AXES WILL APPEAR.                           PLOT  21
C IF ID IS PLUS OR MINUS 2, THE POINTS WILL BE IDENTIFIED BY ROWS       PLOT  22
C IF ID IS NOT EQUAL TO 2 THEN POINTS WILL BE LABELED BY COLS.          PLOT  23
C                                                                       PLOT  24
C IF ID IS PLUS OR MINUS THREE THE PLOT IS ASSUMED TO BE 1 DIMENSIONAL. PLOT  25
C ONE DIMENSIONAL PLOTS ARE PLOTTED WITH NO AXES, AND ALL Y VALUES      PLOT  26
C FIXED AT MAXY-MINY/2                                                  PLOT  27
C POINTS IN THE 1 DIM. PLOTS ARE IDENDTIFIED BY ROW.                    PLOT  28
C                                                                       PLOT  29
C DEFINITIONS                                                           PLOT  30
C  PARAMETERS                                                           PLOT  31
C    Y--A MATRIX (LONG,NX) IT CONTAINS VALUES TO BE PLOTTED AGAINST X   PLOT  32
C    X--A VECTOR (LONG); IT CONTAINS VALUES TO BE PLOTTED AGAINST Y     PLOT  33
C    XA-- X AXIS UPPER BOUND   (MAY BE PASSED AS ZERO)                  PLOT  34
C    YA-- Y AXIS  "     "     "      "     "                            PLOT  35
C     XI-- X AXIS LOWER BOUND  (MAY BE PASSED AS ZERO)                  PLOT  36
C    YI-- Y AXIS LOWER BOUND(MAY BE PASSED AS ZERO )                    PLOT  37
C    NPOI-- NO. OF ELEMENTS IN X, NO. OF ROWS IN Y                      PLOT  38
C    NY--NUMBER OF COLUMNS IN Y MATRIX                                  PLOT  39
C    ID-- PRINT CONTROL PARAMETER, CONTROLS LABELING, AND AXIS GEN.     PLOT  40
C     LONG-- MAXIMUM  NUMBER OF ELEMENTS IN X, MAX. NO. ROWS IN Y       PLOT  41
C  MX  -- WORK WITH ELEMENTS MX THROUGH MX+NPOI-1                       PLOT  42
C  VARIABLES                                                            PLOT  43
C    DELX-- X AXIS INCREMENT                                            PLOT  44
C    DELY-- Y AXIS INCREMENT                                            PLOT  45
C    FINC-- FLOATING INCREMENT FOR XLABEL                               PLOT  46
C    IDAT-- COUNTER OVER ROWS OF X AND Y                                PLOT  47
C    ILAB-- INTEGER LABEL FLAG                                          PLOT  48
C    ITEM(101)--ITEMS TO BE PRINTED INITIALLY BLANK                     PLOT  49
C    LINE-- COUNTER OF CURRENT LINE                                     PLOT  50
C     NLPP-- NUMBER OF LINES PLOTTED PER PAGE                           PLOT  51
C    SPANX-- RANGE OF X                                                 PLOT  52
C    SPANY-- RANGE OF Y                                                 PLOT  53
C    VALUE-- DATUM PRINTED ON Y AXIS                                    PLOT  54
C    XLABEL(11)--X-AXIS LABELS                                          PLOT  55
C     XNOW-- CURRENT VALUE OF X VECTOR NOW BEING WORKED ON              PLOT  56
C     YNOW-- CURRENT VALUE OF Y MATRIX NOW BEING WORKED ON              PLOT  57
C                                                                       PLOT  58
      REAL ITEM(101), MINUS, PTID(100), YMINV(2), YMAXV(2)              PLOT  59
      DIMENSION XLABL(11), X(1), Y(LONG,NY), XS(10), YS(10),            PLOT  60
     *    IND(2)                                                        PLOT  61
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             PLOT  62
      COMMON /ACCUR/ XMAXN                                              PLOT  63
      COMMON /PLTCHR/ PTID, ITEM                                        PLOT  64
      DATA MINUS, DD, STAR, BLANK, AIE /1H-,1H.,1H*,1H ,1HI/            PLOT  65
      DATA XS /0.,0.,0.,0.,0.,4.,2.,4.,1.,0./                           PLOT  66
      DATA YS /0.,0.,0.,0.,2.,4.,3.,4.,7.,2./                           PLOT  67
C                                                                       PLOT  68
      JY = 1                                                            PLOT  69
      KODE = 1                                                          PLOT  70
      INDEX = 0                                                         PLOT  71
      LINE = 1                                                          PLOT  72
      IY0L = 0                                                          PLOT  73
      IX0L = 0                                                          PLOT  74
      IDAT = 1                                                          PLOT  75
      NLPP = 53                                                         PLOT  76
      XMAX = XA                                                         PLOT  77
      YMAX = YA                                                         PLOT  78
      YMIN = YI                                                         PLOT  79
      XMIN = XI                                                         PLOT  80
      NPOI2 = NY*NPOI                                                   PLOT  81
      NPOIRL = MX + NPOI - 1                                            PLOT  82
      LXFLG = 0                                                         PLOT  83
      LYFLG = 0                                                         PLOT  84
C                                                                       PLOT  85
C DETERMINE MINIMUM AND MAXIMUM OF X AXIS, IF NECESSARY                 PLOT  86
C                                                                       PLOT  87
      IF (XMAX-XMIN) 60, 10, 60                                         PLOT  88
  10  XMAX = -XMAXN                                                     PLOT  89
      XMIN = XMAXN                                                      PLOT  90
      DO 50 I=MX,NPOIRL                                                 PLOT  91
        IF (X(I)-XMAX) 30, 20, 20                                       PLOT  92
  20    XMAX = X(I)                                                     PLOT  93
  30    IF (X(I)-XMIN) 40, 40, 50                                       PLOT  94
  40    XMIN = X(I)                                                     PLOT  95
  50  CONTINUE                                                          PLOT  96
C                                                                       PLOT  97
C     CALL SERCH TO INITIALIZE; IF DESIRED SET YMAX AND YMIN            PLOT  98
C                                                                       PLOT  99
C                                                                       PLOT 100
C     IF WE HAVE A 1D PLOT NO NEED FOR SERCH                            PLOT 101
C                                                                       PLOT 102
      IF (IABS(ID)-3) 60, 100, 60                                       PLOT 103
  60  DO 70 I=1,NY                                                      PLOT 104
        IY = I                                                          PLOT 105
        CALL SERCH(Y(1,IY), NPOIRL, INDEX, KODE, IY, MX)                PLOT 106
        YMAXV(I) = Y(INDEX,I)                                           PLOT 107
        YMINV(I) = Y(KODE,I)                                            PLOT 108
        IND(I) = INDEX                                                  PLOT 109
        IF (NY.EQ.1) GO TO 80                                           PLOT 110
        KODE = 1                                                        PLOT 111
  70  CONTINUE                                                          PLOT 112
      YMAX = AMAX1(YMAXV(1),YMAXV(2))                                   PLOT 113
      YMIN = AMIN1(YMINV(1),YMINV(2))                                   PLOT 114
      IF (YMAX.NE.YMAXV(1)) JY = 2                                      PLOT 115
      YNOW = YMAX                                                       PLOT 116
      YNOW1 = AMIN1(YMAXV(1),YMAXV(2))                                  PLOT 117
      GO TO 100                                                         PLOT 118
  80  YNOW = Y(INDEX,1)                                                 PLOT 119
      YNOW1 = -XMAXN                                                    PLOT 120
      IF (YMAX-YMIN) 100, 90, 100                                       PLOT 121
  90  YMAX = Y(INDEX,1)                                                 PLOT 122
      YMIN = Y(KODE,1)                                                  PLOT 123
C                                                                       PLOT 124
C DETERMINE RANGE AND INCREMENT OF BOTH AXES                            PLOT 125
 100  P = ALOG10(XMAX-XMIN+.0001)                                       PLOT 126
      IP = IFIX(P)                                                      PLOT 127
      IF (P.LT.0.) IP = IP - 1                                          PLOT 128
      SPAN = (XMAX-XMIN)/(10.**IP)                                      PLOT 129
      IF (SPAN.LT.5.0) GO TO 110                                        PLOT 130
      DEL = 10.**IP                                                     PLOT 131
      GO TO 130                                                         PLOT 132
 110  IF (SPAN.LT.2.0) GO TO 120                                        PLOT 133
      DEL = (10.**IP)/2.0                                               PLOT 134
      GO TO 130                                                         PLOT 135
 120  DEL = (10.**IP)/5.0                                               PLOT 136
 130  IF (XMAX.LT.0.) GO TO 160                                         PLOT 137
      XMAX = FLOAT(IFIX(XMAX/DEL+.9999))*DEL                            PLOT 138
      IF (XMIN) 150, 140, 140                                           PLOT 139
 140  XMIN = FLOAT(IFIX(XMIN/DEL))*DEL                                  PLOT 140
      GO TO 170                                                         PLOT 141
 150  XMIN = FLOAT(IFIX(XMIN/DEL-.9999))*DEL                            PLOT 142
      GO TO 170                                                         PLOT 143
 160  XMAX = FLOAT(IFIX(XMAX/DEL))*DEL                                  PLOT 144
      XMIN = FLOAT(IFIX(XMIN/DEL+.9999))*DEL                            PLOT 145
 170  SPANX = XMAX - XMIN                                               PLOT 146
      NXUNIT = SPANX/DEL + .00001                                       PLOT 147
      IF (NXUNIT.LE.10) GO TO 180                                       PLOT 148
      LXFLG = LXFLG + 1                                                 PLOT 149
      IF (LXFLG-1) 830, 100, 830                                        PLOT 150
 180  NXP1 = NXUNIT + 1                                                 PLOT 151
      FINC = DEL                                                        PLOT 152
      DELX = SPANX/(100.-XS(NXUNIT))                                    PLOT 153
      IXS = IFIX(XS(NXUNIT)/2.+.1)                                      PLOT 154
      XSF = IXS                                                         PLOT 155
      IF (ID.LT.0) GO TO 190                                            PLOT 156
      IF ((XMIN.GE.0.) .OR. (XMAX.LE.0.)) GO TO 190                     PLOT 157
      IX0L = -XMIN/DELX + 1.5                                           PLOT 158
      IX0L = IX0L + IXS                                                 PLOT 159
C                                                                       PLOT 160
 190  IF (ID.EQ.3) GO TO 280                                            PLOT 161
      P = ALOG10(YMAX-YMIN+.0001)                                       PLOT 162
      IP = IFIX(P)                                                      PLOT 163
      IF (P.LT.0.) IP = IP - 1                                          PLOT 164
      SPAN = (YMAX-YMIN)/(10.**IP)                                      PLOT 165
      IF (SPAN.LT.5.0) GO TO 200                                        PLOT 166
      DEL = 10.**IP                                                     PLOT 167
      GO TO 220                                                         PLOT 168
 200  IF (SPAN.LT.2.0) GO TO 210                                        PLOT 169
      DEL = (10.**IP)/2.0                                               PLOT 170
      GO TO 220                                                         PLOT 171
 210  DEL = (10.**IP)/5.0                                               PLOT 172
 220  IF (YMAX.LT.0.) GO TO 250                                         PLOT 173
      YMAX = FLOAT(IFIX(YMAX/DEL+.9999))*DEL                            PLOT 174
      IF (YMIN) 240, 230, 230                                           PLOT 175
 230  YMIN = FLOAT(IFIX(YMIN/DEL))*DEL                                  PLOT 176
      GO TO 260                                                         PLOT 177
 240  YMIN = FLOAT(IFIX(YMIN/DEL-.9999))*DEL                            PLOT 178
      GO TO 260                                                         PLOT 179
 250  YMAX = FLOAT(IFIX(YMAX/DEL))*DEL                                  PLOT 180
      YMIN = FLOAT(IFIX(YMIN/DEL+.9999))*DEL                            PLOT 181
 260  SPANY = YMAX - YMIN                                               PLOT 182
      NYUNIT = SPANY/DEL + .00001                                       PLOT 183
      IF (NYUNIT.LE.10) GO TO 270                                       PLOT 184
      LYFLG = LYFLG + 1                                                 PLOT 185
      IF (LYFLG-1) 830, 190, 830                                        PLOT 186
 270  DELY = SPANY/(FLOAT(NLPP-1)-YS(NYUNIT))                           PLOT 187
      DELTY = DEL                                                       PLOT 188
      VALUE = YMAX                                                      PLOT 189
      YININC = (NLPP-1)/NYUNIT                                          PLOT 190
      LININC = IFIX(YININC)                                             PLOT 191
      IYLAB = LININC                                                    PLOT 192
      IF (ID.LT.0) GO TO 290                                            PLOT 193
      IF ((YMAX.LE.0.) .OR. (YMIN.GE.0.)) GO TO 290                     PLOT 194
      IY0L = YMAX/DELY + 1.5                                            PLOT 195
      GO TO 290                                                         PLOT 196
 280  VALUE = 0.0                                                       PLOT 197
      DELTY = 0.0                                                       PLOT 198
      LININC = 27                                                       PLOT 199
      IYLAB = 1                                                         PLOT 200
C                                                                       PLOT 201
C LABEL PLOT AT TOP                                                     PLOT 202
C                                                                       PLOT 203
 290  XLABL(1) = XMIN                                                   PLOT 204
      DO 300 I=1,NXUNIT                                                 PLOT 205
        XLABL(I+1) = XLABL(I) + FINC                                    PLOT 206
 300  CONTINUE                                                          PLOT 207
      GO TO (830, 830, 830, 310, 320, 330, 340, 350, 360, 370),         PLOT 208
     *    NXUNIT                                                        PLOT 209
 310  WRITE (LPRINT,9999) (XLABL(I),I=1,NXP1)                           PLOT 210
      WRITE (LPRINT,9998)                                               PLOT 211
9999  FORMAT (1H0, F20.4, 4F25.4)                                       PLOT 212
9998  FORMAT (1H , 14X, 36H*.************************.*********,        PLOT 213
     *    50H***************.************************.*********,        PLOT 214
     *    17H***************.*)                                         PLOT 215
      GO TO 380                                                         PLOT 216
 320  WRITE (LPRINT,9997) (XLABL(I),I=1,NXP1)                           PLOT 217
      WRITE (LPRINT,9996)                                               PLOT 218
9997  FORMAT (1H0, 6F20.4)                                              PLOT 219
9996  FORMAT (1H , 14X, 36H*.*******************.**************,        PLOT 220
     *    50H*****.*******************.*******************.****,        PLOT 221
     *    17H***************.*)                                         PLOT 222
      GO TO 380                                                         PLOT 223
 330  WRITE (LPRINT,9995) (XLABL(I),I=1,NXP1)                           PLOT 224
      WRITE (LPRINT,9994)                                               PLOT 225
9995  FORMAT (1H0, 6X, 7F16.4)                                          PLOT 226
9994  FORMAT (1H , 14X, 36H***.***************.***************.,        PLOT 227
     *    50H***************.***************.***************.**,        PLOT 228
     *    17H*************.***)                                         PLOT 229
      GO TO 380                                                         PLOT 230
 340  WRITE (LPRINT,9993) (XLABL(I),I=1,NXP1)                           PLOT 231
      WRITE (LPRINT,9992)                                               PLOT 232
9993  FORMAT (1H0, 7X, 8F14.4)                                          PLOT 233
9992  FORMAT (1H , 14X, 36H**.*************.*************.*****,        PLOT 234
     *    50H********.*************.*************.*************,        PLOT 235
     *    17H.*************.**)                                         PLOT 236
      GO TO 380                                                         PLOT 237
 350  WRITE (LPRINT,9991) (XLABL(I),I=1,NXP1)                           PLOT 238
      WRITE (LPRINT,9990)                                               PLOT 239
9991  FORMAT (1H0, 10X, 9F12.4)                                         PLOT 240
9990  FORMAT (1H , 14X, 36H***.***********.***********.********,        PLOT 241
     *    50H***.***********.***********.***********.**********,        PLOT 242
     *    17H*.***********.***)                                         PLOT 243
      GO TO 380                                                         PLOT 244
 360  WRITE (LPRINT,9989) (XLABL(I),I=1,NXP1)                           PLOT 245
      WRITE (LPRINT,9988)                                               PLOT 246
9989  FORMAT (1H0, 8X, 10F11.3)                                         PLOT 247
9988  FORMAT (1H , 14X, 36H*.**********.**********.**********.*,        PLOT 248
     *    50H*********.**********.**********.**********.*******,        PLOT 249
     *    17H***.**********.**)                                         PLOT 250
      GO TO 380                                                         PLOT 251
 370  WRITE (LPRINT,9987) (XLABL(I),I=1,NXP1)                           PLOT 252
      WRITE (LPRINT,9986)                                               PLOT 253
9987  FORMAT (1H0, 9X, 11F10.3)                                         PLOT 254
9986  FORMAT (1H , 14X, 36H*.*********.*********.*********.****,        PLOT 255
     *    50H*****.*********.*********.*********.*********.****,        PLOT 256
     *    17H*****.*********.*)                                         PLOT 257
 380  CONTINUE                                                          PLOT 258
C                                                                       PLOT 259
C     PREPARE PLOT                                                      PLOT 260
C                                                                       PLOT 261
      KODE = 0                                                          PLOT 262
 390  IF (LINE-NLPP) 400, 400, 740                                      PLOT 263
 400  IF (IDAT-1) 410, 470, 410                                         PLOT 264
 410  IF (IABS(ID)-3) 420, 470, 420                                     PLOT 265
 420  IF (LFLAG) 430, 430, 470                                          PLOT 266
 430  IF (IDAT-NPOI2) 440, 440, 470                                     PLOT 267
 440  CALL SERCH(Y(1,JY), NPOIRL, INDEX, KODE, JY, MX)                  PLOT 268
      IF (INDEX.EQ.(-1)) GO TO 450                                      PLOT 269
      IND(JY) = INDEX                                                   PLOT 270
      YNOW = Y(INDEX,JY)                                                PLOT 271
      IF (YNOW-YNOW1) 460, 470, 470                                     PLOT 272
 450  JY = NY - JY + 1                                                  PLOT 273
      YNOW = YNOW1                                                      PLOT 274
      YNOW1 = -XMAXN                                                    PLOT 275
      GO TO 470                                                         PLOT 276
 460  JY = NY - JY + 1                                                  PLOT 277
      TEMP = YNOW                                                       PLOT 278
      YNOW = YNOW1                                                      PLOT 279
      YNOW1 = TEMP                                                      PLOT 280
 470  CONTINUE                                                          PLOT 281
C                                                                       PLOT 282
C     CHECK TO SEE IF ANY POINTS IN Y REMAIN                            PLOT 283
C                                                                       PLOT 284
      IF (IDAT-NPOI2) 480, 480, 600                                     PLOT 285
 480  IF (IABS(ID)-3) 490, 520, 490                                     PLOT 286
 490  IF (YMAX-YNOW) 570, 500, 500                                      PLOT 287
 500  IF (YNOW-YMIN) 570, 510, 510                                      PLOT 288
 510  I = (YMAX-YNOW)/DELY + 1.50001                                    PLOT 289
      GO TO 550                                                         PLOT 290
 520  IF (LINE-27) 550, 530, 550                                        PLOT 291
 530  DO 540 II=MX,NPOIRL                                               PLOT 292
        I = II                                                          PLOT 293
        CALL XITEM(I, XMAX, XMIN, ID, DELX, X, JY, XSF)                 PLOT 294
 540  CONTINUE                                                          PLOT 295
      GO TO 590                                                         PLOT 296
 550  IF (I-LINE) 580, 580, 560                                         PLOT 297
 560  LFLAG = 1                                                         PLOT 298
      GO TO 600                                                         PLOT 299
 570  LFLAG = 0                                                         PLOT 300
      GO TO 590                                                         PLOT 301
 580  LFLAG = 0                                                         PLOT 302
      INDEX = IND(JY)                                                   PLOT 303
      CALL XITEM(INDEX, XMAX, XMIN, ID, DELX, X, JY, XSF)               PLOT 304
 590  IDAT = IDAT + 1                                                   PLOT 305
      GO TO 390                                                         PLOT 306
C                                                                       PLOT 307
C PRINT ONE LINE OF PLOT                                                PLOT 308
C                                                                       PLOT 309
C     CHECK ID -- IF MINUS NO AXES,IF + AXES                            PLOT 310
C                IF ID=3 1D PLOT                                        PLOT 311
C                                                                       PLOT 312
 600  IF (ID) 700, 610, 610                                             PLOT 313
 610  IF (IABS(ID)-3) 620, 700, 620                                     PLOT 314
 620  IF (IX0L) 650, 650, 630                                           PLOT 315
 630  IF (ITEM(IX0L)-BLANK) 650, 640, 650                               PLOT 316
 640  ITEM(IX0L) = AIE                                                  PLOT 317
 650  IF (IY0L) 700, 700, 660                                           PLOT 318
 660  IF (LINE-IY0L) 700, 670, 700                                      PLOT 319
 670  DO 690 I=1,101                                                    PLOT 320
        IF (ITEM(I)-BLANK) 690, 680, 690                                PLOT 321
 680    ITEM(I) = MINUS                                                 PLOT 322
 690  CONTINUE                                                          PLOT 323
 700  CONTINUE                                                          PLOT 324
C                                                                       PLOT 325
      IF (IYLAB.NE.LININC) GO TO 710                                    PLOT 326
      WRITE (LPRINT,9985) VALUE, DD, (ITEM(J),J=1,101), DD,             PLOT 327
     *    VALUE                                                         PLOT 328
9985  FORMAT (1H , F12.4, 2X, 103A1, F13.4)                             PLOT 329
      VALUE = VALUE - DELTY                                             PLOT 330
      IYLAB = 0                                                         PLOT 331
      GO TO 720                                                         PLOT 332
 710  WRITE (LPRINT,9984) STAR, (ITEM(J),J=1,101), STAR                 PLOT 333
9984  FORMAT (1H , 14X, 103A1)                                          PLOT 334
 720  IYLAB = IYLAB + 1                                                 PLOT 335
      DO 730 I=1,101                                                    PLOT 336
        ITEM(I) = BLANK                                                 PLOT 337
 730  CONTINUE                                                          PLOT 338
      LINE = LINE + 1                                                   PLOT 339
      GO TO 390                                                         PLOT 340
C                                                                       PLOT 341
C LABEL PLOT ALONG THE BOTTOM                                           PLOT 342
C                                                                       PLOT 343
 740  GO TO (830, 830, 830, 750, 760, 770, 780, 790, 800, 810),         PLOT 344
     *    NXUNIT                                                        PLOT 345
 750  WRITE (LPRINT,9998)                                               PLOT 346
      WRITE (LPRINT,9999) (XLABL(I),I=1,NXP1)                           PLOT 347
      GO TO 820                                                         PLOT 348
 760  WRITE (LPRINT,9996)                                               PLOT 349
      WRITE (LPRINT,9997) (XLABL(I),I=1,NXP1)                           PLOT 350
      GO TO 820                                                         PLOT 351
 770  WRITE (LPRINT,9994)                                               PLOT 352
      WRITE (LPRINT,9995) (XLABL(I),I=1,NXP1)                           PLOT 353
      GO TO 820                                                         PLOT 354
 780  WRITE (LPRINT,9992)                                               PLOT 355
      WRITE (LPRINT,9993) (XLABL(I),I=1,NXP1)                           PLOT 356
      GO TO 820                                                         PLOT 357
 790  WRITE (LPRINT,9990)                                               PLOT 358
      WRITE (LPRINT,9991) (XLABL(I),I=1,NXP1)                           PLOT 359
      GO TO 820                                                         PLOT 360
 800  WRITE (LPRINT,9988)                                               PLOT 361
      WRITE (LPRINT,9989) (XLABL(I),I=1,NXP1)                           PLOT 362
      GO TO 820                                                         PLOT 363
 810  WRITE (LPRINT,9986)                                               PLOT 364
      WRITE (LPRINT,9987) (XLABL(I),I=1,NXP1)                           PLOT 365
 820  CONTINUE                                                          PLOT 366
      RETURN                                                            PLOT 367
 830  WRITE (LPRINT,9983)                                               PLOT 368
9983  FORMAT (22H0ERROR EXIT FROM PLOT.)                                PLOT 369
      STOP                                                              PLOT 370
      END                                                               PLOT 371
C                                                                       SERC   1
      SUBROUTINE SERCH(VEC, LEN, INDEX, KODE, JY, MX)                   SERC   2
C SERCH            FOR KYST                  JANUARY,1973               SERC   3
C     OCTOBER 10,1971     JAY R LEVINSOHN                               SERC   4
C MODIFIED FOR KYST BY J.SEERY               JANUARY,1973               SERC   5
C                                                                       SERC   6
C     SEARCH SERVES TWO FUNCTIONS                                       SERC   7
C        1 FIND THE MIN AND MAX OF SOME VECTOR VEC                      SERC   8
C        2 FIND THE NEXT BIGGEST ITEM IN SOME VECTOR VEC                SERC   9
C      CAN KEEP TRACK OF TWO VECTORS (OF THE SAME LENGTH) AT A TIME     SERC  10
C     PARAMETERS                                                        SERC  11
C        VEC-- VECTOR TO BE SEARCHED                                    SERC  12
C        LEN-- NUMBER OF ELEMENTS IN VEC                                SERC  13
C        MX--      WORK WITH ELEMENTS MX THROUGH LEN                    SERC  14
C        INDEX-- POINTER TO NEXT BIGGEST ITEM IN VEC                    SERC  15
C        KODE-- OP. CODE 1=FIND MIN AND MAX,0= FIND NEXT BIGGEST ITEM   SERC  16
C     IF KODE=1 THEN POINTER TO MIN IS RETURNED IN KODE WHILE POINTER   SERC  17
C     TO MAX IS RETURNED IN INDEX                                       SERC  18
C     ROUTINE ASSUMES IT WILL BE CALLED TO FIRST FIND MIN AND MAX OF    SERC  19
C     VECTOR AND THEN CALLED ITERATIVELY TO FIND THE SECOND LARGEST     SERC  20
C     ELEMENT, THIRD LARGEST ELEMENT,, ETC.... UNTIL  FINAL ELEMENT HAS SERC  21
C     HAS BEEN EXTRACTED.  ONE CAN TEST THE FINAL ELEMENT AGAINST THE   SERC  22
C     MIN TO CHECK FOR ACCURATE COMPLETION.                             SERC  23
C                                                                       SERC  24
      DIMENSION VEC(1)                                                  SERC  25
      DIMENSION FMIN(2), IMIN(2), INDB4(2), KNT(2)                      SERC  26
      COMMON /DEVICE/ LREAD, LPRINT, LPUNCH                             SERC  27
      COMMON /ACCUR/ XMAXN                                              SERC  28
      COMMON /SCOM/ FMIN, IMIN, INDB4, KNT                              SERC  29
C                                                                       SERC  30
      IF (KODE) 170, 60, 10                                             SERC  31
C                                                                       SERC  32
C     OP CODE =1 THEN FIND MIN AND MAX                                  SERC  33
C                                                                       SERC  34
  10  FMAX = -XMAXN                                                     SERC  35
      FMIN(JY) = XMAXN                                                  SERC  36
      KNT(JY) = MX                                                      SERC  37
      DO 50 I=MX,LEN                                                    SERC  38
        IF (VEC(I)-FMAX) 30, 20, 20                                     SERC  39
  20    FMAX = VEC(I)                                                   SERC  40
        INDEX = I                                                       SERC  41
  30    IF (VEC(I)-FMIN(JY)) 40, 40, 50                                 SERC  42
  40    FMIN(JY) = VEC(I)                                               SERC  43
        KODE = I                                                        SERC  44
  50  CONTINUE                                                          SERC  45
C                                                                       SERC  46
C     SAVE MIN IN SAVE,MIN POINTER IN IMIN,RESET FMIN, SAVE MAX         SERC  47
C                                                                       SERC  48
      IMIN(JY) = KODE                                                   SERC  49
      FMIN(JY) = -XMAXN                                                 SERC  50
      INDB4(JY) = INDEX                                                 SERC  51
      RETURN                                                            SERC  52
C                                                                       SERC  53
C     OP. CODE NOT EQUAL 1 FIND NEXT LARGEST ELE IN VEC                 SERC  54
C                                                                       SERC  55
  60  IF (KNT(JY).LT.LEN) GO TO 70                                      SERC  56
      INDEX = -1                                                        SERC  57
      GO TO 160                                                         SERC  58
  70  INDEX = INDB4(JY)                                                 SERC  59
      FMB4 = VEC(INDEX)                                                 SERC  60
      VEC(INDEX) = FMB4 + 1.0                                           SERC  61
      KNT(JY) = KNT(JY) + 1                                             SERC  62
      DO 110 I=MX,LEN                                                   SERC  63
        IF (VEC(I)-FMIN(JY)) 110, 80, 80                                SERC  64
  80    IF (VEC(I)-FMB4) 90, 90, 110                                    SERC  65
  90    IF (I-INDB4(JY)) 100, 110, 100                                  SERC  66
 100    FMIN(JY) = VEC(I)                                               SERC  67
        INDEX = I                                                       SERC  68
 110  CONTINUE                                                          SERC  69
C                                                                       SERC  70
C     SUBTRACT 1 FROM EACH ELEMENT OF VEC TO GET IT BACK AS BEFORE      SERC  71
C     DON'T DO SUBTRACTION UNLESS AT LAST ITEM OF VEC                   SERC  72
C                                                                       SERC  73
      IF (KNT(JY)-LEN) 150, 120, 150                                    SERC  74
 120  DO 140 I=MX,LEN                                                   SERC  75
        IF (I-INDEX) 130, 140, 130                                      SERC  76
 130    VEC(I) = VEC(I) - 1.0                                           SERC  77
 140  CONTINUE                                                          SERC  78
C                                                                       SERC  79
C     SAVE INDEX AND RESET FMIN                                         SERC  80
C                                                                       SERC  81
 150  INDB4(JY) = INDEX                                                 SERC  82
      FMIN(JY) = -XMAXN                                                 SERC  83
 160  RETURN                                                            SERC  84
C                                                                       SERC  85
C     ERROR ROUTINE ENTER HERE IF KODE<0                                SERC  86
C                                                                       SERC  87
 170  WRITE (LPRINT,9999) KODE                                          SERC  88
9999  FORMAT (24H ERROR IN SERCH, KODE= ,I4)                            SERC  89
      RETURN                                                            SERC  90
      END                                                               SERC  91
C                                                                       XITE   1
      SUBROUTINE XITEM(INDEX, XMAX, XMIN, ID, DELX, X, JX, XSF)         XITE   2
C XITEM            FOR KYST                  JANUARY,1973               XITE   3
C     JAY R LEVINSOHN     OCTOBER 31,1971                               XITE   4
C ADAPTED FOR KYST BY J.SEERY                JANUARY,1973               XITE   5
C                                                                       XITE   6
C     ITS FUNCTION IS TO PLACE  ELEMENTS OF THE X VECTOR IN THE         XITE   7
C     PRINT LINE.                                                       XITE   8
C                                                                       XITE   9
C                                                                       XITE  10
      REAL ITEM(101), PTID(100), X(1)                                   XITE  11
      COMMON /PLTCHR/ PTID, ITEM                                        XITE  12
C                                                                       XITE  13
      XNOW = X(INDEX)                                                   XITE  14
      IF (XNOW-XMAX) 10, 10, 60                                         XITE  15
  10  IF (XNOW-XMIN) 60, 20, 20                                         XITE  16
  20  J = (XNOW-XMIN)/DELX + 1.5 + XSF                                  XITE  17
C                                                                       XITE  18
C     DETERMINE LABEL FOR OBSERVATION                                   XITE  19
C                                                                       XITE  20
      IF (IABS(ID)-3) 30, 50, 30                                        XITE  21
  30  IF (IABS(ID)-2) 40, 50, 40                                        XITE  22
  40  IF (ITEM(J).NE.PTID(2)) ITEM(J) = PTID(JX)                        XITE  23
      GO TO 60                                                          XITE  24
  50  ITEM(J) = PTID(INDEX)                                             XITE  25
  60  CONTINUE                                                          XITE  26
      RETURN                                                            XITE  27
      END                                                               XITE  28