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