C nindscal.f C This software implements the method called INDSCAL. C The author of this software is 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----+----@----+----@----+----@----+----@----+----@----+----@----+----@ $ FORTREX C NINDSCAL- PERFORMS INDIVIDUAL DIFFERENCES ANALYSIS USING NONMETRIC C PROCEDURE. NINDSCAL DOES NOT PERFORM CANDECOMP ANALYSIS C THEREFORE EACH SUBJECTS DATA MATRIX MUST BE SYMMETRIC. C 7/10/72 VERSION C BY JIH JIE CHANG C PARAMETER DEFINITIONS (SEE >HOW TO USE INDSCAL...> PAPER FOR C DETAILED DESCRIPTION.) C N - NUMBER OF MATRICES (MAX.=3). IN GENERAL MATRIX 1 IS TREATED C AS SBJECT WEIGHT MATRIX AND MATRICES 2 AND 3 ARE STIMULI. C MAXDIM, MINDIM- NUMBER OF DIMENSIONS. (MAX.=5). UNLIKE INDSCAL C NINDSCAL CAN ONLY PROVIDE ONE SET OF SOLUTION IN A C DESIGNATED NUMBER OF DIMENSIONS FOR EACH SET OF DATA, C THEREFORE MAXDIM AND MINDIM MUST BE THE SAME. 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. =50) 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 DIMENSIONAL SOLUTION SUCCESSIVELY. C NOT AVAILABLE FOR NINDSCAL 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 EVERY MONOTONE ITERATION. 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 MONREG - 0, DO NOT PERFORM NONMETRIC ANALYSIS C MONREG - 1, DO SIMPLE MONOTONE REGRESSION C MONREG - 2, DO BLOCK MONOTONE REGRESSION, USE PRIMARY APPROACH C MONREG - 3, DO BLOCK MONOTONE REGRESSION, USE SECONDARY APPROACH C ITMONO -NUMBER OF NON-METRIC ITERATIONS C IBMONO - 0, BEGIN WITH METRIC ANALYSIS C IBMONO - 1, BEGIN WITH MONOTONE REGRESSION, MUST READ IN INITIAL C SUBJECTS WEIGHT MATRIX AND THE STIMULUS MATRIX. (SEE C SUBROUTINE READA) C MPLOT - 0, NO MONOTONE PLOTS C MPLOT - 1, DO MONOTONE PLOTS FOR EACH INDIVIDUAL AT THE END OF C VARIABLES C Y - DATA. Y IS A 3 SUBSCRIPTED VARIABLE. C Y(N1,N2,N3), WHERE N1, N2 ...REFERS TO ELEMENTS C IN MATRICES 1 TO 3. (MAX. (N1*N2*N3) =6250) C (MAX. N(I)=25) C A - COORDINATE MATRIX. C A(I,J,K), WHERE I - ELEMENTS WITHIN MATRICES (MAX. =25) C J - DIMENSIONS. (MAX. =5) C K - MATRIX. (MAX. =3) C MAXIMUM ((MAX.I)*J)=125 C NWT - SIZE OF MATRICES (MAX.=25) C TIT - TITLE, USE NO MORE THAN ONE CARD. DIMENSION Y(6250),A(375),AA(375),NWT(3),TIT(14),BNUM(25) DIMENSION KC(3000),DATA(300),DHAT(625),LBLOCK(3000),YY(3000) DATA TIT(1),TIT(14)/6H(79H ,6H )/ C SORTFL IS A SYSTEM SORT PACKAGE ROUTINE. IT SETS THE MODE TO SORT C FLOATING POINT ARRAYS IN ASCENDING ORDER. CALL SORTFL 90 READ106,N,MAXDIM,MINDIM,IRDATA,ITMAX,ISET,IOY,IDR,ISAM,IPUNSP,IRN IF(N.EQ.0) GOTO999 READ106,MONREG,ITMONO,IBMONO,MPLOT C IN NINDSCAL, N IS ALWAYS 3, IOY MUST BE 0 AND MAXDIM MUST EQUAL TO C MINDIM. HOWEVER, THEY ARE KEPT HERE TO BE IN CORRESPONDENCE WITH C THE INDSCAL PARAMETERS. MAXDIM=MINDIM NF=MAXDIM N=3 IOY=0 NA=2 NONM=MONREG-1 IMONO=0 IF(ISAM.EQ.1) ISET=0 READ101,(TIT(I),I=2,13) 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.125) 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. 6250) GOTO6 PRINT110 STOP 6 NDIM=0 MONOIT=0 N3=NWT(3) N2=NWT(2) N1=NWT(1) ND=(N2*(N2-1))/2 13 INORM=ISET MAXIT=ITMAX ISTI=ISAM PRINT104 PRINT102,(TIT(I),I=2,13) PRINT107 PRINT108 PRINT105,N,NF,IRDATA,MAXIT,INORM,IOY,IDR,ISTI,IPUNSP,IRN,CRIT PRINT113 PRINT114,MONREG,ITMONO,IBMONO PRINT109,(NWT(I),I=1,N) PRINT107 4 IF(ISTI.NE.0) MAXIT=0 8 IF (IRN.NE.0) CALL RAND1(IRN) 44 IF(IRDATA)41,42,41 41 CALL RDATA(Y,N1,N2,N3,IRDATA,IPUNSP,NONM,MONOIT,KC,ND,DATA,DHAT, 1LBLOCK,YY) IF(IBMONO)43,42,43 43 CALL READA(A,MNWT,NF,N1,N2) GOTO52 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 THE REMAINING C MATRICES 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) 52 NDIM=1 IF(NONM.LT.0) GOTO90 CALL NONMET(A,MNWT,NF,N1,N3,DATA,DHAT,NONM,KC,ND,MONOIT,ITMONO,Y, 1 LBLOCK,AA,YY,TIT,MPLOT) MONOIT=MONOIT&1 IF(MONOIT.GT.ITMONO) GOTO90 IBMONO=0 GOTO13 100 FORMAT(18I4) 101 FORMAT(12A6) 102 FORMAT(1H012A6) 103 FORMAT(10F7.3) 104 FORMAT(69H1NINDSCAL- NONMETRIC VERSION OF INDSCALS, ONLY UP TO 3 1WAY 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(49H0THE PRODUCT OF MATRIX SIZE EXCEEDS LIMIT OF 6250) 111 FORMAT(52H0(MAX. MATRIX SIZE * DIMENSION) EXCEEDS LIMIT OF 125) 113 FORMAT(21H0NONM ITMONO IBMONO) 114 FORMAT(I4,2I7) 999 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 NONMET(A,MNWT,NF,N1,N3,DATA,DHAT,LFITSW,KC,ND,MONOIT,IT 1MONO,Y, LBLOCK,NAS,YY,TIT,MPLOT) DIMENSION A(MNWT,NF,3),DATA(1),DHAT(1),KC(1),LBLOCK(ND,N1) DIMENSION Y(N1,N3,N3),NAS(300),YY(ND,N1),FMTS(9),TIT(1) DATA FMTS(1)/45H(21HDISTANCES SUBJECTI3,5X,9HITERATIONI3)/ NPASS=0 C COMPUTE DISTANCES FROM A3 PRINT100,MONOIT PRINT102 MM=0 MA=0 SSTR1=0. SSTR2=0. XND=ND XN1=N1 DO10I1=1,N1 L=0 DO2I=2,N3 II=I-1 DO2J=1,II L=L&1 DHAT(L)=0. DO2M=1,NF IF(A(I1,M,1))2,2,9 9 DHAT(L)=DHAT(L)&A(I1,M,1)*(A(I,M,3)-A(J,M,3))**2 2 CONTINUE DO3M=1,ND MM=MM&1 I=KC(MM) NAS(M)=I 3 DATA(M)=SQRT(DHAT(I)) IF(LFITSW.EQ.1) NPASS=1 CALL MFIT(DATA,DHAT,LBLOCK,ND,LFITSW,NPASS, NAS,PAS2) IF(MPLOT.EQ.0) GOTO16 CALL PXY(DATA,DHAT,YY(1,I1),ND,I1,TIT,FMTS,MONOIT) C MONOTONE FUNCTION IS STORED IN DHAT C COMPUTE STRESS 16 SUM=0. DO12I=1,ND 12 SUM=SUM&DATA(I) DMEAN=SUM/XND 11 SUMN=0. SUMD=0. SUMB=0. DO13I=1,ND SUMN=SUMN&(DATA(I)-DHAT(I))**2 15 SUMB=SUMB&DATA(I)**2 14 SUMD=SUMD&(DATA(I)-DMEAN)**2 13 CONTINUE STR1=SUMN/SUMB STR2=SUMN/SUMD PRINT103,I1,STR1,STR2 SSTR1=SSTR1&STR1 SSTR2=SSTR2&STR2 C REARRANGE DHAT TO ORIGINAL ORDER AND STORE IN DATA IF(MONOIT.EQ.ITMONO) GOTO10 6 DO4M=1,ND MA=MA&1 IF(LFITSW.EQ.1) GOTO67 I=KC(MA) GOTO4 67 I=NAS(M) 4 DATA(I)=DHAT(M) L=0 DO7I=2,N3 II=I-1 DO7J=1,II L=L&1 7 Y(I1,I,J)=DATA(L) 10 CONTINUE SSQ1=SSTR1/XN1 SSQ2=SSTR2/XN1 PRINT105 PRINT101,SSQ1 PRINT104,SSQ2 IF(MONOIT.EQ.0) GOTO20 IF(SSQ1.GE.PSQ1.AND.SSQ2.GE.PSQ2) GOTO21 20 PSQ1=SSQ1 PSQ2=SSQ2 GOTO99 21 PRINT106 ITMONO=MONOIT 106 FORMAT(37H0STRESS VALUES GO UP, ITERATION STOPS) 99 RETURN 100 FORMAT(19H0MONOTONE ITERATIONI3) 101 FORMAT(6H SSQ1=E15.6) 102 FORMAT(25H0INDIVIDUAL STRESS VALUES/5H SUBJ,6X,4HSSQ1,11X,4HSSQ2) 103 FORMAT(I3,2X,2E15.6) 104 FORMAT(6H SSQ2=E15.6) 105 FORMAT(12H0MEAN STRESS) 107 FORMAT(I2,10F7.4/(2X,10F7.4)) END $ FORTRAN SUBROUTINE RDATA(Y,N1,N2,N3,IRDATA,IPUNSP,NONM,MONOIT,KC,ND,DATA, 1D,LBLOCK,YY) DIMENSION Y(N1,N2,N3),D(N3,N3),FMT(12),KC(1),DATA(1) DIMENSION LBLOCK(ND,N1),YY(ND,N1) 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 IF(MONOIT.EQ.0) READ100,FMT MM=1 4 DO50I1=1,N1 IF(MONOIT.EQ.0) GOTO51 DO53I=2,N3 II=I-1 DO53J=1,II D(I,J)=Y(I1,I,J) 53 D(J,I)=D(I,J) GOTO43 51 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 IF(IRDATA.EQ.7) GOTO42 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 42 IF(NONM.LT.0) GOTO44 MA=MM L=0 DO32I=2,N3 II=I-1 DO32J=1,II L=L&1 KC(MM)=L MM=MM&1 32 DATA(L)=D(I,J) C SORT IS A SYSTEM SORT PACKAGE ROUTINE. IT SORTS THE ARRAY DATA OF C ND ELEMENTS IN ASCENDING ORDER. KC IS TO BE ARRANGED IN THE SAME C ORDER AS DATA. CALL SORT(DATA(1),DATA(2),ND,KC(MA)) DO52I=1,ND 52 YY(I,I1)=DATA(I) IF(NONM.EQ.0) GOTO44 CALL BMONO(DATA,LBLOCK(1,I1),ND) 44 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 READA(A,MNWT,NF,N1,N2) DIMENSION A(MNWT,NF,3),FMT(12) READ100,FMT 100 FORMAT(12A6) DO2I=1,NF 2 READFMT,(A(J,I,1),J=1,N1) DO3I=1,NF 3 READFMT,(A(J,I,2),J=1,N2) DO4I=1,NF DO4J=1,N2 4 A(J,I,3)=A(J,I,2) RETURN END $ FORTRAN SUBROUTINE PXY(X,Z,Y,N,K,TIT,FMTH) DIMENSION X(1),Y(1),Z(1),FMTH(1),TIT(1) CALL MM(X,N,XMAX,XMIN) CALL MM(Y,N,YMAX,YMIN) CALL FRAME(TIT) CALL HORIZ(FMTH,K) CALL LABEL(XMAX,XMIN,YMAX,YMIN) CALL DRAW(Z,Y,N,1) CALL DRAW(X,Y,N,0,1H*) 20 RETURN END $ FORTRAN SUBROUTINE BMONO(SCR,LBLOCK,N2) DIMENSION SCR(1),LBLOCK(1) L=1 LL=1 L1=1 DO17I=2,N2 II=I-1 IF(SCR(I).EQ.SCR(II))GOTO19 DO18J=L1,LL 18 LBLOCK(J)=L L1=L1&L L=1 GOTO16 19 L=L&1 16 LL=LL&1 17 CONTINUE DO20J=L1,LL 20 LBLOCK(J)=L RETURN END $ FORTRAN $ INCODE IBMF SUBROUTINE MFIT(DIST,DHAT,LBLOCK,MM,LFITSW,NPASS,PAS1,PAS2) 403 C FIT PERFORMS WEIGHTED-LEAST-SQUARES MONOTONE REGRESSION C THIS ROUTINE FINDS THE VALUES FOR DHAT WHICH ARE MONOTONIC C AND WHICH BEST APPROXIMATE THE VALUES OF DIST, C IN THE SENSE THAT THE SUM OF THE SQUARED DEVIATIONS, C EACH ONE WEIGHTED BY THE CORRESPONDING VALUE IN WW , C IS A MINIMUM. C IT DIRECTLY USES SORT. 407 408 DIMENSION PAS1(1),PAS2(1),DIST(1),DHAT(1),LBLOCK(1) 411 C FORM FIRST APPROXIMATION TO CORRECT PARTITON 412 C IF LFITSW=0, DO PLAIN REGRESSION C IF LFITSW GT 0, DO BLOCK REGRESSION C IF LFITSW=1, USE PRIMARY APPROACH. THUS WE SORT EACH 414 C BLOCK OF EQUAL VALUES OF DATA ACCORDING TO DIST VALUES. 415 C THEN WE CREATE PARTITON INTO BLOCKS OF SIZE 1 TO START. 416 417 C IF LFITSW=2, USE SECONDARY APPROACH. THUS WE START WITH 418 C PARTITION INTO BLOCKS OF EQUAL DATA VALUES. 419 420 C WITHIN EACH BLOCK OF WHATEVER SIZE, THE FIRST DHAT VALUE 421 C GIVES THE WEIGHTED TOTAL OF THE DIST VALUES IN THAT BLOCK, 423 C SIMILARLY, WITHIN EACH BLOCK, THE FIRST LBLOCK 424 C VALUE AND THE LAST LBLOCK VALUE BOTH GIVE THE SIZE OF THE 425 C BLOCK. C BLOCKS OF SIZE 1 FORM A PARTIAL EXCEPTION. EVERYTHIN IS THE SAME C EXCEPT THAT THE SUM OF THE W VALUES IS NOT FOUND IN THE C LAST DHAT VALUE IN THE BLOCK. 427 IF(LFITSW.GT.0)GOTO100 DO90M=1,MM DHAT(M)=DIST(M) 90 LBLOCK(M)=1 GOTO500 100 MA=1 430 110 K=LBLOCK(MA) MB=MA&K-1 432 GO TO ( 200, 300 ), LFITSW 433 434 C PRIMARY APPROACH 435 436 200 IF(K-1) 9999, 220, 210 437 C IF K=1, SAVE SORTING TIME 438 210 CALL SORTFL(&1) IF(NPASS.EQ.0)CALL SORT(DIST(MA),DIST(MA&1),K) IF(NPASS.EQ.1)CALL SORT(DIST(MA),DIST(MA&1),K,PAS1(MA)) IF(NPASS.EQ.2)CALL SORT(DIST(MA),DIST(MA&1),K,PAS1(MA),PAS2(MA)) 440 C 'SORT' WILL SORT THE K ELEMENTS OF 'DIST' IN ALGEBRAIC 441 C ORDER. 442 C BECAUSE THE FINAL ARGUMENT IS ZERO, SORTING WILL BE 443 C ASCENDING OR DESCENDING ACCORDING TO THE VALUE OF SDSWIT 444 C SUPPLIED DURING THE CALL FOR SORT IN MDSCAL. 445 C THE ELEMENTS OF'IJ' AND 'DATA' WILL BE REARRANGED IN 446 C EXACTLY THE SAME ORDER AS THOSE OF 'DIST'. 447 C ALSO THE ELEMENTS OF 'WW'. 448 220 DO 230 M=MA,MB 449 DHAT(M)=DIST(M) LBLOCK(M)=1 230 CONTINUE GO TO 400 452 453 C SECONDARY APPROACH 454 455 300 CONTINUE TEMP1=0.0 458 DO 310 M=MA,MB 459 TEMP1=TEMP1&DIST(M) 310 CONTINUE DHAT(MA)=TEMP1 461 464 C PROCEED TO NEXT BOCK. QUERY. IS IT THE LAST 465 466 400 MA=MA&K 467 IF(MA-MM-1) 110, 500, 9999 468 469 C START MAIN COMPUTATION ************************************* 470 471 500 MA=1 472 510 LUD=2 473 NSATIS=0 474 520 K=LBLOCK(MA) 475 MB=MA&K-1 476 WT=FLOAT(K) 550 DAV = DHAT(MA) / WT GO TO ( 600, 700 ), LUD 478 479 C IS BLOCK DOWN-SATISFIED. IF NOT, MERGE 480 481 600 IF(MA-1) 9999, 630, 610 482 610 MBD=MA-1 483 KD=LBLOCK(MBD) 484 MAD=MBD-KD&1 485 WTD=FLOAT(KD) 617 DAVD = DHAT(MAD) / WTD IF( DAVD-DAV ) 630, 620, 620 487 620 KNEW=K&KD 488 LBLOCK(MAD)=KNEW 489 LBLOCK(MB)=KNEW 490 DTONEW = DHAT(MAD) & DHAT(MA) DHAT(MAD) = DTONEW NSATIS=0 494 MA=MAD 495 GO TO 800 496 630 NSATIS=NSATIS&1 497 GO TO 800 498 499 C IS BLOCK UP-SATISFIED. IF NOT, MERGE 500 501 700 IF(MB-MM) 710, 730, 9999 502 710 MAU=MB&1 503 KU=LBLOCK(MAU) 504 MBU=MAU&KU-1 505 WTU=FLOAT(KU) 717 DAVU = DHAT(MAU) / WTU IF( DAV-DAVU ) 730, 720, 720 507 720 KNEW=K&KU 508 LBLOCK(MA)=KNEW 509 LBLOCK(MBU)=KNEW 510 DTONEW = DHAT(MA) & DHAT(MAU) DHAT(MA) = DTONEW NSATIS=0 514 GO TO 800 515 730 NSATIS=NSATIS&1 516 GO TO 800 517 518 C PROCEED TO NEXT BLOCK IF READY. 519 520 800 LUD = 3-LUD 521 C QUERY. IS BLOCK BOTH UP AND DOWN SATISFIED. IF NOT, RETUR 522 IF(NSATIS-1) 520, 520, 810 523 C QUERY. IS THIS LAST BLOCK. IF NOT, GO ON TO NEXT BLOCK. 524 810 IF(MB-MM) 820, 900, 9999 525 820 MA=MB&1 526 GO TO 510 527 528 C MAIN COMPUTATION COMPLETE. PLACE ANSWERS IN DHAT. 529 530 900 MA=1 531 910 K=LBLOCK(MA) 532 MB=MA&K-1 533 IF(K-1) 9999, 940, 920 534 920 TEMP1=DHAT(MA)/FLOAT(K) DO 930 M=MA,MB 536 930 DHAT(M)=TEMP1 537 GO TO 945 940 DHAT(MA) = DIST(MA) 945 MA = MB & 1 IF(MA-MM-1) 910, 950, 9999 539 950 RETURN 540 541 C TROUBLE EXIT 542 543 9999 PRINT99 99 FORMAT(58H0MFIT DIAGNOSTIC. IMPOSSIBLE BRANCH TAKEN ON IF STATEMEN 1T.) STOP 546 547 END 584 $ FORTRAN SUBROUTINE MM(X,N,XMAX,XMIN) C FINDS THE MAXIMUM AND MINIMUM OF ARRAY X DIMENSION X(100) XMAX=X(1) XMIN=XMAX DO20I=2,N IF(X(I)-XMAX)15,20,10 10 XMAX=X(I) GOTO20 15 IF(X(I)-XMIN)16,20,20 16 XMIN=X(I) 20 CONTINUE RETURN 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 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 SSQ=0. $ 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) MTPK1115 DIMENSION A(I1,I2), B(I3,I4) MTPK1116 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 DOUBLE PRECISION B,BMAX,BMULT MTPK1122 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 SUBMATRIX. 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 MTPK1163 401 DO 425 I = 1,N MTPK1164 C IGNORE ROW K. MTPK1165 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 999 END MTPK1215