C indscal.f C This software implements the method called INDSCAL. C The authors of this software are J D Carroll and Jih-Jie Chang; C it was adapted by Vithala R Rao and revised by Arun K Jain. 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 INSCAL-REVISED ON JUNE 28,1971 BY ARUN K. JAIN REVISED JAN. 1972 BY ARUN K. JAIN DIMENSION Y(18000),NWT(7),TIT(20),AA(2400),A(2400) COMMON NF,NQ WRITE(6,120) C 90 READ 106,N,MAXDIM,MINDIM,IRDATA,MAXIT,INORM,IOY, C XIDR,ISTI,IPUNSP,IRN,IVEC,IT 90 READ 106,N,MAXDIM,MINDIM,IRDATA,ITMAX,ISET,IOY,IDR,ISAM,IPUNSP, 1IRN,IVEC,IT IF(N.EQ.0) GOTO999 IF(IOY.EQ.1) IDR=0 IF(ISAM.EQ.1) ISET=0 READ101,(TIT(I),I=1,20) READ100,(NWT(I),I=1,N) READ103,CRIT IF(IT.EQ.0) IT=5 C FIND MAX. OF NWT 40 MNWT=NWT(1) DO 3 I=2,N IF(NWT(I).LE.MNWT)GOTO3 MNWT=NWT(I) 3 CONTINUE NTAU=MNWT*MAXDIM IF(NTAU.LE.2400) GOTO9 PRINT111 STOP 9 NTAU=1 DO 2 I=1,N 2 NTAU=NTAU*NWT(I) IF (NTAU.LE.18000) GO TO 6 PRINT110 STOP 6 NDIM=0 MMDIM=MAXDIM+MINDIM DO 500 IDIM=MINDIM,MAXDIM MAXIT=ITMAX INORM=ISET ISTI=ISAM NF=MMDIM-IDIM PRINT104,N,NF PRINT102,(TIT(I),I=1,20) PRINT107 PRINT108 PRINT105,N,NF,IRDATA,MAXIT,INORM,IOY,IDR,ISTI,IPUNSP,IRN,CRIT, -IVEC,IT PRINT109,(NWT(I),I=1,N) PRINT107 CALL IDENT 4 IF(ISTI.EQ.1.AND.N.EQ.3) MAXIT=0 C INITIALIZE A AND NWT IF(N.EQ.7)GOTO40 NONE=N+1 DO 35 I=NONE,7 35 NWT(I)=1 5 N7 = NWT(7) N6=NWT(6) N5=NWT(5) N4=NWT(4) N3=NWT(3) N2=NWT(2) N1=NWT(1) TAU=N1*N2*N3*N4*N5*N6*N7 IF(NDIM.EQ.0) GO TO 44 NFO=NF+1 NQ=NFO*N CALL DIMA(A,MNWT,NFO,N) GO TO 42 44 IF(IRDATA) 41,42,41 41 CALL RDATA(Y,N1,N2,N3,N4,N5,N6,N7,IRDATA,A,IPUNSP,IVEC,IT) 42 CALL CANDE(Y,N1,N2,N3,N4,N5,N6,N7,A,MNWT,NF,N,NWT,IRDATA, ISTI,MAX 1IT,CRIT,IRN,IOY,INORM,IT,NDIM,TAU) IF(IOY.EQ.1.OR.INORM.EQ.0) GOTO51 C SET MATRIX 2 TO MATRIX 3, ITERATE AGAIN TO COMPUTE THE REMAINING C MATRICES PRINT107 PRINT112 CALL SETA(A,MNWT,NF,N,N3) MAXIT=0 INORM=0 ISTI=2 IF(N.GT.3) MAXIT=5 GOTO42 C NORMALIZE A MATRICES. 51 CALL NORMA(A,AA,MNWT,NF,N,NWT,TIT) 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,N4,N5,N6,N7) NDIM=1 500 CONTINUE GOTO90 100 FORMAT(18I4) 101 FORMAT(20A4) 102 FORMAT(1H020A4) 103 FORMAT(10F7.3) 104 FORMAT(73H1INDIFF- INDIVIDUAL DIFFERENCES ANALYSIS USING CANONICAL 1 DECOMPOSITION OFI3,13H WAY TABLE INI3,11H DIMENSIONS) 105 FORMAT(I2,I3,I6,3I7,I5,3I7,3X,F7.3,2I4) 106 FORMAT(10I4,I6,2I4) 107 FORMAT(51H0**************************************************) 108 FORMAT(11H0PARAMETERS/77HON NF IRDATA MAXIT INORM IOY IDR 1ISTI IPUNSP IRN CRIT IVEC IT) 109 FORMAT(13H0MATRIX SIZES/2X,7I4) 110 FORMAT(50H0THE PRODUCT OF MATRIX SIZE EXCEEDS LIMIT OF 6) 0000 111 FORMAT(52H0(MAX. MATRIX SIZE * DIMENSION) EXCEEDS LIMIT OF2400) 112 FORMAT(44H0EQUATE MATRIX 2 AND MATRIX 3, ITERATE AGAIN) 120 FORMAT('1',2X,'*************************************************** 1********************'//2X,'** INDSCAL 2 ***',//2X,'** PROGRAM ORIGINALLY WRITTEN BY J.DO 3UGLAS CARROLL AND JIE JIH CHANG OF BELL LABS. , NEW JERSEY **'/ 42X,'** ADAPTED TO SYSTEM IBM 360/75 BY VITHALA R. RAO , UNIV. OF 5 PENNSYLVANIA , AUGUST 1969 **'/ 2X,'**************************** 7*********************************************'/) 999 CALL EXIT END SUBROUTINE CANDE(Y,N1,N2,N3,N4,N5,N6,N7,A,MNWT,ND,N,NWT,IRDATA,JST 1I ,NAXIT,CRIT,IRN,IOY,INORM,NTAPE,NDIM,TAU) C PERFORMS CANONICAL DECOMPOSITION ANALYSIS C LIMITS 5,4,10000 0000 REAL*8 SCR(201) DIMENSION Y(N1,N2,N3,N4,N5,N6,N7),A(MNWT,ND,7),S(100,10),R(10,10) DIMENSION Q5(10),Q4(10),Q3(10),Q2(10),ESQ(50),CORY(50) DIMENSION NI(7) DIMENSION NWT(7),SAY(10),FMT(18) DIMENSION TRV(50),VAF(50) C ITZ AND IRSD ARE CONTROL PARAMETERS FOR IOY = 1 ITZ=0 IRSD=0 IZERO=0 REV=1. MAXIT=NAXIT ISTI=JSTI NF=ND L=1 IF(N.EQ.7) GOTO54 NONE=N+1 DO 35 I=NONE,7 DO35J=1,NF DO35M=1,MNWT 35 A(M,J,I)=1. 54 IF(ISTI.GT.1.OR.NDIM.EQ.1) GO TO 48 46 IF(IRDATA)600,65,301 C READ IN DATA (Y MATRIX) 65 READ101,(FMT(I),I=1,18) DO2I7=1,N7 DO2I6=1,N6 DO2I5=1,N5 DO2I4=1,N4 DO2I1=1,N1 DO 2 I2=1,N2 2 READ(NTAPE,FMT)(Y(I1,I2,I3,I4,I5,I6,I7),I3=1,N3) C READ IN A MATRIX 301 IF(IRN)61,300,61 61 CALL ARAN(A,MNWT,NF,N,NWT,IRN) GOTO47 300 READ101,(FMT(I),I=1,18) DO4I=2,N NN=NWT(I) DO 4 J=1,NN 4 READ FMT,(A(J,K,I),K=1,NF) 47 DO53J=1,NF DO53M=1,MNWT 53 A(M,J,1)=1. 48 IT=0 PRINT1020 DO 8 I=1,N PRINT1001,I NN=NWT(I) DO8J=1,NF 8 PRINT1002,J,(A(M,J,I),M=1,NN) GOTO40 55 IF(IOY.EQ.0) GOTO302 IF(INORM.EQ.0) IRSD=1 CRIT = -100.0 NFF=NF NF=1 303 PRINT1017,NF 302 DO18I=1,7 18 NI(I)=NWT(I) 36 IT=1 90 PRINT1010,IT C DO LOOP ON N WAY TABLE 62 DO100JB=1,N IF(ISTI.EQ.0) GOTO52 IF(JB.EQ.2.OR.JB.EQ.3) GOTO100 C COMPUTE R 52 DO 12 K1=L,NF DO12K2=L,K1 SUM1=1. DO3I=1,N IF(I.EQ.JB) GOTO3 MM=NWT(I) SUM = 0.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 IF(IOY.EQ.0) GOTO68 R(L,L) = 1.0/R(L,L) GOTO6 68 MX=1 CALL MTINV(R,10,10,SCR,10,10,NF,MX) IF(MX.EQ.1)GOTO6 C PRINT MESSAGE PRINT1004 STOP 6 PRINT1001,JB NI(JB)=1 C COMPUTE S NN=NWT(JB) DO9I=1,NN DO9K=L,NF 9 A(I,K,JB)=1. DO50J=1,NN DO5K=L,NF 5 S(J,K)=0. K7=NI(7) DO10I7=1,K7 J7=INDY(JB,7,I7,J) K6=NI(6) DO10I6=1,K6 J6=INDY(JB,6,I6,J) K5=NI(5) DO10I5=1,K5 DO20K=L,NF 20 Q5(K)=A(I7,K,7)*A(I6,K,6)*A(I5,K,5) J5=INDY(JB,5,I5,J) K4=NI(4) DO10I4=1,K4 DO21K=L,NF 21 Q4(K)=Q5(K)*A(I4,K,4) J4=INDY(JB,4,I4,J) K3=NI(3) DO10I3=1,K3 2001 DO22K=L,NF 22 Q3(K)=Q4(K)*A(I3,K,3) J3 = INDY(JB,3,I3,J) K2=NI(2) DO10I2=1,K2 DO23K=L,NF 23 Q2(K)=Q3(K)*A(I2,K,2) J2=INDY(JB,2,I2,J) K1=NI(1) DO16K=L,NF 16 SAY(K)=0. DO11I1=1,K1 J1=INDY(JB,1,I1,J) DO11K=L,NF 11 SAY(K)=SAY(K)+A(I1,K,1)*Y(J1,J2,J3,J4,J5,J6,J7) DO13K=L,NF 13 S(J,K)=S(J,K)+Q2(K)*SAY(K) 10 CONTINUE 50 CONTINUE C COMPUTE A DO15I=1,NN DO15J=L,NF A(I,J,JB)=0. DO15K=L,NF 15 A(I,J,JB)=A(I,J,JB)+S(I,K)*R(K,J) DO7I=L,NF 7 PRINT1002,I,(A(J,I,JB),J=1,NN) 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. DO200I7=1,N7 DO200I6=1,N6 DO200I5=1,N5 DO30I=L,NF 30 Q5(I)=A(I7,I,7)*A(I6,I,6)*A(I5,I,5) DO200I4=1,N4 DO31I=L,NF 31 Q4(I)=Q5(I)*A(I4,I,4) DO200I3=1,N3 DO32I=L,NF 32 Q3(I)=Q4(I)*A(I3,I,3) DO200I2=1,N2 DO33I=L,NF 33 Q2(I)=Q3(I)*A(I2,I,2) DO201I1=1,N1 SUM=0. DO34I=L,NF 34 SUM=SUM+Q2(I)*A(I1,I,1) YS=Y(I1,I2,I3,I4,I5,I6,I7) DIF=YS-SUM SY=SY+YS SYH=SYH+SUM SSQY=SSQY+YS**2 SSQYH=SSQYH+SUM**2 YY=YY+YS*SUM IF(IRSD.EQ.1.AND.IT.GE.MAXIT) Y(I1,I2,I3,I4,I5,I6,I7)=DIF 201 ERR=ERR+DIF**2 200 CONTINUE YMEAN=SY/TAU YHMEAN=SYH/TAU CORYY=YY/SQRT(SSQY*SSQYH) RSQ=CORYY**2 EVA=1.-CORYY**2 XTRV=EVA*REV IF(IOY.EQ.1) RSQ=1.0-XTRV IF(IT.EQ.0)GOTO56 CORY(IT)=CORYY ESQ(IT)=EVA TRV(IT) = XTRV VAF(IT) = RSQ PRINT1007,ERR GOTO57 56 FCOR=CORYY FESQ=EVA FTRV=XTRV FVAF=RSQ IF(ITZ.EQ.1.AND.IT.EQ.0) GO TO 303 57 IF(IT-1)55,203,204 204 ETEST=PERR-ERR IF(ETEST.LE.CRIT) GO TO 49 203 PERR=ERR IF(IT.GE.MAXIT)GOTO51 IT=IT+1 GOTO90 49 PRINT1008,ETEST,CRIT GOTO502 51 PRINT1009,ETEST 502 PRINT1006 DO17I=1,N PRINT1001,I NN=NWT(I) DO17J=L,NF 17 PRINT1002,J,(A(M,J,I),M=1,NN) IF(IOY.EQ.0) GO TO 64 IF(IRSD.EQ.1) GO TO 63 IF(IOY.EQ.1.AND.INORM.EQ.1) GOTO59 GO TO 63 59 NN=NWT(2) DO58J=L,NF DO58I=1,NN 58 A(I,J,2)=A(I,J,3) MAXIT=0 IRSD=1 ISTI=2 IF(N.GT.3) MAXIT=5 GOTO62 63 PRINT 1021 PRINT1018 ,IZERO,FCOR,FVAF,FESQ,FTRV PRINT 1018,(I,CORY(I),VAF(I),ESQ(I),TRV(I),I=1,IT) GO TO 501 64 PRINT1019 PRINT1022 ,IZERO,FCOR,FVAF,FESQ PRINT 1022, (I,CORY(I),VAF(I),ESQ(I),I=1,IT) 501 IF(IOY.EQ.0.OR.NF.EQ.NFF) GOTO600 NF=NF+1 L=NF IRSD=0 MAXIT=NAXIT ISTI=JSTI REV=REV*EVA IT=0 ITZ=1 GO TO 40 600 RETURN 101 FORMAT(18A4) 1001 FORMAT(7H0MATRIXI2) 1002 FORMAT(I3,10F12.4/(3X,10F12.4)) 1004 FORMAT(21H0R CANNOT BE INVERTED) 1006 FORMAT(11H0A MATRICES) 1007 FORMAT(7H0ERROR=F13.5) 1008 FORMAT(18H0REACHED CRITERION5X,6HETEST=F12.8,3X,5HCRIT=F12.8) 1009 FORMAT(25HREACHED MAXIMUM ITERATION,5X,6HETEST=F12.8) 1010 FORMAT(11H0ITERATION=I3) 1011 FORMAT(11H(2X,5F13.5)) 1012 FORMAT(2I1,5F13.5/(2X,5F13.5)) 1016 FORMAT(18I4) 1017 FORMAT(18H0COMPUTE DIMENSION,I4) 1018 FORMAT(I7,3X,4F20.6) 1019 FORMAT(23H1HISTORY OF COMPUTATION//10H0ITERATION,4X,20HCORRELATION 1S BETWEEN,6X,28HNORMALIZED RESIDUAL VARIANCE/16X,16HY(DATA) AND YH 1AT,8X,8H(1-R**2)) 1020 FORMAT(19H0INITIAL A MATRICES) 1021 FORMAT(23H1HISTORY OF COMPUTATION//10H0ITERATION,4X,20HCORRELATION 1S BETWEEN,6X,3HVAF,12X,17HRESIDUAL VARIANCE,4X,19HTOTAL RES. VARIA 2NCE/16X,16HY(DATA) AND YHAT,6X,7H(1-TRV),13X,8H(1-R**2),15X,5H(TRV 3)) 1022 FORMAT(I7,3X,3F20.6) END SUBROUTINE RDATA(Y,N1,N2,N3,N4,N5,N6,N7,IRDATA,D,IPUNSP,IVEC,IT) DIMENSION Y(N1,N2,N3,N4,N5,N6,N7) ,D(N2,N3),FMT(18) DIMENSION V(400) 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 IF IVEC=1 THE PROGRAM WILL READ TRICON TYPE OUTPUT C IF A MATRIX IS READ IN SET IVEC = 0 C IADDC=1, COMPUTE ADDITIVE CONSTANT. C IADDC=0, DO NOT COMPUTE ADDITIVE CONSTANT. NUMB = N2*(N2-1)/2 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,18) 100 FORMAT(18A4) DO50I7=1,N7 DO50I6=1,N6 DO50I5=1,N5 DO50I4=1,N4 4 DO50I1=1,N1 IF(IVEC.NE.1) GO TO 62 60 READ(IT,FMT)(V(K),K=1,NUMB) LS = 0 MN = N2-1 DO 61 M=1,MN MM = M+1 DO 61 J=MM,N2 LS = LS + 1 61 D(J,M) = V(LS) GO TO 63 62 CONTINUE DO40J=IFIRST,N2 IF(IRDATA-5) 11,12,13 11 II=J-1 GOTO13 12 II=J 13 READ(IT,FMT)(D(J,M),M=1,II) 40 CONTINUE 63 CONTINUE IF(IRDATA.EQ.7) GOTO43 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,I4,I5,I6,I7)=D(I,J) 50 CONTINUE 101 FORMAT(8H0SUBJECTI4) 102 FORMAT(5H0DATA) 103 FORMAT(23HSCALAR PRODUCT MATRICES) 104 FORMAT(11H(6X,5F13.5)) 105 FORMAT(2I3,5F13.5/(6X,5F13.5)) RETURN END SUBROUTINE NORMA(A,AA,MNWT,K,N,NWT,TIT) C PROGRAM TO NORMALIZE A MATRICES FROM THE CANNONICAL DECOMPOSITION C OF N WAY TABLE PROGRAM DIMENSION TIT(20),NWT(7),A(MNWT, K,N),S(10,10),BNUM(100) DIMENSION Z(100,10) DIMENSION AA(MNWT,K),DIAG(10),LPAS(10) DIMENSION PAS(10) DO 6 I=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 OF SQUARES OF WEIGHTS FOR EACH DIMENSION NN=NWT(1) DO 11 I=1,K PAS(I) = I DIAG(I)=0. DO 11 J=1,NN 11 DIAG(I)=DIAG(I)+A(J,I,1)**2 CALL SORT(DIAG,K,PAS,BNUM,BNUM,1,-1) C PERMUTE DIMENSIONS ACCORDING TO SUM OF SQUARES OF WEIGHTS IN C DESCENDING ORDER DO 15 I=1,N NN=NWT(I) DO 16 J=1,K DO 16 M=1,NN 16 AA(M,J)=A(M,J,I) DO 17 J=1,K JJ=PAS(J) DO 17 M=1,NN 17 A(M,J,I)=AA(M,JJ) 15 CONTINUE PRINT 103,(TIT(I),I=1,20) PRINT102 PUNCH102 PUNCH109 DO10I=1,N NN=NWT(I) DO 55 M=1,NN 55 PUNCH110,M,I,(A(M,J,I),J=1,K) 33 PRINT106,I 34 DO 10 M=1,NN 10 PRINT 105,M,(A(M,J,I),J=1,K) C COMPUTE SUMS OF PRODUCTS AND SUMS OF SQUARES OF A MATRIX DO25I=1,N PRINT106,I SSQ=0. 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) DO21L=1,K DO21M=1,NN 21 SSQ=SSQ+A(M,L,I)**2 PRINT107 DO22L=1,K 22 PRINT105,L,( S(L,J),J=1,K) PRINT108,SSQ 25 CONTINUE C PLOT TABLES DO35I=1,N NN=NWT(I) DO 40 J=1,K DO 40 M=1,NN 40 Z(M,J)=A(M,J,I) CALL PLOTX(Z,K,NN,I) 35 CONTINUE 102 FORMAT(22H0NORMALIZED A MATRICES) 103 FORMAT(1H120A4) 105 FORMAT(I3,3X,10F12.5/(6X,10F12.5)) 106 FORMAT(7H0MATRIXI3) 108 FORMAT(16H0SUM OF SQUARES=F13.5) 107 FORMAT(17H0SUMS OF PRODUCTS) 109 FORMAT(11H(2X,5F13.5)) 110 FORMAT(2I3,5F13.5/(6X,5F13.5)) RETURN END SUBROUTINE SUBR(TIT,A,MNWT,NF,N,Y,N1,N2,N3,N4,N5,N6,N7) DIMENSION TIT(20),A(MNWT,NF,N),Y(N1,N2,N3,N4,N5,N6,N7) SR=0.0 SRR=0.0 XN=N2*N3 PRINT100,(TIT(I),I=1,20) PRINT101 DO5I7=1,N7 DO5I6=1,N6 DO5I5=1,N5 DO5I4=1,N4 IF(N.LE.3)GOTO6 GOTO(11,12,13,14),N 11 PRINT103,I4 GOTO6 12 PRINT103,I4,I5 GOTO6 13 PRINT103,I4,I5,I6 GOTO6 14 PRINT103,I4,I5,I6,I7 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,I4,I5,I6,I7) SUM=0. DO2I=1,NF 2 SUM=SUM+A(I1,I,1)*A(I2,I,2)*A(I3,I,3)*A(I4,I,4)*A(I5,I,5)*A(I6,I,6 1)*A(I7,I,7) 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 SR=SR+R SRR=SRR+R*R 5 CONTINUE XN1 = N1 SR=SR/XN1 SRR=SRR/XN1 WRITE(6,104) SR,SRR 100 FORMAT(1H120A4) 101 FORMAT (67H0CORRELATION BETWEEN COMPUTED SCORES AND ORIGIONAL DATA 1FOR SUBJECTS/1H0) 102 FORMAT(I4,3X,F10.6) 103 FORMAT(9H0VARIABLE,4I4/9H SUBJECT,5X,1HR) 104 FORMAT(//10X,'AVERAGE SUBJECT CORR. COEFT.',F8.5/10X, 'MEAN SQ XUARE CORR. COEFT.' ,F9.5//) RETURN END SUBROUTINE ARAN(A,MNWT,NF,N,NWT,IRN) DIMENSION A(MNWT,NF,N),NWT(N) DO10I=2,N NN=NWT(I) DO2K=1,NF DO2J=1,NN CALL RANDU(IRN,IY,X) IRN=IY 2 A(J,K,I)=X-.5 10 CONTINUE RETURN END SUBROUTINE PLOTX(V,K,N,I) DIMENSION V(100,10),X(100),Y(100) KK=K-1 DO 10 IM=1,KK KM=IM+1 DO 10 JM=KM,K WRITE (6,20) IM,JM,I 20 FORMAT('1',10X,'THIS IS PLOT OF DIMENSION ',I2 ,'VS.DIMENSION ',I2 X,'FOR TABLE NO. ',I2/) DO 15 M=1,N X(M) = V(M,IM) 15 Y(M) = V(M,JM) CALL GRAPH(+1.2,-1.2,+1.2,-1.2) CALL AXES CALL LOC2(X,Y,N,1) CALL PLOT 10 CONTINUE RETURN END SUBROUTINE GRAPH(XMX,XMN,YMX,YMN) GRPH 10 DIMENSION Z(1),X(1),Y(1),SMALL(13),CHAR(50),HOLL(11) GRPH 20 REAL ITEM(55,121) GRPH 30 DATA CHAR /'1','2','3','4','5','6','7','8','9','A','B','C', GRPH 50 1 'D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S', GRPH 60 2 'T','U','V','W','X','Y','Z','+','/','=','*','&','$','@','%','?', GRPH 70 3 '<','(',')','"',';',' '/ GRPH 80 DATA HOLL /' ','X','2','3','4','5','6','7','8','9','M'/ GRPH 90 DATA BLANK,DD,PLUS,AIE,AMINUS,ORIG,AMULT,GT50/' ','.','+','\', GRPH 100 1 '-','0','#','>'/ GRPH 110 DO 115 I = 1,55 GRPH 120 DO 115 J = 1,121 GRPH 130 115 ITEM(I,J) = BLANK GRPH 140 XMAX = XMX GRPH 150 XMIN = XMN GRPH 160 YMAX = YMX GRPH 170 YMIN = YMN GRPH 180 RETURN GRPH 190 ENTRY LIMIT(Z,NP,IS) GRPH 200 ZMAX = -1.0E9 GRPH 210 ZMIN = +1.0E9 GRPH 220 DO 10 I = 1,NP GRPH 230 IF(Z(I) .GT. ZMAX) ZMAX = Z(I) GRPH 240 10 IF(Z(I) .LT. ZMIN) ZMIN = Z(I) GRPH 250 IF(IS .EQ. 2) GO TO 11 GRPH 260 IF(ZMAX .GT. XMAX) XMAX = ZMAX GRPH 270 IF(ZMIN .LT. XMIN) XMIN = ZMIN GRPH 280 GO TO 12 GRPH 290 11 IF(ZMAX .GT. YMAX) YMAX = ZMAX GRPH 300 IF(ZMIN .LT. YMIN) YMIN = ZMIN GRPH 310 12 RETURN GRPH 320 ENTRY SCALE(IP) GRPH 330 XEXP = 0.05 * (XMAX - XMIN) GRPH 340 XMAX = XMAX + XEXP GRPH 350 XMIN = XMIN - XEXP GRPH 360 YEXP = 0.05 * (YMAX - YMIN) GRPH 370 YMAX = YMAX + YEXP GRPH 380 YMIN = YMIN - YEXP GRPH 390 IF(IP .NE. 1) RETURN GRPH 400 IF(YMAX - XMAX) 20,21,22 GRPH 410 20 YMAX = XMAX GRPH 420 GO TO 21 GRPH 430 22 XMAX = YMAX GRPH 440 21 CONTINUE GRPH 450 IF(YMIN - XMIN) 23,24,25 GRPH 460 23 XMIN = YMIN GRPH 470 GO TO 24 GRPH 480 25 YMIN = XMIN GRPH 490 24 EXP = 0.1667 * (XMAX - XMIN) GRPH 500 XMAX = XMAX + EXP GRPH 510 XMIN = XMIN - EXP GRPH 520 RETURN GRPH 530 ENTRY AXES GRPH 540 DELX = (XMAX - XMIN)/120.0 GRPH 550 DELY = (YMAX - YMIN)/54.0 GRPH 560 K = 0 GRPH 570 M = 0 GRPH 580 IF(XMIN .LT. 0.0 .AND. XMAX .GT. 0.0) GO TO 30 GRPH 590 GO TO 31 GRPH 600 30 K = (-XMIN/DELX) + 1.5 GRPH 610 DO 32 I = 1,55 GRPH 620 32 ITEM(I,K) = AIE GRPH 630 31 CONTINUE GRPH 640 IF(YMIN .LT. 0.0 .AND. YMAX .GT. 0.0) GO TO 33 GRPH 650 GO TO 34 GRPH 660 33 M = (YMAX/DELY) + 1.5 GRPH 670 DO 35 J = 1,121 GRPH 680 35 ITEM(M,J) = AMINUS GRPH 690 IF(K .EQ. 0 .OR. M .EQ. 0) GO TO 34 GRPH 700 ITEM(M,K) = ORIG GRPH 710 34 RETURN GRPH 720 ENTRY LOC1(X,Y,NPOI) GRPH 730 DELX = (XMAX - XMIN)/120.0 GRPH 740 DELY = (YMAX - YMIN)/54.0 GRPH 750 DO 50 II = 1,NPOI GRPH 760 I = (YMAX - Y(II))/DELY + 1.5 GRPH 770 J = (X(II) - XMIN)/DELX + 1.5 GRPH 780 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 50 GRPH 790 DO 52 JJ = 1,10 GRPH 800 IF(ITEM(I,J) .EQ. HOLL(JJ)) GO TO 51 GRPH 810 52 CONTINUE GRPH 820 IF(ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS.OR.ITEM(I,J).EQ.ORIG) GRPH 830 1 ITEM(I,J) = HOLL(2) GRPH 840 GO TO 50 GRPH 850 51 ITEM(I,J) = HOLL(JJ+1) GRPH 860 50 CONTINUE GRPH 870 RETURN GRPH 880 ENTRY LOC2(X,Y,NPOI,ISTART) GRPH 890 DELX = (XMAX - XMIN)/120.0 GRPH 900 DELY = (YMAX - YMIN)/54.0 GRPH 910 ISPN = ISTART + NPOI - 1 GRPH 920 DO 53 II = ISTART,ISPN GRPH 930 I = (YMAX - Y(II))/DELY + 1.5 GRPH 940 J = (X(II) - XMIN)/DELX + 1.5 GRPH 950 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 53 GRPH 960 IF(ITEM(I,J).EQ.BLANK.OR.ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS GRPH 970 1 .OR.ITEM(I,J).EQ.ORIG) GO TO 54 GRPH 980 ITEM(I,J) = AMULT GRPH 990 GO TO 53 GRPH1000 54 CONTINUE GRPH1010 IF(II .GT. 50) GO TO 55 GRPH1020 ITEM(I,J) = CHAR(II) GRPH1030 GO TO 53 GRPH1040 55 ITEM(I,J) = GT50 GRPH1050 53 CONTINUE GRPH1060 RETURN GRPH1070 ENTRY LOC3(X,Y,NPOI,POINT) GRPH1080 DELX = (XMAX - XMIN)/120.0 GRPH1090 DELY = (YMAX - YMIN)/54.0 GRPH1100 DO 56 II = 1,NPOI GRPH1110 I = (YMAX - Y(II))/DELY + 1.5 GRPH1120 J = (X(II) - XMIN)/DELX + 1.5 GRPH1130 IF(I.GT.55.OR.I.LT.1.OR.J.GT.121.OR.J.LT.1) GO TO 56 GRPH1140 IF(ITEM(I,J).EQ.BLANK.OR.ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS GRPH1150 1 .OR. ITEM(I,J) .EQ. ORIG .OR. ITEM(I,J) .EQ. POINT) GO TO 57 GRPH1160 ITEM(I,J) = AMULT GRPH1170 GO TO 56 GRPH1180 57 ITEM(I,J) = POINT GRPH1190 56 CONTINUE GRPH1200 RETURN GRPH1210 ENTRY PLOT GRPH1220 SMALL(1) = XMIN GRPH1230 DO 120 I = 1,12 GRPH1240 120 SMALL(I+1) = SMALL(I) + 10.0 * DELX GRPH1250 VALUE = YMAX + DELY GRPH1260 WRITE(6,301) GRPH1270 301 FORMAT(9X,' +.........+.........+.........+.........+.........+...GRPH1280 1......+.........+.........+.........+.........+.........+.........GRPH1290 2+') GRPH1300 DO 330 I = 1,55 GRPH1310 VALUE = VALUE - DELY GRPH1320 II = I - 1 GRPH1330 L = II + 6 GRPH1340 IF(L/6*6 - L) 320,310,320 GRPH1350 310 A = PLUS GRPH1360 WRITE(6,302) VALUE,A,(ITEM(I,J), J=1,121),A GRPH1370 302 FORMAT(1X,F8.3,123A1) GRPH1380 GO TO 330 GRPH1390 320 A = DD GRPH1400 WRITE(6,303) A,(ITEM(I,J), J=1,121),A GRPH1410 303 FORMAT(9X,123A1) GRPH1420 330 CONTINUE GRPH1430 WRITE(6,301) GRPH1440 WRITE(6,304) (SMALL(K), K=1,12) GRPH1450 304 FORMAT(4X,12F10.3) GRPH1460 WRITE(6,305) SMALL(13) GRPH1470 305 FORMAT(1H+,123X,F8.2/) GRPH1480 RETURN GRPH1490 ENTRY IDENT GRPH1500 WRITE(6,400) (J, J=1,25) GRPH1510 400 FORMAT(/35X,'*****IDENTIFICATION KEY FOR PLOTS WITH IDENTIFIED POIGRPH1520 1NTS*****'//' PT #',25I5) GRPH1530 WRITE(6,401) (CHAR(I), I=1,25) GRPH1540 401 FORMAT(' CHAR',25(4X,A1)/) GRPH1550 WRITE(6,402) (J, J=26,50) GRPH1560 402 FORMAT(' PT #', 25I5) GRPH1570 WRITE(6,403) (CHAR(I), I=26,50) GRPH1580 403 FORMAT(' CHAR', 25(4X,A1)//) GRPH1590 WRITE(6,404) GT50, AMULT GRPH1600 404 FORMAT(' POINT NUMBERS ABOVE 50 IDENTIFIED AS ',A1/' MULTIPLE POIGRPH1610 1NTS IDENTIFIED AS ',A1/) GRPH1620 RETURN GRPH1630 END GRPH1640 FUNCTION INDY(JB,M,I,J) IF(JB.EQ.M) GOTO20 INDY=I GOTO5 20 INDY=J 5 RETURN END SUBROUTINE SETA(A,MNWT,NF,N,N3) DIMENSION A(MNWT,NF,N) DO2I=1,NF DO 2 J=1,N3 2 A(J,I,2)=A(J,I,3) RETURN END 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 SUBROUTINE TORGB(A,N,IADDC) C INCODE IBMF C TORGERSON METHOD TO COMPUTE SCALAR PRODUCT MATRICES 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.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 SUBROUTINE MTINV (A,I1,I2,B,I3,I4,N,M) MTPK1115 C INCODE IBMF C FORTRAN DECK,LSTOU,STAB C MTINV REAL*8 B(I3,I4),BMAX,BMULT DIMENSION A(I1,I2) C MTPK1117 C THIS SUBROUTINE CALCULATES THE INVERSE OF MATRIX A BY THE GAUSS- MTPK1118 C JORDAN ELIMINATION SCHEME. ALL CALCULATIONS ARE DONE IN DOUBLE- MTPK1119 C PRECISION ARITHMETIC MTPK1120 C MTPK1121 EQUIVALENCE (BMAX,BMULT) MTPK1123 C MTPK1124 C SINGLE TO DOUBLE PRECISION TRANSFER MTPK1125 C MTPK1126 101 DO 103 I = 1,N MTPK1127 102 DO 103 J = 1,N MTPK1128 103 B(I,J) = A(I,J) MTPK1129 C MTPK1130 C FIND THE LARGEST ELEMENT IN THE N-K BY N-K LOWER RIGHT SUBMARTIX. MTPK1133 C MTPK1134 200 DO 430 K = 1,N MTPK1131 A(K,1) = FLOAT(K) MTPK1135 A(K,2) = A(K,1) 203 BMAX = DABS(B(K,K)) MTPK1137 210 DO 219 I = K,N MTPK1138 211 DO 219 J = K,N MTPK1139 IF (BMAX - DABS(B(I,J)) ) 213,219,219 213 BMAX = DABS ( B(I,J)) MTPK1141 214 A(K,1) = FLOAT(I) MTPK1142 215 A(K,2) = FLOAT(J) MTPK1143 219 CONTINUE MTPK1145 216 IF (BMAX - 1.D-38) 801,301,301 C MTPK1146 C EXCHANGE ROWS AND COLUMNS TO PUT B(I,J) ON DIAGONAL. MTPK1147 C MTPK1148 301 I = IFIX(A(K,1)) MTPK1149 302 J = IFIX(A(K,2)) MTPK1150 313 DO 314 M = 1,N MTPK1151 BMULT = B(I,M) MTPK1152 B(I,M) = B(K,M) MTPK1153 314 B(K,M) = BMULT MTPK1155 303 DO 304 M = 1,N MTPK1156 BMULT = B(M,J) MTPK1157 B(M,J) = B(M,K) MTPK1158 304 B(M,K) = BMULT MTPK1159 C MTPK1160 C ACTUAL INVERSION STEP. MTPK1161 C BEGIN ROW ITERATION MTPK1162 C IGNORE ROW K 401 DO 425 I = 1,N MTPK1164 C 402 IF (I - K) 405,425,405 C SET ROW MULTIPLIER. MTPK1167 405 BMULT = B(I,K) / B(K,K) C MODIFY ROW ELEMENTS. MTPK1169 410 DO 420 J = 1,N MTPK1170 C IGNORE COLUMN K MTPK1171 411 IF (J - K) 414,412,414 412 B(I,J) = -BMULT MTPK1173 413 GO TO 420 MTPK1174 414 B(I,J) = B(I,J) - B(K,J) * BMULT MTPK1175 420 CONTINUE MTPK1176 425 CONTINUE MTPK1177 C MTPK1178 C DIVIDE PIVOT ROW BY PIVOT ELEMENT MTPK1179 C MTPK1180 426 BMULT = 1.D0/B(K,K) 427 B(K,K) = 1.D0 C ACTUAL DIVISION MTPK1183 428 DO 429 J = 1,N MTPK1184 429 B(K,J) = B(K,J) * BMULT MTPK1185 430 CONTINUE MTPK1186 C MTPK1187 C REARRANGE AND REASSEMBLE INVERSE MATRIX. MTPK1188 C MTPK1189 501 DO 512 IJK = 1,N MTPK1190 502 I = N-IJK+1 MTPK1191 503 L = IFIX(A(I,2)) MTPK1192 504 J = IFIX(A(I,1)) MTPK1193 505 DO 508 M = 1,N MTPK1194 506 BMULT = B(I,M) MTPK1195 507 B(I,M) = B(L,M) MTPK1196 508 B(L,M) = BMULT MTPK1197 509 DO 512 M = 1,N MTPK1198 510 BMULT = B(M,I) MTPK1199 511 B(M,I) = B(M,J) MTPK1200 512 B(M,J) = BMULT MTPK1201 C MTPK1202 C MOVE A INVERSE BACK INTO A MTPK1203 C MTPK1204 601 DO 603 I = 1,N MTPK1205 602 DO 603 J = 1,N MTPK1206 603 A(I,J) = B(I,J) 701 M = 1 MTPK1208 702 RETURN MTPK1209 C MTPK1210 C ERROR RETURN AT SOME TIME BMAX WAS .LT. 0.1D-38 C MTPK1212 801 M = 2 MTPK1213 802 RETURN MTPK1214 END SUBROUTINE DIMA(A,MNWT,NFO,N) DIMENSION A(MNWT,NQ) COMMON NF,NQ NN=NFO L=NF DO 8 M=2,N DO 6 I=1,NF L=L+1 NN=NN+1 DO 6 J=1,MNWT 6 A(J,L)=A(J,NN) NN=NN+1 8 CONTINUE RETURN END SUBROUTINE SORT (A, N, B, C, D, K, SWITCH ) 75800000 CSORT 75700000 C 75900000 C THIS ROUTINE SORTS INPUT ARRAY 'A' AND REARRANGES, OPTIONALLY, 76000000 C ARRAYS 'B', 'C', AND 'D', IN ORDER CORRESPONDING TO 'A'. 76100000 C N = NUMBER OF ITEMS IN 'A' (AND 'B', 'C', 'D', IF USED) 76200000 C K = 0--SORT 'A' ONLY, 1--REARRANGE 'B', 2--REARRANGE 'B' AND 'C', 76300000 C 3--REARRANGE 'B', 'C', AND 'D'. 76400000 C IF 'SWITCH' IS POSITIVE, SORT WILL BE IN ASCENDING ORDER, 76500000 C IF ZERO OR NEGATIVE, IN DESCENDING ORDER. 76600000 C ALGORITHM FROM CACM, JULY 1959, PAGE 30 BY D. L. SHELL 76700000 C 76800000 DIMENSION A(1800), B(1800), C(1800), D(1800) 76900000 INTEGER SWITCH 77000000 105 KP1=K+1 77100000 IF(N.LE.1) GO TO 999 77200000 M = 1 77300000 106 M = M + M 77400000 IF( M .LT. N ) GO TO 106 77500000 M = M - 1 77600000 994 M = M/2 77700000 IF(M.EQ.0) GO TO 999 77800000 KK = N-M 77900000 J = 1 78000000 992 I = J 78100000 996 IM = I + M 78200000 IF(SWITCH) 810,810,800 78300000 800 IF(A(I).GT.A(IM)) GO TO 110 78400000 GO TO 995 78500000 810 IF(A(I).LT.A(IM)) GO TO 110 78600000 995 J = J+1 78700000 IF(J.GT.KK) GO TO 994 78800000 GO TO 992 78900000 110 TEMP=A(I) 79000000 A(I) = A(IM) 79100000 A(IM) = TEMP 79200000 GO TO ( 140, 130, 120, 115), KP1 79300000 115 TEMP = D(I) 79400000 D(I) = D(IM) 79500000 D(IM) = TEMP 79600000 120 TEMP=C(I) 79700000 C(I) = C(IM) 79800000 C(IM) = TEMP 79900000 130 TEMP=B(I) 80000000 B(I) = B(IM) 80100000 B(IM) = TEMP 80200000 140 I = I-M 80300000 IF(I.LT.1) GO TO 995 80400000 GO TO 996 80500000 999 RETURN 80600000 END 80700000 SUBROUTINE RANDU(IX,IY,YFL) IY=IX*65539 IF(IY)5,6,6 5 IY=IY+2147483647+1 6 YFL=YFL*0.4656613E-9 RETURN END