C indscals.f C This software implements the method called INDSCAL. 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 Unless you have some reason to do otherwise, it is recommended that you C start using SINDSCAL as your initial choice among the 4 programs C in this collection that carry out the INDSCAL method, namely, C INDSCAL, NINDSCAL, INDSCALS, SINDSCAL. 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 INDSCALS C (A SHORT VERSION OF INDSCAL) C INDIVIDUAL DIFFERENCES SCALING C BY C J. D. CARROLL AND J. J. CHANG C BELL TELEPHONE LABORATORIES C MURRAY HILL, NEW JERSEY C C THIS PROGRAM MATERIAL CONTAINS PROPRIETARY C INFORMATION OF BELL TELEPHONE LABORATORIES, C INCORPORATED AND IS NOT INTENDED FOR C PUBLICATION. THIS PRIVATE DOCUMENT IS THE C PROPERTY OF BELL LABORATORIES AND IS C LOANED TO THE RECIPIENT SOLELY TO FACILITATE C ENHANCEMENT OF THE ARTS TO WHICH THE C DOCUMENT PERTAINS. THE DOCUMENT IS TO BE C RETURNED TO BELL LABORATORIES UPON REQUEST. $ FORTREX C INDSCALS- SAME AS INDSCAL EXCEPT THE MAXIMUM WAYS IS 3. C 6/72 C BY CARROLL AND CHANG C PARAMETER DEFINITIONS C N - NUMBER OF MATRICES (MAX.=3). MATRIX 1 IS SUBJECTS WEIGHT C MATRIX, MATRIX 2 IS STIMULUS MATRIX AND MATRIX 3 CAN BE EITHER C STIMULUS OR CONDITION MATRIX. C MAXDIM- MAXIMUM NUMBER OF DIMENSIONS (MAX. =10) C MINDIM- MINIMUM NUMBER OF DIMENSIONS C IRDATA- 0, READ IN N WAY TABLES. C IRDATA- 1, READ IN LOWERHALF OF SIMILARITY MATRIX WITHOUT DIAGS. C IRDATA- 2, READ IN DISSIMILARITIES, LOWERHALF WITHOUT DIAGS. C IRDATA- 3, READ IN EUCLIDEAN DISTANCES, LOWERHALF WITHOUT DIAGS. C IRDATA- 4, READ IN LOWERHALF CORRELATION MATRIX WITHOUT DIAGS. C IRDATA- 5, READ IN LOWERHALF COVARIANCE MATRIX WITH DIAGS. C IRDATA- 6, READ IN FULL SYMMETRIC SIMILARITY MATRIX, DIAG. IGNORED C IRDATA- 7, READ IN FULL SYMMETRIC DISSIMILARITY MATRIX,DIAG. IGNORED C ITMAX - MAXIMUM NUMBER OF ITERATIONS. (MAX. =200) C ISET - 1, AT THE END OF ITERATIVE PROCEDURE, SET MATRIX 2 = C MATRIX 3, THEN GO BACK TO ITERATE AGAIN WHILE KEEPING C MATRICES 2 AND 3 CONSTANT C ISET - 0, DO NOT SET MATRIX 2 = MATRIX 3. C IOY - 0, COMPUTE NF DIMENSIONAL SOLUTION SIMULTANEOUSLY. C IOY - 1, COMPUTE EACH DIM. SOLUTION SUCCESSIVELY. (NOT AVAILABLE) C IDR - 1, COMPUTE CORRELATIONS BETWEEN SOLUTION AND DATA FOR EACH C SUBJECT. C IDR - 0, DO NOT COMPUTE CORRELATIONS FOR EACH SUBJECT. C ISAM - 1, KEEP MATRICES 2 AND 3 (STIMULI) UNCHANGED AND SOLVE FOR C THE REMAINING MATRICES. C ISAM - 0, ITERATE AND SOLVE FOR ALL MATRICES. C IPUNSP - 1, PUNCH ON CARDS SCALAR PRODUCT MATRICES. C IPUNSP - 0, DO NOT PUNCH ON CARDS SCALAR PRODUCT MATRICES. C IRN - 0, READ IN INITIAL CONFIGURATION. C IRN - A FIVE DIGIT ODD INTEGER FOR GENERATING A RANDOM STARTING C CONFIGURATION. C CRIT - CRITERION FOR TERMINATING THE ITERATIVE PROCEDURE. C VARIABLES C Y - DATA. Y IS A 3 SUBSCRIPTED VARIABLE. C Y(N1,N2,N3), WHERE N1,N2,N3 REFERS TO ELEMENTS C IN MATRICES 1 TO 3. (MAX. (N1*N2*N3)=18K). C (MAX. N(I)=100) C A - COORDINATE MATRIX. C A(I,J,K), WHERE I - ELEMENTS WITHIN MATRICES (MAX. =100) C J - DIMENSIONS. (MAX. =10) C K - MATRIX. (MAX. =3) C MAXIMUM ((MAX.I)*J)=700 C NWT - SIZE OF MATRICES (MAX.=100) C TIT - TITLE, USE NO MORE THAN ONE CARD. DIMENSION Y(18000),A(2100),AA(700),NWT(3),TIT(12),BNUM(100) 90 READ106,N,MAXDIM,MINDIM,IRDATA,ITMAX,ISET,IOY,IDR,ISAM,IPUNSP,IRN IF(N.EQ.0) GOTO999 N=3 IOY=0 C IN INDSCALS, N IS ALWAYS 3 AND IOY IS 0, BUT THEY ARE KEPT IN C ORDER TO BE THE SAME AS IN INDSCAL. NA=2 IF(ISAM.EQ.1) ISET=0 READ101,(TIT(I),I=1,12) READ100,(NWT(I),I=1,N) READ103,CRIT C FIND MAX. OF NWT 40 MNWT=NWT(1) IF(NWT(2).GT.MNWT) MNWT=NWT(2) IF(NWT(3).GT.MNWT) MNWT=NWT(3) NTAU=MNWT*MAXDIM IF(NTAU.LE.700) GOTO9 PRINT111 STOP C SET UP LABELS FOR USE ON MICROFILM. C PIBCI AND FORCNV ARE MICROFILM PACKAGE ROUTINES. 9 DO30I=1,MNWT 30 CALL PIBCI(I,BNUM(I),4H(I2)) NTAU=1 DO2I=1,N 2 NTAU=NTAU*NWT(I) IF(NTAU.LE.18000) GOTO6 PRINT110 STOP 6 NDIM=0 MMDIM=MAXDIM&MINDIM DO500IDIM=MINDIM,MAXDIM NF=MMDIM-IDIM MAXIT=ITMAX INORM=ISET ISTI=ISAM PRINT104 PRINT102,(TIT(I),I=1,12) PRINT107 PRINT108 PRINT105,N,NF,IRDATA,MAXIT,INORM,IOY,IDR,ISTI,IPUNSP,IRN,CRIT PRINT109,(NWT(I),I=1,N) PRINT107 4 IF(ISTI.NE.0) MAXIT=0 8 IF (IRN.NE.0) CALL RAND1(IRN) C INITIALIZE A AND NWT 5 N3=NWT(3) N2=NWT(2) N1=NWT(1) TAU=N1*N2*N3 IF(NDIM.EQ.0) GOTO44 NFO=NF&1 NQ=NFO*N CALL DIMA(A,MNWT,NQ,NFO,NF,N) GOTO42 44 IF(IRDATA)41,43,41 41 CALL RDATA(Y,N1,N2,N3, IRDATA,A,IPUNSP) GOTO42 43 NA=3 42 CALL CANDE(Y,N1,N2,N3,NA, A,MNWT,NF,N,NWT,IRDATA, ISTI,MAX 1IT,CRIT,IRN, INORM,NDIM,FVAF) IF(INORM.EQ.0) GOTO51 C SET MATRIX 2 TO MATRIX 3, ITERATE AGAIN TO COMPUTE SUBJECTS WEIGHTS. CALL SETA(A,MNWT,NF,N,N3) MAXIT=0 INORM=0 ISTI=2 GOTO42 C NORMALIZE A MATRICES. 51 CALL NORMA(A,AA,MNWT,NF,N,NWT,TIT,BNUM,NA,FVAF) C COMPUTE CORRELATIONS BETWEEN DATA AND SOLUTION FOR EACH N2 BY N3 C MATRIX. IF(IDR.NE.0) CALL SUBR(TIT,A,MNWT,NF,N,Y,N1,N2,N3) NDIM=1 500 CONTINUE GOTO90 100 FORMAT(18I4) 101 FORMAT(12A6) 102 FORMAT(1H012A6) 103 FORMAT(10F7.3) 104 FORMAT(65H1INDSCALS - A SHORT VERSION OF INDSCAL, ONLY UP TO 3 WAY 1 ANALYSIS) 105 FORMAT(I2,I3,I6,3I7, I5,3I7,3X,E12.4) 106 FORMAT(10I4,I6) 107 FORMAT(51H0**************************************************) 108 FORMAT(11H0PARAMETERS/65H0N DIM IRDATA ITMAX ISET IOY IDR 1ISAM IPUNSP IRN CRIT) 109 FORMAT(13H0MATRIX SIZES/2X,7I4) 110 FORMAT(50H0THE PRODUCT OF MATRIX SIZE EXCEEDS LIMIT OF 18000) 111 FORMAT(52H0(MAX. MATRIX SIZE * DIMENSION) EXCEEDS LIMIT OF 700) 999 CALL CLEAN CALL EXIT END $ FORTREX SUBROUTINE CANDE(Y,N1,N2,N3,NA, A,MNWT,ND,N,NWT,IRDATA,JST 1I ,NAXIT,CRIT,IRN, INORM,NDIM,FVAF) C PERFORMS CANONICAL DECOMPOSITION ANALYSIS DIMENSION Y(N1,N2,N3), A(MNWT,ND,3),S(100,10),R(10,10) DIMENSION ESQ(200),CORY(200),VAF(200),ER(200) DIMENSION NWT(3),SAY(10),FMT(12),SCR(201),NI(3) IZERO=0 MAXIT=NAXIT ISTI=JSTI NF=ND 54 IF(ISTI.EQ.2) GOTO303 IF(NDIM.EQ.1) GOTO302 46 IF(IRDATA)600,65,301 C READ IN DATA (Y MATRIX) 65 READ101,FMT DO2I1=1,N1 DO2I2=1,N2 2 READFMT,(Y(I1,I2,I3),I3=1,N3) C READ IN A MATRIX 301 IF(IRN)61,300,61 61 CALL ARAN(A,MNWT,NF,N,NWT) GOTO302 300 READ101,FMT DO4I=2,NA NN=NWT(I) DO4K=1,NF 4 READ FMT,(A(J,K,I),J=1,NN) IF(NA.EQ.3) GOTO302 DO304K=1,NF DO304J=1,N2 304 A(J,K,3)=A(J,K,2) 302 DO18I=1,3 18 NI(I)=NWT(I) 48 IT=0 PRINT1020 DO8I=2,NA NN=NWT(I) DO8J=1,NF 8 PRINT1002,J,(A(M,J,I),M=1,NN) IF(ISTI.EQ.1) GOTO303 GOTO40 36 IT=1 JJ1=ILAPSZ(0) C DO LOOP ON N WAY TABLE GOTO90 303 IT=0 90 DO100JB=1,N IF(ISTI.EQ.0) GOTO52 IF(JB.EQ.2.OR.JB.EQ.3) GOTO100 C COMPUTE R 52 DO12K1=1,NF DO12K2=1,K1 SUM1=1. DO3I=1,N IF(I.EQ.JB) GOTO3 MM=NWT(I) SUM=0. DO14J=1,MM 14 SUM=SUM&A(J,K1,I)*A(J,K2,I) SUM1=SUM*SUM1 3 CONTINUE R(K1,K2)=SUM1 12 R(K2,K1)=SUM1 68 MX=1 CALL MTINV(R,10,10,SCR,10,10,NF,MX) IF(MX.EQ.1)GOTO6 C PRINT MESSAGE PRINT1004 STOP 6 NI(JB)=1 C COMPUTE S NN=NWT(JB) DO9K=1,NF 9 A(1,K,JB)=1. K1=NI(1) K2=NI(2) K3=NI(3) DO50J=1,NN DO5K=1,NF 5 S(J,K)=0. DO10I3=1,K3 IF(JB.EQ.3) GOTO220 J3=I3 GOTO225 220 J3=J 225 DO10I2=1,K2 IF(JB.EQ.2) GOTO221 J2=I2 GOTO222 221 J2=J 222 DO16K=1,NF 16 SAY(K)=0. IF(JB.NE.1) GOTO210 AY=Y(J,J2,J3) DO211K=1,NF 211 SAY(K)=AY GOTO213 210 DO11I1=1,K1 YAL=Y(I1,J2,J3) DO11K=1,NF 11 SAY(K)=SAY(K)&A(I1,K,1)*YAL 213 DO13K=1,NF 13 S(J,K)=S(J,K)&A(I2,K,2)*A(I3,K,3)*SAY(K) 10 CONTINUE 50 CONTINUE C COMPUTE A DO15I=1,NN DO15J=1,NF SUM=0. DO17K=1,NF 17 SUM=SUM&S(I,K)*R(K,J) 15 A(I,J,JB)=SUM NI(JB)=NWT(JB) 100 CONTINUE C COMPUTE Y HAT AND ERROR 40 ERR=0. SY=0. SYH=0. SSQY=0. SSQYH=0. YY=0. DO200I3=1,N3 DO200I2=1,N2 DO201I1=1,N1 SUM=0. DO34I=1,NF 34 SUM=SUM&A(I1,I,1)*A(I2,I,2)*A(I3,I,3) YS=Y(I1,I2,I3) DIF=YS-SUM SY=SY&YS SYH=SYH&SUM SSQY=SSQY&YS**2 SSQYH=SSQYH&SUM**2 YY=YY&YS*SUM 201 ERR=ERR&DIF**2 200 CONTINUE CORYY=YY/SQRT(SSQY*SSQYH) RSQ=CORYY**2 EVA=1.-RSQ IF(IT.EQ.0)GOTO56 CORY(IT)=CORYY ESQ(IT)=EVA VAF(IT)=RSQ ER(IT)=ERR GOTO57 56 FCOR=CORYY FESQ=EVA FVAF=RSQ FER=ERR IF(ISTI-1) 57,69,67 69 PRINT1019 67 PRINT1023, FCOR,FVAF,FESQ,FER GOTO600 57 IF(IT-1)36,203,204 204 ETEST=PERR-ERR IF(ETEST.LE.CRIT) GOTO49 203 PERR=ERR IF(IT.GE.MAXIT)GOTO51 IT=IT&1 GOTO90 49 PRINT1008,IT,ETEST,CRIT GOTO64 51 PRINT1009,ETEST 64 JJ2=ILAPSZ(0) LAPSE=JJ2-JJ1 JT=LAPSE/64000 PRINT1003,LAPSE,JT 1003 FORMAT(7H0LAPSE=2I20) 66 PRINT1019 PRINT1022,IZERO,FCOR,FVAF,FESQ,FER PRINT1022,(I,CORY(I),VAF(I),ESQ(I),ER(I),I=1,IT) FVAF=VAF(IT) PUNCH 1005 DO226I=1,N NN=NWT(I) DO226J=1,NF 226 PUNCH1007,I,J,(A(M,J,I),M=1,NN) 600 RETURN 101 FORMAT(12A6) 1002 FORMAT(I3,10E12.4/(3X,10E12.4)) 1004 FORMAT(21H0R CANNOT BE INVERTED) 1005 FORMAT(23HUNNORMALIZED A MATRICES) 1007 FORMAT(2I1,5E13.5/(2X,5E13.5)) 1008 FORMAT (31H0REACHED CRITERION ON ITERATIONI3,3X,6HETEST=E12.4,3X, 15HCRIT=E12.4) 1009 FORMAT(25HREACHED MAXIMUM ITERATION,5X,6HETEST=E12.4) 1011 FORMAT(11H(2X,5E13.5)) 1012 FORMAT(2I1,5E13.5/(2X,5E13.5)) 1016 FORMAT(18I4) 1018 FORMAT(I7,3X,4E20.6) 1019 FORMAT(23H0HISTORY OF COMPUTATION//10H0ITERATION,4X,20HCORRELATION 1S BETWEEN,6X,3HVAF,12X,17HRESIDUAL VARIANCE/16X,16HY(DATA) AND YHA 2T,7X,6H(R**2)13X,8H(1-R**2),12X,11H(Y-YHAT)**2) 1020 FORMAT(24H0INITIAL STIMULUS MATRIX) 1021 FORMAT(24H0REACHED RUN TIME LIMITS) 1022 FORMAT(I7,3X,4E20.6) 1023 FORMAT(4X,5HFINAL,1X,4E20.6) END $ FORTRAN SUBROUTINE NORMA(A,AA,MNWT,K,N,NWT,TIT,BNUM,NA,FVAF) C PROGRAM TO NORMALIZE A MATRICES FROM THE CANNONICAL DECOMPOSITION C OF N WAY TABLE PROGRAM DIMENSION TIT(12),NWT(7),A(MNWT, K,N),S(10,10),BNUM(100),V(10) DIMENSION FTB(2),FATB(2) DIMENSION AA(MNWT,K),DIAG(10),LPAS(10) DATA FTB(1)/11H(5HTABLEI2)/ XN1=NWT(1) DO6I=2,N NN=NWT(I) DO6J=1,K SUM=0. DO5M=1,NN 5 SUM=SUM&A(M,J,I)**2 SUM=SQRT(SUM) DO3M=1,NN 3 A(M,J,I)=A(M,J,I)/SUM S(J,I)=SUM 6 CONTINUE NN=NWT(1) DO7J=1,K SS=1. DO8JJ=2,N SS=SS*S(J,JJ) 8 CONTINUE DO7M=1,NN 7 A(M,J,1)=A(M,J,1)*SS C COMPUTE SUM SQUARES OF WEIGHTS FOR EACH DIM. NN=NWT(1) DO11I=1,K LPAS(I)=I DIAG(I)=0. DO11J=1,NN 11 DIAG(I)=DIAG(I)&A(J,I,1)**2 CALL SORTFL(-1) CALL SORT(DIAG(1),DIAG(2),K,LPAS(1)) C PERMUTE DIMENSIONS ACCORDING TO SUM SQUARES IN DESCENDING ORDER DO15I=1,N NN=NWT(I) DO16J=1,K DO16M=1,NN 16 AA(M,J)=A(M,J,I) DO17J=1,K JJ=LPAS(J) DO17M=1,NN 17 A(M,J,I)=AA(M,JJ) 15 CONTINUE PRINT102 PUNCH102 PUNCH109 DO10I=1,NA NN=NWT(I) DO55J=1,K 55 PUNCH110,J,I,(A(M,J,I),M=1,NN) GOTO(41,42,43),I 41 PRINT106 GOTO34 42 PRINT111 GOTO34 43 PRINT112 34 DO10J=1,K 10 PRINT105,J,(A(M,J,I),M=1,NN) C COMPUTE SUMS OF PRODUCTS AND SUMS OF SQUARES OF A MATRIX DO25I=1,NA GOTO(45,46,47),I 45 PRINT101 GOTO48 46 PRINT103 GOTO48 47 PRINT104 48 NN=NWT(I) DO20L=1,K DO20J=1,K S(L,J)=0. DO20M=1,NN 20 S(L,J)= S(L,J)&A(M,L,I)*A(M,J,I) SUM=0. 23 DO22L=1,K IF(I.GT.1) GOTO22 V(L)=S(L,L) SUM=SUM&V(L) 22 PRINT105,L,( S(L,J),J=1,K) IF(I.GT.1) GOTO25 DO26II=1,K DO26JJ=1,II 26 S(II,JJ)=S(II,JJ)/SQRT(V(II)*V(JJ)) PRINT114 DO27L=1,K 27 PRINT105,L,(S(L,J),J=1,L) DO24II=1,K 24 V(II)=(V(II)*FVAF)/SUM PRINT113,(II,II=1,K) PRINT107, (V(II),II=1,K) 25 CONTINUE C PLOT TABLES DO35I=1,NA NN=NWT(I) CALL FORCNV(I,1,FATB,FTB,DUM1,DUM2) CALL PLOTX(A(1,1,I),MNWT, K,NN,TIT,FATB,2,BNUM) 35 CONTINUE 101 FORMAT(27H0SUM OF PRODUCTS (SUBJECTS)) 102 FORMAT(20H0NORMALIZED SOLUTION) 103 FORMAT(26H0SUM OF PRODUCTS (STIMULI)) 104 FORMAT(29H0SUM OF PRODUCTS (CONDITIONS)) 105 FORMAT(I3,3X,10E12.4/(6X,10E12.4)) 106 FORMAT(23H0SUBJECTS WEIGHT MATRIX) 107 FORMAT(4X,10E12.4) 109 FORMAT(11H(2X,5E13.5)) 110 FORMAT(2I1,5E13.5/(2X,5E13.5)) 111 FORMAT(16H0STIMULUS MATRIX) 112 FORMAT(17H0CONDITION MATRIX) 113 FORMAT(73H0APPROXIMATE PROPORTION OF TOTAL VARIANCE ACCOUNTED FOR 1BY EACH DIMENSION/4X,I7,9I12) 114 FORMAT(38H0NORMALIZED SUM OF PRODUCTS (SUBJECTS)) RETURN END $ FORTRAN SUBROUTINE RDATA(Y,N1,N2,N3, IRDATA,D,IPUNSP) DIMENSION Y(N1,N2,N3), D(N2,N3),FMT(12) C IRDATA- 1, READ IN LOWERHALF OF SIMILARITY MATRIX WITHOUT DIAGS. C IRDATA- 2, READ IN DISSIMILARITIES, LOWERHALF WITHOUT DIAGS. C IRDATA- 3, READ IN EUCLIDEAN DISTANCES, LOWERHALF WITHOUT DIAGS. C IRDATA- 4, READ IN LOWERHALF CORRELATION MATRIX WITHOUT DIAGS. C IRDATA- 5, READ IN LOWERHALF COVARIANCE MATRIX WITH DIAGS. C IRDATA- 6, READ IN FULL SYMMETRIC SIMILARITY MATRIX, DIAG. IGNORED C IRDATA- 7, READ IN FULL SYMMETRIC DISSIMILARITY MATRIX,DIAG. IGNORED C IADDC=1, COMPUTE ADDITIVE CONSTANT. C IADDC=0, DO NOT COMPUTE ADDITIVE CONSTANT. IF(IPUNSP.EQ.0) GOTO3 PUNCH 103 PUNCH 104 3 IADDC=1 IFIRST=2 II=N2 IF(IRDATA.EQ.3) IADDC=0 IF(IRDATA.GE.5) IFIRST=1 READ100,(FMT(I),I=1,12) 4 DO50I1=1,N1 DO40J=IFIRST,N2 IF(IRDATA-5)11,12,13 11 II=J-1 GOTO13 12 II=J 13 READFMT,(D(J,M),M=1,II) 40 CONTINUE DO41J=2,N2 JJ=J-1 DO41M=1,JJ IF(IRDATA.EQ.1.OR.IRDATA.EQ.6) GOTO14 GOTO15 14 D(J,M)=-D(J,M) 15 D(M,J)= D(J,M) 41 CONTINUE GOTO(43,43,43,47,46,43,43),IRDATA C COMPUTE ADDITIVE CONSTANT AND SCALAR PRODUCT MATRICES. 43 CALL TORGB(D,N2,IADDC) IF(IPUNSP.EQ.0) GOTO46 DO20I=1,N2 20 PUNCH 105,I1,I,(D(I,J),J=1,I) GOTO46 47 DO16I=1,N2 16 D(I,I)=1. 46 DO45I=1,N2 DO45J=1,N3 45 Y(I1,I,J) =D(I,J) 50 CONTINUE 100 FORMAT(12A6) 103 FORMAT(23HSCALAR PRODUCT MATRICES) 104 FORMAT(11H(6X,5E13.5)) 105 FORMAT(2I3,5E13.5/(6X,5E13.5)) RETURN END $ FORTRAN SUBROUTINE SUBR(TIT,A,MNWT,NF,N,Y,N1,N2,N3) DIMENSION TIT(12),A(MNWT,NF,N),Y(N1,N2,N3) XN=N2*N3 PRINT101 6 DO5I1=1,N1 SUMX=0. SUMY=0. SSQX=0. SSQY=0. SXY=0. DO3I2=1,N2 DO3I3=1,N3 YY=Y(I1,I2,I3) SUM=0. DO2I=1,NF 2 SUM=SUM&A(I1,I,1)*A(I2,I,2)*A(I3,I,3) SUMX=SUMX&SUM SUMY=SUMY&YY SSQX=SSQX&SUM**2 SSQY=SSQY&YY**2 3 SXY=SXY&SUM*YY R=(XN*SXY-SUMX*SUMY)/SQRT((XN*SSQX-SUMX**2)*(XN*SSQY-SUMY**2)) PRINT102,I1,R 5 CONTINUE 101 FORMAT (67H0CORRELATION BETWEEN COMPUTED SCORES AND ORIGINAL DATA 1FOR SUBJECTS/1H0) 102 FORMAT(I4,3X,F10.6) 103 FORMAT(9H0VARIABLE,4I4/9H SUBJECT,5X,1HR) RETURN END $ FORTRAN SUBROUTINE PLOTX(X,NL,NF,N,TIT,SUBT,NW,BNUM) DIMENSION FMXS(3),FXS(3) C PROGRAM TO PLOT N STIMULUS POINTS X IN NF DIMENSIONAL SPACE DIMENSION X(NL,NF),BNUM(1),TIT(1) DIMENSION FFM(2),FI(2),FJ(2),FDM(4),FFDM(3) DATA FMXS(1)/14H(6HSCALE=F7.2)/ 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) 31 DO260I=1,NNF JJ=I&1 DO260J=JJ,NF 33 CALL MAXI(X(1,I),X(1,J),XMAX,N) XSCALE=460./XMAX 34 CALL ROLL CALL REFRNC(-512,512,1024,1024) C PRINT TITLE CALL TEXT(-460,512,TIT,0,12) C PRINT SUBTITLE CALL TEXT1(-460,500,SUBT,0,NW) C LABEL NF DIMENSIONAL SOLUTION CALL TEXT(-460,460,FFDM,0,3) 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 TSP(IAX,IAY,1H*,1) C LABEL POINTS NAX=IAX-7 NAY=IAY&7 CALL TEXT1(NAX,NAY,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) C PRINT SCALE CALL FORCNV(XMAX,1,FXS,FMXS,DUM1,DUM2) CALL TEXT1(-460,-490,FXS,0,3) IF(NF-1)50,50,260 260 CONTINUE 50 RETURN END $ FORTRAN SUBROUTINE ARAN(A,MNWT,NF,N,NWT) DIMENSION A(MNWT,NF,N),NWT(N) MIN=0 MAX=1000 DO10I=2,N NN=NWT(I) DO2K=1,NF DO2J=1,NN CALL RAND(MIN,MAX,NUMB) 2 A(J,K,I)=(FLOAT(NUMB)-500.)*.001 10 CONTINUE RETURN END $ FORTRAN SUBROUTINE SETA(A,MNWT,NF,N,N3) DIMENSION A(MNWT,NF,N) DO2I=1,NF DO2J=1,N3 2 A(J,I,2)=A(J,I,3) RETURN END $ FORTRAN SUBROUTINE DIMA(A,MNWT,NQ,NFO,NF,N) DIMENSION A(MNWT,NQ) NN=NFO L=NF DO8M=2,N DO6I=1,NF L=L&1 NN=NN&1 DO6J=1,MNWT 6 A(J,L)=A(J,NN) NN=NN&1 8 CONTINUE 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 END $ GMAP * CALL RAND1(START) * CALL RAND(MIN,MAX,NUMB) * GENERATE RANDOM NUMBER AND RETURN IN NUMB SYMDEF RAND1 RAND1 LDQ 2,1* MPY MN ANQ =O377777777777 STQ RN TRA 0,1 SYMDEF RAND RAND LDQ 3,1* SBQ 2,1* ADQ 1,DL STQ T LDQ RN MPY MN ANQ =O377777777777 STQ RN MPY T LLS 1 ADA 2,1* STA 4,1* TRA 0,1 T BSS 1 RN OCT 532471462257 MN OCT 343277244615 END $ FORTRAN $ INCODE IBMF C TORGERSON METHOD TO COMPUTE SCALAR PRODUCT MATRICES SUBROUTINE TORGB(A,N,IADDC) DIMENSION DC(64),DR(64),A(N,N) XN=N DO8I=1,N DC(I)=0. DR(I)=0. 8 A(I,I)=0. SUMD=0. 51 IF(IADDC)27,27,26 26 ADDC=-100. N1=N-1 N2=N-2 DO30I=1,N2 JJ=I&1 DO30J=JJ,N1 KK=J&1 DO30K =KK,N DO31M=1,3 GOTO(21,22,23),M 21 CIJK=A(I,K)-A(J,K)-A(I,J) GOTO14 22 CIJK=A(J,K)-A(I,K)-A(I,J) GOTO14 23 CIJK=A(I,J)-A(J,K)-A(I,K) 14 IF(CIJK.GT.ADDC) ADDC=CIJK 31 CONTINUE 30 CONTINUE 27 DO40I=2,N II=I-1 DO40J=1,II IF(IADDC)29,40,29 29 A(I,J)=A(I,J)&ADDC A(J,I)=A(I,J) 40 SUMD=SUMD&A(I,J)**2 ADIJ=2.*SUMD/XN**2 DO7I=1,N DO7J=1,N DC(I)=DC(I)&A(J,I)**2 7 DR(I)=DR(I)&A(I,J)**2 C COMPUTE B - SCALAR PRODUCT MATRIX SUM=0. SUMD=0. DO10I=1,N DO10J=1,I QQ=(DR(I)&DC(J))/XN A (I,J)=(QQ-A(I,J)**2-ADIJ)*(.5) IF(I.EQ.J) GOTO11 SUM=SUM&A(I,J)**2 GOTO10 11 SUMD=SUMD&A(I,J)**2 10 CONTINUE C NORMALIZE SCALAR PRODUCT MATRIX Q=1./SQRT(SUMD&2.*SUM) DO12I=1,N DO12J=1,I A(I,J)=A(I,J)*Q 12 A(J,I)=A(I,J) 500 RETURN END $ FORTRAN $ INCODE IBMF CMTINV MTINV SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M) DIMENSION A(I1,I2), B(I3,I4) C C THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS- C JORDAN ELIMINATION SCHEME. ALL CALCULATIONS ARE DONE IN DOUBLE- C PRECISION ARITHMETIC C DOUBLE PRECISION B,BMAX,BMULT EQUIVALENCE (BMAX,BMULT) C C SINGLE TO DOUBLE PRECISION TRANSFER C 101 DO 103 I = 1,N 102 DO 103 J = 1,N 103 B(I,J) = A(I,J) C C FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMATRIX. C 200 DO 430 K = 1,N A(K,1) = FLOAT(K) A(K,2) = A(K,1) 203 BMAX = DABS(B(K,K)) 210 DO 219 I = K,N 211 DO 219 J = K,N IF (BMAX - DABS(B(I,J)) ) 213,219,219 213 BMAX = DABS ( B(I,J)) 214 A(K,1) = FLOAT(I) 215 A(K,2) = FLOAT(J) 219 CONTINUE 216 IF (BMAX - 1.D-38) 801,301,301 C C EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL. C 301 I = IFIX(A(K,1)) 302 J = IFIX(A(K,2)) 313 DO 314 M = 1,N BMULT = B(I,M) B(I,M) = B(K,M) 314 B(K,M) = BMULT 303 DO 304 M = 1,N BMULT = B(M,J) B(M,J) = B(M,K) 304 B(M,K) = BMULT C C ACTUAL INVERSION STEP. C BEGIN ROW ITERATION C 401 DO 425 I = 1,N C IGNORE ROW K. 402 IF (I - K) 405,425,405 C SET ROW MULTIPLIER. 405 BMULT = B(I,K) / B(K,K) C MODIFY ROW ELEMENTS. 410 DO 420 J = 1,N C IGNORE COLUMN K 411 IF (J - K) 414,412,414 412 B(I,J) = -BMULT 413 GO TO 420 414 B(I,J) = B(I,J) - B(K,J) * BMULT 420 CONTINUE 425 CONTINUE C C DIVIDE PIVOT ROW BY PIVOT ELEMENT C 426 BMULT = 1.D0/ B(K,K) 427 B(K,K) = 1.D0 C ACTUAL DIVISION 428 DO 429 J = 1,N 429 B(K,J) = B(K,J) * BMULT 430 CONTINUE C C REARRANGE AND REASSEMBLE INVERSE MATRIX. C 501 DO 512 IJK = 1,N 502 I = N-IJK&1 503 L = IFIX(A(I,2)) 504 J = IFIX(A(I,1)) 505 DO 508 M = 1,N 506 BMULT = B(I,M) 507 B(I,M) = B(L,M) 508 B(L,M) = BMULT 509 DO 512 M = 1,N 510 BMULT = B(M,I) 511 B(M,I) = B(M,J) 512 B(M,J) = BMULT C C MOVE A INVERSE BACK INTO A C 601 DO 603 I = 1,N 602 DO 603 J = 1,N 603 A(I,J) = B(I,J) 701 M = 1 702 RETURN C C ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38 C 801 M = 2 802 RETURN 999 END