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