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

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

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

C For explanation of the method this software implements, see
C Carroll,J.D. & Chang,J.J. (1972). SIMULES: SIMUltaneous Linear
C Equation Scaling. Proceedings of the 80th Annual Convention of the
C American Psychological Association, 7, 11-12.
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

$      FORTRAN
C     SIMULES- SIMULTANEOUS LINEAR EQUATION SCALING
C     6/28/72, BY CARROLL AND CHANG
C     NSUB- NUMBER OF SUBJECTS.  (NO LIMIT)
C     N- NUMBER OF STIMULI.  (MAX.=60)
C     IDA- 0-ANALOGY DATA, 4 ITEMS
C          1- BISECTION DATA,  3 ITEMS.
C          2-DIFFERENCE DATA, 2 ITEMS
C     KPU  =0, DO NOT PUNCH COORDINATES
C          NE 1, PUNCH COORDINATES IN KPU DIMENSIONS.
C     KPL  =0, DO NOT PLOT
C          NE 1, PLOT IN KPL DIMENSIONS.
C     THERE IS NO LIMIT ON NUMBER OF JUDGMENTS.
      DIMENSION R(60,60),Q(60,60),COE(4),IDEN(4),A(60),TIT(12)
      DIMENSION FSUB(2),FB(3),BNUM(60)
      DATA FB(1)/13H(7HSUBJECTI3)/
 98   READ100,NSUB,N,KPU,KPL,IDA,LAB
      IF(NSUB.EQ.0) GOTO99
      N1=N-1
      READ101,TIT
      IF(LAB.EQ.0) GOTO18
      READ107,(BNUM(I),I=1,N)
      GOTO19
C     PIBCI IS A SYSTEM ROUTINE WHICH CONVERTS INTEGER TO BCD CHARACTERS
C     , FOR LABELLING ON PLOTS.
 18   DO16I=1,N
 16   CALL PIBCI(I,BNUM(I),4H(I2))
C     READ IN DATA BY ROW- EACH CARD GIVES THE STIMULUS NUMBERS
 19   NL=60
      PRINT112,TIT
      PRINT111,NSUB,N, KPU,KPL,IDA,LAB
      IF(KPU.NE.0) PUNCH110
      PUNCH104
      COE(1)=1.
      IF(IDA-1)2,3,4
 2    ITEM=4
      COE(2)=-1.
      COE(3)=-1.
      COE(4)=1.
      GOTO5
 3    ITEM=3
      COE(2)=1.
      COE(3)=-2.
      GOTO5
 4    ITEM=2
      COE(2)=1.
 5    DO50IS=1,NSUB
      PRINT108,IS
      DO7I=1,N
      DO7J=1,N
 7    R(I,J)=0.
      L=0
 1    DO11J=1,N
 11   A(J)=0.
      READ100,(IDEN(I),I=1,ITEM)
      IF(IDEN(1).EQ.0)GOTO49
      L=L&1
      DO6I=1,ITEM
      M=IDEN(I)
 6    A(M)=A(M)&COE(I)
      DO9I=1,ITEM
      M1=IDEN(I)
      DO9J=1,I
      M2=IDEN(J)
      R(M1,M2)=R(M1,M2)&A(M1)*A(M2)
 9    R(M2,M1)=R(M1,M2)
      GOTO1
 49   MX=1
      CALL MTEVV(R,NL,NL,Q,NL,NL,N,MX)
      DO17I=1,N1
      M=N-I
 17   A(I)=R(M,M)
      PRINT105
      PRINT109,(A(I),I=1,N1)
      CALL SVAF(A,N1)
      DO10I=1,N1
      M=N-I
      P=1./SQRT(A(I))
      DO10J=1,N
 10   R(J,I)=Q(J,M)*P
      PRINT106
      DO62J=1,N1
 62   PRINT102,J,(R(I,J),I=1,N)
      IF(KPL.EQ.0) GOTO15
      CALL FORCNV(IS,1,FSUB,FB,DUM1,DUM2)
      CALL DSUBJ(R,N,KPL,NL,TIT,BNUM,ISCALE,FSUB,2)
 15   IF(KPU.EQ.0) GOTO12
      DO63J=1,KPU
 63   PUNCH103,IS,J,(R(I,J),I=1,N)
