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