C     COMPUTE SCALAR PRODUCT FOR EACH SUBJECT
 12   DO13I=1,N
      DO13J=1,I
      Q(I,J)=0.
      DO13M=1,N1
 13   Q(I,J)=Q(I,J)&R(I,M)*R(J,M)
      DO14I=1,N
 14   PUNCH103,IS,I,(Q(I,J),J=1,I)
 50   CONTINUE
      GOTO98
 100  FORMAT(20I4)
 101  FORMAT(12A6)
 102  FORMAT(I3,3X,10E12.4/(6X,10E12.4))
 103  FORMAT(I2,I3,5E13.5/(5X,5E13.5))
 104  FORMAT(21HSCALAR PRODUCT MATRIX/11H(5X,5E13.5))
 105  FORMAT(18H0EIGEN VALUES OF R)
 106  FORMAT(23H0WEIGHTED EIGEN VECTORS)
 107  FORMAT(20A4)
 108  FORMAT(8H0SUBJECTI3)
 109  FORMAT(6X,10E12.4)
 110  FORMAT(11HCOORDINATES)
 111  FORMAT(28H0NSUB  N  KPU  KPL  IDA  LAB/3I4,3I5)
 112  FORMAT(1H112A6)
 99   STOP
      END
$      FORTRAN
      SUBROUTINE SVAF(A,N)
      DIMENSION A(1),P(60)
      SUM=0.
      DO2I=1,N
      P(I)=1./(A(I)**2)
 2    SUM=SUM&P(I)
      DO3I=1,N
 3    P(I)=P(I)/SUM
      PRINT100,(I,I=1,N)
      PRINT101,(P(I),I=1,N)
      SUM=P(1)
      DO4I=2,N
      SUM=SUM&P(I)
 4    P(I)=SUM
      PRINT102,(I,I=1,N)
      PRINT101,(P(I),I=1,N)
 100  FORMAT(73H0PROPORTION OF VARIANCE ACCOUNTED FOR BY EACH FACTOR (IN
     1 SCALAR PRODUCTS)/(I9,9I12))
 101  FORMAT(6X,10E12.4)
 102  FORMAT(48H0CUMULATIVE PROPORTION OF VARIANCE ACCOUNTED FOR/(I9,9I1
     12))
      RETURN
      END
$      FORTRAN
      SUBROUTINE DSUBJ(X,N,NF,NL,TIT,BNUM,FSUB,NW)
C     PROGRAM TO PLOT N STIMULUS POINTS X IN NF DIMENSIONAL SPACE
      DIMENSION X(NL,NF),BNUM(1),TIT(1),FSUB(1)
      DIMENSION FFM(2),FI(2),FJ(2),FDM(4),FFDM(3)
      DATA FDM(1)/20H(I2,12H DIMENSIONAL)/
      DATA FFM(1)/10H(4HDIM.I3)/
C     PLOTTING AREA IS 460 BY 460, ORIGIN AT CENTER OF SQUARE
      NNF=NF-1
      CALL FORCNV(NF,1,FFDM,FDM,DUM1,DUM2)
 5    XMAX=0.
      DO10I=1,NF
      CALL MAXIS(X(1,I),SCL,N)
      IF(SCL.GT.XMAX) XMAX=SCL
 10   CONTINUE
      XSCALE=460./XMAX
 6    DO260I=1,NNF
      JJ=I&1
      DO260J=JJ,NF
 34   CALL ROLL
      CALL REFRNC(-512,512,1024,1024)
C     PRINT TITLE
      CALL TEXT(-460,512,TIT,0,12)
      CALL TEXT(-460,460,FFDM,0,3)
C     LABEL NF DIMENSIONAL SOLUTION
C     PRINT SUBTITLE
      CALL TEXT (-460,500,FSUB,0,NW)
C     DRAW POINTS
      DO 22 K=1,N
      IAX=X(K,I)*XSCALE
      IF(NF.EQ.1)GOTO23
      IAY=X(K,J)*XSCALE
      GOTO24
 23   IAY=0
 24   CALL TEXT (IAX,IAY,BNUM(K),0,1)
 22   CONTINUE
      IF(NF.EQ.1)GOTO50
C     DRAW AXES
      CALL DVR1(-460,0,460,0,1)
      CALL DVR1(0,-460,0,460,1)
C     LABEL AXES
      CALL FORCNV(J,1,FJ,FFM,DUM1,DUM2)
      CALL TEXT(-10,467,FJ,0,2)
      CALL FORCNV(I,1,FI,FFM,DUM1,DUM2)
      CALL TEXT(440,7,FI,0,2)
 260  CONTINUE
 50   RETURN
      END
$      FORTRAN
      SUBROUTINE MAXI(X,Y,Z,N)
C     MAXI FINDS THE MAGNITUDE OF ARRAYS X AND Y EACH OF LENGTH N.
C     ON RETURN ABSOLUTE MAXIMUM IS STORED IN Z.
      DIMENSION X(1),Y(1)
      Z=ABS(X(1))
      DO10I=1,N
      AXO=ABS(X(I))
      IF(AXO.GT.Z)  Z=AXO
      AYO=ABS(Y(I))
      IF(AYO.GT.Z)  Z=AYO
 10   CONTINUE
      RETURN
      ENTRY MAXIS(X,Z,N)
      Z=ABS(X(1))
      DO11I=2,N
      AXO=ABS(X(I))
      IF(AXO.GT.Z) Z=AXO
 11   CONTINUE
      RETURN
      END
$      FORTRAN
$      INCODE  IBMF
*EIGENJ      EIGENVALUE AND EIGENVECTOR SUBROUTINE                      EIGJ0020
*                              CD600D4.009      DATE 05/05/65           EIGJ0030
      SUBROUTINE EIGENJ (A,V,NN,IV,E,I3,I4)
      DIMENSION A(5050),V( I3, I4),E(NN)
      LOGICAL IV                                                        EIGJ0060
*     EIGENJ FINDS THE EIGENVALUES AND EIGENVECTORS OF                  EIGJ0070
*     A SYMMETRIC MATRIX (A) USING A MODIFIED THRESHOLD JACOBI METHOD   EIGJ0080
*     ONLY THE UPPER TRIANGLE IS STORED                                 EIGJ0090
*     DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER         EIGJ0100
*     DIMENSION A(K=M(M&1)/2),V(M,M),E(M) WHERE M=MAXIMUM ORDER         EIGJ0110
      N=NN                                                              EIGJ0120
      IND=1                                                             EIGJ0130
      NM2=N-2                                                           EIGJ0140
      AMAX=0.0                                                          EIGJ0150
      NM=N-1                                                            EIGJ0160
      IF(.NOT.IV)GO TO 5                                                EIGJ0170
*     SET UP VECTOR MATRIX                                              EIGJ0180
    1 DO 3 I=1,N                                                        EIGJ0190
      DO 2 J=1,N                                                        EIGJ0200
    2 V(I,J)=0.0                                                        EIGJ0210
    3 V(I,I)=1.0                                                        EIGJ0220
    5 K=2                                                               EIGJ0230
      DO 28 I=1,NM                                                      EIGJ0240
      IP=I&1                                                            EIGJ0250
      DO 27 J=IP,N                                                      EIGJ0260
      Y=A(K)                                                            EIGJ0270
      X=ABS(Y)                                                          EIGJ0280
      IF(X-AMAX)27,27,6                                                 EIGJ0290
    6 GO TO (7,8),IND                                                   EIGJ0300
    7 AMAX=X                                                            EIGJ0310
      GO TO 27                                                          EIGJ0320
*     TRANSFORMATION SETUP FOLLOWS                                      EIGJ0330
    8 II=K&I-J                                                          EIGJ0340
      JJ=(J*(N&N-J&3))/2-N                                              EIGJ0350
      ITEST=1                                                           EIGJ0360
      X=A(II)-A(JJ)                                                     EIGJ0370
      IF(X)10,11,10                                                     EIGJ0380
   10 X=Y/X                                                             EIGJ0390
      Y=.5*X                                                            EIGJ0400
      IF(ABS(Y)-.414213562)14,11,11                                     EIGJ0410
   11 C=.707106781                                                      EIGJ0420
      IF(Y)12,12,13                                                     EIGJ0430
   12 S=-C                                                              EIGJ0440
      GO TO 15                                                          EIGJ0450
   13 S=C                                                               EIGJ0460
      GO TO 15                                                          EIGJ0470
   14 X=Y*Y                                                             EIGJ0480
      C=1.&X                                                            EIGJ0490
      S=2.*Y/C                                                          EIGJ0500
      C=(1.-X)/C                                                        EIGJ0510
   15 X=S*S                                                             EIGJ0520
      Y=C*C                                                             EIGJ0530
      XY=S*C                                                            EIGJ0540
      AXY=2.*A(K)*XY                                                    EIGJ0550
      Z=A(II)*Y&AXY&A(JJ)*X                                             EIGJ0560
      W=A(II)*X-AXY&A(JJ)*Y                                             EIGJ0570
      A(K)=A(K)*(Y-X)&XY*(A(JJ)-A(II))                                  EIGJ0580
      A(II)=Z                                                           EIGJ0590
      A(JJ)=W                                                           EIGJ0600
      IF(NM2)24,24,16                                                   EIGJ0610
   16 IT=I                                                              EIGJ0620
      JT=J                                                              EIGJ0630
*     TRANSFORM A                                                       EIGJ0640
      DO 23 M=1,NM2                                                     EIGJ0650
      NP=N-M                                                            EIGJ0660
      IF(JT-K)20,17,18                                                  EIGJ0670
   17 IT=IT&1                                                           EIGJ0680
      JT=JT&NP                                                          EIGJ0690
      ITEST=2                                                           EIGJ0700
   18 IF(IT-K)20,19,19                                                  EIGJ0710
   19 IT=IT&1                                                           EIGJ0720
      JT=JT&1                                                           EIGJ0730
      ITEST=3                                                           EIGJ0740
   20 X=A(IT)*C&A(JT)*S                                                 EIGJ0750
      A(JT)=A(JT)*C-A(IT)*S                                             EIGJ0760
      A(IT)=X                                                           EIGJ0770
      GO TO (21,22,23),ITEST                                            EIGJ0780
   21 IT=IT&NP                                                          EIGJ0790
      JT=JT&NP                                                          EIGJ0800
      GO TO 23                                                          EIGJ0810
   22 IT=IT&1                                                           EIGJ0820
      JT=JT&NP-1                                                        EIGJ0830
   23 CONTINUE                                                          EIGJ0840
   24 IF(.NOT.IV)GO TO 27                                               EIGJ0850
*     TRANSFORM V                                                       EIGJ0860
   25 DO 26 M=1,N                                                       EIGJ0870
      X=V(M,I)*C&V(M,J)*S                                               EIGJ0880
      V(M,J)=V(M,J)*C-V(M,I)*S                                          EIGJ0890
   26 V(M,I)=X                                                          EIGJ0900
   27 K=K&1                                                             EIGJ0910
   28 K=K&1                                                             EIGJ0920
*     SEQUENCE TO ADJUST THRESHOLD FOR EACH SWEEP                       EIGJ0930
      GO TO (29,31),IND                                                 EIGJ0940
   29 IF(AMAX)30,35,30                                                  EIGJ0950
   30 AT=AMAX                                                           EIGJ0960
      IND=2                                                             EIGJ0970
      F=.25                                                             EIGJ0980
      AMAX=F*AT                                                         EIGJ0990
      MAX=6                                                             EIGJ1000
      GO TO 5                                                           EIGJ1010
   31 MAX=MAX-1                                                         EIGJ1020
      IF(MAX)34,33,32                                                   EIGJ1030
   32 F=2.0*F*F                                                         EIGJ1040
      AMAX=AT*F                                                         EIGJ1050
   33 C=0.0                                                             EIGJ1060
      GO TO 5                                                           EIGJ1070
   34 IF(C)33,35,33                                                     EIGJ1080
*     COPY EIGENVALUES INTO E                                           EIGJ1090
   35 MM=1                                                              EIGJ1100
      NM=N&1                                                            EIGJ1110
      DO 36 M=1,N                                                       EIGJ1120
      E(M)=A(MM)                                                        EIGJ1130
   36 MM=MM&NM-M                                                        EIGJ1140
      RETURN                                                            EIGJ1150
      END
$      FORTRAN
$      INCODE  IBMF
      SUBROUTINE MTEVV(X,I1,I2,D,I3,I4,N,M)
      DIMENSION A(5050),X(I1,I2),D(I3,I4),E(100),NS(100),MS(100)
      K=1
      DO1J=1,N
      DO1I=J,N
      A(K)=X(I,J)
 1    K=K&1
      CALL EIGENJ(A,D,N,M,E,I3,I4)
C     THE FOLLOWING SYSTEM SORT PACKAGE ROUTINES SORT THE FLOATING POINT
C     ARRAY E OF N ELEMENTS IN DESCENDING ORDER.
      CALL SORTFL(-1)
      IF(M)4,3,4
 3    CALL SORT(E(1),E(2),N)
      DO 50 I=1,N
 50   X(I,I)=E(I)
      RETURN
 4    DO 5 I=1,N
       MS(I)=I
 5    NS(I)=I
C     THE FOLLOWING ROUTINE SORTS THE FLOATING POINT ARRAY E OF N
C     ELEMENTS IN DESCENDING ORDER AND ARRAY MS IS THE PASSIVE LIST
C     WHICH WILL BE REARRANGED IN THE SAME ORDER AS E.
      CALL SORT(E(1),E(2),N,MS(1))
C     THE FOLLOWING ROUTINES SORT THE FIXED POINT ARRAY MS OF N ELEMENTS
C     IN ASCENDING ORDER WITH ARRAY NS AS THE PASSIVE LIST.
       CALL SORTFX(1)
       CALL SORT(MS(1),MS(2),N,NS(1))
      CALL ARREV(N,NS,D,I3,I4)
      DO 6 I=1,N
 6    X(I,I)=E(I)
      RETURN
      END
$      FORTRAN
$      INCODE  IBMF
      SUBROUTINE ARREV(N,NS,D,I3,I4)
      DIMENSION NS(N),D(I3,I4)
      DO 10 I=1,N
 1    IF(NS(I)-I)2,10,2
 2    CALL XCH(N,I,NS(I),D,NS,I3,I4)
      GOTO1
 10   CONTINUE
      RETURN
      END
$      FORTRAN
$      INCODE  IBMF
      SUBROUTINE XCH(N,L,M,D,NS,I3,I4)
      DIMENSION        NS(N),D(I3,I4)
      I=L
       J=M
      DO 1 K=1,N
      T=D(K,I)
      D(K,I)=D(K,J)
 1    D(K,J)=T
      NS(I)=NS(J)
      NS(J)=J
      RETURN
      END