Trailing-Edge
-
PDP-10 Archives
-
decuslib10-12
-
43,50550/forml1.for
There are no other files named forml1.for in the archive.
SUBROUTINE DAPPR(A,RNP,PJ,APP,G,VEC,S1,S2,HHH,NAT,NF,MCASE,
1MPJ,MAPP,MPRIN,PCON,WCON,COSLT,NCYCLT)
C FACTOR ROTATION, DIRECT METHOD WITH APP WEIGHTING
C APP IS ARTIFICIAL PERSONAL PROBABILITY OF BEING IN A HYPERPLANE
C INPUT OR: APP = 1/(1 + ACON*(ABS(B) TO POWER PCON)
C PROGRAMMED BY LEDYARD R TUCKER, NOVEMBER 1977
C REVISED JANUARY 1978
C USES FORMAL SUBROUTINES: TITLE,INVRS,MATML,XTX,GENR8,SIMLN,COPY
C OTHER SUBROUTINES USED: XMATML, XPRINT, XXTX, XVARMX
C A IS ORIGINAL FACTOR MATRIX: INPUT/OUTPUT
C RNP IS FORMAL MATRIX N' OF HYPERPLANE NORMALS (COLUMNS): OUTPUT
C PJ IS MATRIX P DURING SOLUTION, IS CONVERTED TO B FOR OUTPUT
C SEE NOTE CONCERNING MPJ
C APP IS MATRIX OF APP: OUTPUT
C SEE NOTE CONCERNING MAPP
C G IS FORMAL MATRIX (NN') INVERSE DURING SOLUTION,
C IS CONVERTED TO PHI FOR OUTPUT
C VEC IS FORMAL MATRIX HAVING ROW VECTORS AS FOLLOWS:
C 1 MEAN ABSOLUTE PJ FOR EACH ROTATED FACTOR
C 2 SUM APP*(PJ**2) OVER ATTRIBUTES FOR EACH ROTATED FACTOR
C 3 SHIFT VECTOR IN ROTATION OF EACH FACTOR IN TURN
C 4 SQUARE ROOTS OF DIAGONAL G
C 5-I1 SCRATCH
C S1 AND S2 ARE SCRATCH FORMAL MATRICES
C ON OUTPUT S1 CONTAINS CORRELATION OF ITEMS AND FACTORS,
C S2 CONTAINS T'
C HHH IS A SCRATCH VECTOR OF LENGTH I1
C NAT IS NUMBER OF ATTRIBUTES: INPUT/OUTPUT
C NF IS NUMBER OF FACTORS: INPUT/OUTPUT
C MCASE IS CASE OF APP FUNCTION: INPUT/OUTPUT
C MUST BE SPECIFIED WHEN APP IS TO BE ITERATED
C = 1 FOR ONE SIDED FUNCTION
C = 2 FOR TWO SIDED FUNCTION
C MPJ IS INTEGER FLAG FOR INITIAL WEIGHTS HYPOTHESIZED: INPUT/OUTPUT
C = 0 FOR INITIAL WEIGHTS OBTAINED BY A NORMAL VARIMAX ROTATION
C = 1 FOR INITIAL HYPOTHESIS FACTOR WEIGHTS INPUT IN PJ
C DEFAULT = 0
C MAPP IS INTEGER FLAG CONCERNING ITERATION OF APP
C = 0 FOR APP ITERATION
C = 1 FOR FIXED APP AS INPUT IN APP
C DEFAULT = 0
C MPRIN IS INTEGER FLAG CONCERNING PRINTED OUTPUT
C = 0 FOR NO PRINTED OUTPUT
C = 1 FOR PRINTING MATRICES P AND APP FOR EACH CYCLE
C DEFAULT = 0
C PCON IS POWER CONSTANT IN APP FUNCTION:INPUT/OUTPUT
C DEFAULT = 10.0
C WCON IS PROPORTION OF MEAN ABSOLUTE FACTOR WEIGHT FOR APP = 0.5:
C INPUT/OUTPUT; DEFAULT = 0.5
C COSLT IS LIMIT FOR COS BETWEEN SUCCESSIVE TRIAL NORMALS IN TEST FOR
C CONVERGENCE: INPUT/OUTPUT; DEFAULT = 0.99999
C NCYCLT IS LIMIT ON NUMBER OF ITERATIVE CYCLES: INPUT/OUTPUT
C DEFAULT = 20
C NRMA IS THE NUMBER OF ROWS OF A IN CALLING PROGRAM SET TO I1
C NRMRNP IS THE NUMBER OF ROWS OF RNP IN CALLING PROGRAM SET TO I1
C NRMPJ IS THE NUMBER OF ROWS OF PJ IN CALLING PROGRAM SET TO I1
C NRCPJ IS THE NUMBER OF COL. OF PJ IN CALLING PROGRAM SET TO I1
C NRMAPP IS THE NUMBER OF ROWS OF APP IN CALLING PROGRAM SET TO I1
C NRMG IS THE NUMBER OF ROWS OF G IN CALLING PROGRAM SET TO I1
C NRMVEC IS THE NUMBER OF ROWS OF VEC IN CALLING PROGRAM SET TO I1
C NRMS1 IS THE NUMBER OF ROWS OF S1 IN CALLING PROGRAM SET TO I1
C NRMS2 IS THE NUMBER OF ROWS OF S2 IN CALLING PROGRAM SET TO I1
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION A(I1,I1),RNP(I1,I1),PJ(I1,I1),APP(I1,I1),G(I1,I1)
DIMENSION VEC(I1,I1),S1(I1,I1),S2(I1,I1),HHH(I1)
COMMON/B/I1
COMMON/ZINCR/IZINCR
CALL QINCR ('DAPPR')
NRMA = I1
NRMRNP = I1
NRMPJ = I1
NCMPJ = I1
NRMAPP = I1
NRMG = I1
NCMG = I1
NRMVEC = I1
NRMS1 = I1
NCMS1 = I1
NRMS2 = I1
NCMS2 = I1
CALL XDAPPR(NRMA,NRMRNP,NRMPJ,NCMPJ,NRMAPP,NRMG,NCMG,NRMVEC,
2 NRMS1,NCMS1,NRMS2,NCMS2,A,RNP,PJ,APP,G,VEC,S1,S2,HHH,NAT,NF,
3 MCASE,MPJ,MAPP,MPRIN,PCON,WCON,COSLT,NCYCLT)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XDAPPR(NRMA,NRMRNP,NRMPJ,NCMPJ,NRMAPP,NRMG,NCMG,
2 NRMVEC,NRMS1,NCMS1,NRMS2,NCMS2,A,RNP,PJ,APP,G,VEC,S1,S2,HHH,
3 NAT,NF,MCASE,MPJ,MAPP,MPRIN,PCON,WCON,COSLT,NCYCLT)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION A(NRMA,1),RNP(NRMRNP,1),PJ(NRMPJ,1),APP(NRMAPP,1)
DIMENSION G(NRMG,1),VEC(NRMVEC,1),S1(NRMS1,1),S2(NRMS2,1),HHH(1)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
CALL QINCR ('XDAPPR')
CALL XDIMCH(NRMA,NAT,'DAPPR',0)
CALL XDIMCH(NRMRNP,NF,'DAPPR',0)
CALL XDIMCH(NCMPJ,NAT,'DAPPR',0)
CALL XDIMCH(NRMAPP,NAT,'DAPPR',0)
CALL XDIMCH(NRMG,NF,'DAPPR',0)
CALL XDIMCH(NRMVEC,NAT,'DAPPR',0)
CALL XDIMCH(NRMS1,NAT,'DAPPR',0)
CALL XDIMCH(NRMS2,NF,'DAPPR',0)
CALL TITLE('DIRECT FACTOR ROTATION WITH APP WEIGHTING',41)
IF(NF.LE.1) GO TO 310
IF(MAPP.EQ.1) GO TO 3
IF(MCASE.EQ.1) GO TO 1
IF(MCASE.EQ.2) GO TO 2
CALL TITLE('FORM OF APP FUNCTION MUST BE SPECIFIED',38)
CALL TITLE('CALCULATIONS TERMINATED',23)
STOP
1 CALL TITLE('APP FUNCTION IS ONE SIDED',25)
GO TO 3
2 CALL TITLE('APP FUNCTION IS TWO SIDED',25)
3 CONTINUE
IF(MPJ.NE.0) GO TO 4
CALL TITLE('INITIAL HYPOTHESIZED FACTOR WEIGHTS OBTAINED BY VARIMA
1X ROTATION',64)
GO TO 5
4 CALL TITLE('INITIAL HYPOTHESIZED FACTOR WEIGHTS ARE INPUT',45)
5 CONTINUE
IF(MAPP.NE.0) GO TO 6
CALL TITLE('APP MATRIX IS TO BE ITERATED',28)
GO TO 7
6 CALL TITLE('MATRIX APP IS FIXED AS INPUT',28)
7 CONTINUE
PW = PCON
IF(PW.LT.1.0D0) PW = 1.0D1
WRITE(LUWN,910) PW
910 FORMAT('0','POWER CONSTANT EQUALS',F10.5)
WW = WCON
IF(WW.LE.0.0D0) WW = 0.5D0
WRITE(LUWN,911) WW
911 FORMAT('0','PROPORTION OF MEAN ABSOLUTE FACTOR WEIGHT FOR APP = 1/
12 EQUALS',F10.5)
WW = WW*2.0D0
COSLTW = COSLT
IF(COSLT.LE.0.0D0) COSLTW = 0.99999D0
T1 = 1.0D0 - COSLTW
WRITE(LUWN,912)T1
912 FORMAT('0','LIMIT ON COSINE BETWEEN SUCCESSIVE NORMALS FOR CONVERG
1ENCE EQUALS 1.0 -',E8.1)
NCYLTW = NCYCLT
IF(NCYCLT.LE.0) NCYLTW = 20
WRITE(LUWN,913) NCYLTW
913 FORMAT('0','LIMIT ON NUMBER OF CYCLES OF ITERATION EQUALS',I5)
NFM1 = NF - 1
C OBTAIN VARIMAX ROTATION WHEN MPJ = 0
MRAW = 1
IF(MPJ.EQ.0) CALL XVARMX(NRMA,NRMPJ,NCMPJ,A,PJ,NAT,NF,
1HHH,MRAW)
C OBTAIN INITIAL RNP AND G
CALL XXTX(NRMA,NRMS2,A,S2,NAT,NF)
CALL XINVRS(NRMS1,NCMS1,NRMS2,NCMS2,S2,S1,NF)
DO 10 J = 1,NF
DO 10 K = 1,NF
S2(J,K) = 0.0D0
DO 10 I = 1,NAT
10 S2(J,K) = S2(J,K) + A(I,J)*PJ(I,K)
CALL XMATML(NRMS1,NRMS2,NRMRNP,S1,S2,RNP,NF,NF,NF)
DO 20 J = 1,NF
T1 = 0.0D0
DO 15 I = 1,NF
15 T1 = T1 + RNP(I,J)**2
T1 = SQRT(T1)
DO 20 I = 1,NF
20 RNP(I,J) = RNP(I,J)/T1
CALL XXTX(NRMRNP,NRMS1,RNP,S1,NF,NF)
CALL XINVRS(NRMS1,NCMS1,NRMG,NCMG,S1,G,NF)
DO 22 J = 1,NF
22 VEC(4,J) = SQRT(G(J,J))
C OBTAIN INITIAL ACTUAL PJ
IF(MPJ.NE.0)CALL XMATML(NRMA,NRMRNP,NRMPJ,A,RNP,PJ,NAT,NF,NF)
IF(MAPP.NE.0) GO TO 42
C OBTAIN FOR EACH ROTATED FACTOR MEAN ABSOLUTE PJ
DO 30 J = 1,NF
VEC(1,J) = 0.0D0
DO 25 I = 1,NAT
25 VEC(1,J) = VEC(1,J) + ABS(PJ(I,J))
30 VEC(1,J) = VEC(1,J)/NAT
NTT = 0
23 CONTINUE
C OBTAIN AWL (LOG OF ACON)
T1 = 0.0D0
DO 35 J = 1,NF
35 T1 = T1 + VEC(1,J)*VEC(4,J)
T1 = WW*T1/NF
AWL = -PW*ALOG(T1)
C OBTAIN MATRIX APP
IF(MCASE.EQ.2) GO TO 39
DO 38 I = 1,NAT
DO 38 J = 1,NF
IF(PJ(I,J).LE.0) GO TO 37
T1 = PJ(I,J)*VEC(4,J)
T1 = AWL + PW*ALOG( ABS(T1))
APP(I,J) = 1.0D0/(1.0D0 + EXP(T1))
GO TO 38
37 APP(I,J) = 1.0D0
38 CONTINUE
GO TO 42
39 CONTINUE
DO 40 I = 1,NAT
DO 40 J = 1,NF
T1 = PJ(I,J)*VEC(4,J)
T1 = AWL + PW*ALOG( ABS(T1))
40 APP(I,J) = 1.0D0/(1.0D0 + EXP(T1))
C OBTAIN SUM OF WEIGHTED SQUARES VECTOR
42 DO 45 J = 1,NF
VEC(2,J) = 0.0D0
DO 45 I = 1,NAT
T1 = PJ(I,J)**2
45 VEC(2,J) = VEC(2,J) + APP(I,J)*T1
C INITIALIZE COUNTERS
MF = 1
NCYC = 0
NCOSH = 0
C START OF CYCLE
C TRIAL FOR EACH ROTATED FACTOR IN TURN
100 MFM1 = MF - 1
MFP1 = MF + 1
C PRINT WHEN MPRIN.EQ.1
IF(MF.GT.1) GO TO 1000
IF(MPRIN.NE.1) GO TO 1000
ITEM = NCYC + 1
WRITE(LUWN,1900) ITEM
1900 FORMAT('0','CYCLE',I6)
CALL TITLE('INITIAL MATRIX P',16)
CALL XPRINT(NRMPJ,PJ,NAT,NF,'F')
CALL TITLE('INITIAL MATRIX APP',18)
CALL XPRINT(NRMAPP,APP,NAT,NF,'F')
1000 CONTINUE
C WEIGHTED PRODUCT MATRIX
CALL XGENR8(NRMS2,0.0D0,S2,NF,NF)
DO 105 I = 1,NAT
DO 105 J = 1,NF
VEC(5,J) = PJ(I,J)*APP(I,MF)
DO 105 K = 1,J
105 S2(J,K) = S2(J,K) + VEC(5,J)*PJ(I,K)
DO 110 J = 1,NF
DO 110 K = 1,J
110 S2(K,J) = S2(J,K)
C MATRIX R, INCLUDE CLOSING OF RANKS TO ELIMINATE ROW AND COLUMN MF
IF(MF.EQ.1) GO TO 125
DO 120 J = 1,MFM1
DO 113 K = J,MFM1
S1(J,K) = S2(J,K)*G(MF,MF)
113 S1(K,J) = S1(J,K)
IF(MF.EQ.NF) GO TO 118
DO 116 K = MFP1,NF
KM = K - 1
S1(J,KM) = S2(J,K)*G(MF,MF)
116 S1(KM,J) = S1(J,KM)
118 CONTINUE
120 S1(J,J) = S1(J,J) + G(MF,MF)*VEC(2,J)
125 IF(MF.EQ.NF) GO TO 140
DO 135 J = MF,NFM1
JP = J + 1
DO 130 K = MF,NFM1
KP = K + 1
130 S1(J,K) = S2(JP,KP)*G(MF,MF)
135 S1(J,J) = S1(J,J) + G(MF,MF)*VEC(2,JP)
140 CONTINUE
C COLUMN VECTOR C, CLOSING RANKS TO EXCLUDE ROW MF
IF(MF.EQ.1) GO TO 150
DO 145 J = 1,MFM1
145 S1(J,NF) = G(J,MF)*VEC(2,J) - G(MF,MF)*S2(J,MF)
150 IF(MF.EQ.NF) GO TO 160
DO 155 J = MF,NFM1
JP = J + 1
155 S1(J,NF) = G(JP,MF)*VEC(2,JP) - G(MF,MF)*S2(JP,MF)
160 CONTINUE
IF(NF.GT.2) GO TO 165
S2(1,1) = S1(1,2)/S1(1,1)
GO TO 170
165 CALL XSIMLN(NRMS1,NRMS2,S1,S2,NFM1,NF)
C PLACE S ENTRIES IN VEC(3,J)
170 IF(MF.EQ.1) GO TO 180
DO 175 J = 1,MFM1
175 VEC(3,J) = S2(J,1)
180 VEC(3,MF) = 1.0D0
IF(MF.EQ.NF) GO TO 190
DO 185 J = MF,NFM1
JP = J + 1
185 VEC(3,JP) = S2(J,1)
190 CONTINUE
C OBTAIN NEW NORMAL
DO 210 J = 1,NF
VEC(6,J) = RNP(J,MF)
IF(MF.EQ.1) GO TO 200
DO 195 K = 1,MFM1
195 VEC(6,J) = VEC(6,J) + RNP(J,K)*VEC(3,K)
200 IF(MF.EQ.NF) GO TO 210
DO 205 K = MFP1,NF
205 VEC(6,J) = VEC(6,J) + RNP(J,K)*VEC(3,K)
210 CONTINUE
RL = 0.0D0
DO 215 J = 1,NF
215 RL = RL + VEC(6,J)**2
RL = SQRT(RL)
DO 220 J = 1,NF
220 VEC(6,J) = VEC(6,J)/RL
C OBTAIN COS
COS = 0.0D0
DO 225 J = 1,NF
225 COS = COS + RNP(J,MF)*VEC(6,J)
C CHECK COS AND DO NOT MAKE ROTATION FOR LARGE COS
IF(COS.LT.COSLTW) GO TO 230
NCOSH = NCOSH + 1
IF(NCOSH.GE.NF) GO TO 305
GO TO 295
230 NCOSH = 0
C MOVE NORMAL INTO PLACE
DO 235 J = 1,NF
235 RNP(J,MF) = VEC(6,J)
C UPDATE G
DO 245 J = 1,NF
IF(J.EQ.MF) GO TO 245
DO 240 K = 1,NF
240 G(J,K) = G(J,K) - VEC(3,J)*G(MF,K)
245 CONTINUE
DO 255 K = 1,NF
IF(K.EQ.MF) GO TO 255
DO 250 J = 1,K
250 G(J,K) = G(J,K) - G(J,MF)*VEC(3,K)
255 CONTINUE
DO 260 K = 1,NF
DO 260 J = 1,K
260 G(K,J) = G(J,K)
G(MF,MF) = G(MF,MF)*RL
DO 265 J = 1,NF
G(MF,J) = G(MF,J)*RL
265 G(J,MF) = G(MF,J)
DO 267 J = 1,NF
267 VEC(4,J) = SQRT(G(J,J))
C OBTAIN NEW COLUMN MF OF PJ
DO 270 I = 1,NAT
PJ(I,MF) = 0.0D0
DO 270 J = 1,NF
270 PJ(I,MF) = PJ(I,MF) + A(I,J)*RNP(J,MF)
IF(MAPP.NE.0) GO TO 287
VEC(1,MF) = 0.0D0
DO 275 I = 1,NAT
275 VEC(1,MF) = VEC(1,MF) + ABS(PJ(I,MF))
VEC(1,MF) = VEC(1,MF)/NAT
C NEW AWL AND COLUMN MF OF APP
T1 = 0.0D0
DO 280 J = 1,NF
280 T1 = T1 + VEC(1,J)*VEC(4,J)
AWL = -PW*ALOG(WW*T1/NF)
IF(MCASE.EQ.2) GO TO 284
DO 283 I = 1,NAT
IF(PJ(I,MF).LE.0) GO TO 282
T1 = PJ(I,MF)*VEC(4,MF)
T1 = AWL + PW*ALOG( ABS(T1))
APP(I,MF) = 1.0D0/(1.0D0 + EXP(T1))
GO TO 283
282 APP(I,MF) = 1.0D0
283 CONTINUE
GO TO 287
284 CONTINUE
DO 285 I = 1,NAT
T1 = PJ(I,MF)*VEC(4,MF)
T1 = AWL + PW*ALOG( ABS(T1))
285 APP(I,MF) = 1.0D0/(1.0D0 + EXP(T1))
C NEW SUM OF WEIGHTED PJ SQUARED
287 VEC(2,MF) = 0.0D0
DO 290 I = 1,NAT
T1 = PJ(I,MF)**2
290 VEC(2,MF) = VEC(2,MF) + APP(I,MF)*T1
C PREPARE FOR NEXT FACTOR
295 MF = MF + 1
IF(MF.LE.NF) GO TO 100
NCYC = NCYC + 1
IF(NCYC.GE.NCYLTW) GO TO 300
MF = 1
GO TO 100
300 IF(NTT.EQ.1) WRITE(LUWN,914) NCYC
914 FORMAT('0','CONVERGENCE WAS NOT OBTAINED IN',I4,2X,'CYCLES')
305 IF(NTT.EQ.1) GO TO 315
NTT = NTT + 1
WW = WW/2.0D0
GO TO 23
310 IF(NF.NE.1) GO TO 340
RNP(1,1) = 1.0D0
G(1,1) = 1.0D0
VEC(4,1) = 1.0D0
DO 312 I = 1,NAT
312 PJ(I,1) = A(I,1)
C PRINT RESULTS
315 CALL TITLE('RESULTS OBTAINED',16)
CALL TITLE (' ',1)
CALL XXTX(NRMRNP,NRMS1,RNP,S1,NF,NF)
CALL TITLE('MATRIX P PROJECTIONS ON NORMALS (STRUCTURE LOADINGS O
1N REFERENCE VECTORS)',74)
CALL TITLE ('PARTIAL CORRELATIONS OF VARIABLES AND FACTOR - HOLDING
1G OTHER FACTORS CONSTANT', 77)
CALL XPRINT(NRMPJ,PJ,NAT,NF,'F')
CALL XGENR8(NRMS1,0.0D0,S1,NF,NF)
DO 320 I = 1,NF
320 S1(I,I) = 1.0D0/VEC(4,I)
DO 325 I = 1,NAT
DO 325 J = 1,NF
325 PJ(I,J) = PJ(I,J)*VEC(4,J)
CALL TITLE('FACTOR WEIGHT MATRIX B (PATTERN LOADINGS ON PRIMARY V
1ECTORS)',61)
CALL XPRINT(NRMPJ,PJ,NAT,NF,'F')
CALL XMATML(NRMRNP,NRMG,NRMS2,RNP,G,S2,NF,NF,NF)
DO 335 I = 1,NF
DO 335 J = 1,NF
335 S2(I,J) = S2(I,J)*S1(J,J)
DO 330 I = 1,NF
DO 330 J = 1,NF
330 G(I,J) = S1(I,I)*G(I,J)*S1(J,J)
CALL TITLE('FACTOR INTERCORRELATION MATRIX PHI',34)
CALL XPRINT(NRMG,G,NF,NF,'F')
CALL TITLE('MATRIX T'' PRIMARY VECTORS COLUMNS',34)
CALL XPRINT(NRMS2,S2,NF,NF,'F')
CALL TITLE ('ZERO ORDER CORRELATIONS OF VARIABLES AND FACTOR', 47)
CALL XMATML(NRMPJ,NRMG,NRMS1,PJ,G,S1,NAT,NF,NF)
CALL XPRINT (NRMS1, S1, NAT, NF, 'F')
340 IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE DISCR(SH,SD,GPM,F,G,CV,TM1,TM2,NAT,NG,NDFH,NDFD,MOD,
1NPS,NPE,ISWP)
C COMPUTES DISCRIMINANT COMPOSITE SOLUTION
C SOLUTION SPACE RESTRICTED TO NP PRINCIPAL AXES OF ESTIMATED
C TOTAL COVARIANCE MATRIX, THETA, DERIVED FROM COMPONENTS OF
C COVARIANCE.
C THETA IS A REPARAMETERIZATION OF THE TOTAL SSCP MATRIX.
C E.G., FOR MOD = 1, THETA = D(SH + SD)D,
C WHERE D IS A DIAGONAL MATRIX CONTAINING THE DIAGONAL ENTRIES
C OF (SH + SD) RAISED TO THE -1/2 POWER. THETA WILL APPEAR TO
C BE IN A FORM SIMILAR TO A CORRELATION MATRIX.
C PROGRAM WRITTEN BY LEDYARD R TUCKER, MAY 1978
C USES FORMAL; MATRICES MUST BE AT LEAST (10,10), IF NAT IS
C LARGER THAN 9, THE MATRICES MUST BE AT LEAST (NAT+1,NAT+1)
C PROBLEM COEFFICIENTS; INPUT/OUTPUT
C NAT IS NUMBER OF ATTRIBUTES (VARIABLES)
C NG IS NUMBER OF GROUPS (ELEMENTS IN HYPOTHESIS EFFECT)
C NDFH IS NUMBER OF DEGREES OF FREEDOM FOR HYPOTHESIS
C NDFD IS NUMBER OF DEGREES OF FREEDOM FOR DENOMINATOR
C MOD IS MODEL IDENTIFICATION NUMBER
C = 1 WHEN TOTAL SAMPLE IS DIVIDED INTO SAMPLES OF
C SUBPOPULATIONS IN PROPORTION TO RELATIVE FREQUENCIES
C OF THE SUBPOPULATIONS IN THE TOTAL POPULATION
C = 2 WHEN THE TOTAL SAMPLE IS DRAWN AT RANDOM FROM THE TOTAL
C POPULATION
C NPS IS STARTING NUMBER OF PRINCIPAL AXES OF ESTIMATED THETA
C NPE IS ENDING NUMBER OF PRINCIPAL AXES OF ESTIMATED THETA
C SEPARATE SOLUTIONS ARE OBTAINED USING NP PRINCIPAL AXES
C OF ESTIMATED THETA FROM NPS TO NPE IN STEPS OF 1
C WORKING VALUES, NPSW AND NPEW, OF NPS AND NPE ARE ADJUSTED
C WHEN NECESSARY AS FOLLOWS:
C LET NPAM BE THE MAXIMUM USABLE NUMBER OF PRINCIPAL AXES
C IF NPS EQUALS 0 (OR LESS), NPSW IS SET EQUAL TO NPAM
C IF NPS IS GREATER THAN NPAM, NPSW IS SET EQUAL TO NPAM
C IF NPE IS LESS THAN NPSW, NPEW IS SET EQUAL TO NPSW
C IF NPE IS GREATER THAN NPAM, NPEW IS SET EQUAL TO NPAM
C MATRICES
C SH IS HYPOTHESIS EFFECT SSCP MATRIX; INPUT/OUTPUT
C SD IS DENOMINATOR EFFECT SSCP MATRIX; INPUT/OUTPUT
C GPM IS GROUP MEANS MATRIX (GROUPS BY HYPOTHESIS EFFECT);
C INPUT/OUTPUT
C ROW FOR EACH GROUP PLUS ONE ROW FOR TOTAL MEAN
C F, G, CV, TM1, TM2 ARE STORAGE MATRICES
C PRINTING CONTROL
C ISWP = 0 FOR STANDARD RESULTS
C = 1 FOR ADDITIONAL INFORMATION:
C MATRICES THETA-TILDA , L , F , G , A
C COLUMNS OF CV ARE USED AS FOLLOWS
C 1 = 1/D
C 2 = 1/SQRT(DELTA)
C 3 = PHI
C 4 = 1/SQRT(1 - PHI)
C 5 = LAMBDA
C 6 = LN(1 + LAMBDA)
C 7 = V STATISTIC
C 8 = SQRT(SD(J,J)/NDFD)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION SH(I1,I1),SD(I1,I1),GPM(I1,I1),F(I1,I1),G(I1,I1)
DIMENSION CV(I1,I1),TM1(I1,I1),TM2(I1,I1)
COMMON/B/I1
COMMON/ZINCR/IZINCR
CALL QINCR ('DISCR')
NRMSH = I1
NRMSD = I1
NRMGPM = I1
NRMF = I1
NRMG = I1
NRMCV = I1
NRMTM1 = I1
NRMTM2 = I1
CALL XDISCR(NRMSH,NRMSD,NRMGPM,NRMF,NRMG,NRMCV,NRMTM1,NRMTM2,
2SH,SD,GPM,F,G,CV,TM1,TM2,NAT,NG,NDFH,NDFD,MOD,NPS,NPE,ISWP)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XDISCR(NRMSH,NRMSD,NRMGPM,NRMF,NRMG,NRMCV,NRMTM1,
2NRMTM2,SH,SD,GPM,F,G,CV,TM1,TM2,NAT,NG,NDFH,NDFD,MOD,NPS,NPE,
3ISWP)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION SH(NRMSH,1),SD(NRMSD,1),GPM(NRMGPM,1),F(NRMF,1)
DIMENSION G(NRMG,1),CV(NRMCV,1),TM1(NRMTM1,1),TM2(NRMTM2,1)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
CALL QINCR ('XDISCR')
CALL XDIMCH(NRMSH,NAT,'DISCR',0)
CALL XDIMCH(NRMSD,NAT,'DISCR',0)
CALL XDIMCH(NRMTM1,NAT,'DISCR',0)
CALL XDIMCH(NRMTM2,NAT,'DISCR',0)
CALL XDIMCH(NRMF,NAT,'DISCR',0)
CALL XDIMCH(NRMGPM,NPS,'DISCR',0)
CALL XDIMCH(NRMG,NAT,'DISCR',0)
CALL XDIMCH(NRMCV,NAT,'DISCR',0)
WRITE(LUWN,900) NAT
900 FORMAT('0','DISCRIMINANT SOLUTION FOR PROBLEM HAVING',I4,2X,'ATTRI
1BUTES')
WRITE(LUWN,901) NDFH,NDFD
901 FORMAT('0',5X,'DEGREES OF FREEDOM: HYPOTHESIS',I4,4X,'DENOMINATOR'
1,I4)
NGP = NG + 1
NPSW = NPS
NPEW = NPE
EPS = 1.0D-3
C OBTAIN NT: FOR MOD = 1, NT = NDFH + NDFD + 1; FOR MOD = 2, NT = NDFH + NDFD
NT = NDFH + NDFD
IF(MOD.EQ.1) NT = NT + 1
RNT = NT
RTNT = SQRT(RNT)
C OBTAIN NT*THETA-HAT, THETA-TILDA, EIGENSOLUTION
IF(MOD.EQ.1) GO TO 5
CALL XADD(NRMSH,NRMSD,NRMTM1,SH,SD,TM1,NAT,NAT)
GO TO 7
5 T1 = NDFD + 1.0D0
T1 = T1/NDFD
CALL XSCAML(NRMSD,NRMTM2,T1,SD,TM2,NAT,NAT)
CALL XADD(NRMSH,NRMTM2,NRMTM1,SH,TM2,TM1,NAT,NAT)
CALL TITLE ('TOTAL SSCP MATRIX', 17)
CALL XPRINT (NRMTM1, TM1, NAT, NAT, 'F')
7 DO 10 J = 1,NAT
CV(J,1) = 1.0D0/SQRT(TM1(J,J))
CV(J,8) = SQRT(SD(J,J)/NDFD)
DO 10 K = 1,J
TM1(J,K) = CV(J,1)*TM1(J,K)*CV(K,1)
10 TM1(K,J) = TM1(J,K)
CALL XEIGEN(NRMTM1,NRMF,NRMTM2,TM1,F,TM2,NAT,NAT)
C PRINT RESULTS
IF(ISWP.EQ.0) GO TO 20
CALL TITLE('MATRIX THETA-TILDA',19)
CALL XPRINT(NRMTM1,TM1,NAT,NAT,'F')
CALL TITLE('MATRIX L OF EIGENVECTORS OF THETA-TILDA',41)
CALL XPRINT(NRMF,F,NAT,NAT,'F')
20 CONTINUE
CALL TITLE('EIGENVALUES OF THETA-TILDA = DELTAS',36)
CALL PRDIF(TM2,NAT,I2)
C DETERMINE MAXIMUM NUMBER OF PRINCIPAL AXES TO USE
C ALSO 1/SQRT(DELTA)
NPAM = 0
DO 25 J = 1,NAT
IF(TM2(J,J).LT.EPS) GO TO 30
NPAM = NPAM + 1
25 CV(J,2) = 1.0D0/SQRT(TM2(J,J))
30 WRITE(LUWN,902) NPAM
902 FORMAT('0','MAXIMUM POSSIBLE NUMBER OF PRINCIPAL AXES OF THETA-TIL
1DA',I5)
C RESET NPSW AND NPEW WHEN NECESSARY
IF(NPSW.LT.1.OR.NPSW.GT.NPAM) NPSW = NPAM
IF(NPEW.LT.NPSW) NPEW = NPSW
IF(NPEW.GT.NPAM) NPEW = NPAM
C DETERMINE MATRICES F AND G
DO 35 J = 1,NAT
DO 35 K = 1,NPAM
35 F(J,K) = CV(J,1)*F(J,K)*CV(K,2)
T1 = NDFH
T1 = T1/NDFD
CALL XSCAML(NRMSD,NRMTM1,T1,SD,TM1,NAT,NAT)
CALL XSUB(NRMSH,NRMTM1,NRMTM2,SH,TM1,TM2,NAT,NAT)
CALL XMATML(NRMTM2,NRMF,NRMTM1,TM2,F,TM1,NAT,NAT,NPAM)
DO 40 J = 1,NPAM
DO 40 K = 1,J
G(J,K) = 0.D0
DO 37 L = 1,NAT
37 G(J,K) = G(J,K) + F(L,J)*TM1(L,K)
40 G(K,J) = G(J,K)
C PRINT MATRICES F AND G IF ISWP = 1
IF(ISWP.EQ.0) GO TO 45
CALL TITLE('MATRIX F',9)
CALL XPRINT(NRMF,F,NAT,NPAM,'E')
CALL TITLE('MATRIX G',9)
CALL XPRINT(NRMG,G,NPAM,NPAM,'E')
45 CONTINUE
C OBTAIN DISCRIMINANT SOLUTION FOR NPSW,NPEW PRINCIPAL AXES
DO 95 NPA = NPSW,NPEW
NPAAAA = NPA
WRITE(LUWN,903)NPA
903 FORMAT('0','DISCRIMINANT SOLUTION IN',I4,3X,'PRINCIPAL AXES')
C EIGENSOLUTION FOR SECTION OF G
IF(NPA.GT.1) GO TO 50
TM1(1,1) = 1.0D0
TM2(1,1) = G(1,1)
GO TO 55
50 NPAAAA=NPA
CALL XEIGEN(NRMG,NRMTM1,NRMTM2,G,TM1,TM2,NPAAAA,NPAAAA)
55 DO 60 J = 1,NPA
60 CV(J,3) = TM2(J,J)
C COMPUTATION OF LAMBDA AND BARLETT'S V STATISTIC
C ND = TOTAL NUMBER OF DISCRIMINANT DIMENSIONS
C RM = M PARAMETER
ND = NPA
IF(NDFH.LT.NPA) ND = NDFH
IRM = 2*(NDFH + NDFD) - (NPA + NDFH + 1)
RM = IRM/2.0D0
VTOT = 0.D0
DO 65 J = 1,ND
CV(J,4) = 1.0D0 - CV(J,3)
T1 = NDFD*CV(J,4)
T2 = RNT
IF(MOD.EQ.1) T2 = T2 - CV(J,4)
CV(J,5) = T2/T1
CV(J,4) = 1.0D0/SQRT(CV(J,4))
CV(J,6) = ALOG(CV(J,5))
CV(J,5) = CV(J,5) - 1.0D0
CV(J,7) = RM*CV(J,6)
65 VTOT = VTOT + CV(J,7)
C PRINT ND AND RM
WRITE(LUWN,904) ND
904 FORMAT('0','NUMBER OF DIMENSIONS, J, IN DISCRIMINANT SPACE',I4)
WRITE(LUWN,905)RM
905 FORMAT('0','PARAMETER M',F9.1)
C PRINT PHI AND LAMBDA
906 FORMAT('0',15X,'PHI',10X,'LAMBDA')
907 FORMAT('0',I6,2F14.5)
WRITE(LUWN,906)
DO 70 J = 1,ND
70 WRITE(LUWN,907)J,CV(J,3),CV(J,5)
C COMPUTE AND PRINT BARTLETT'S VT, DFT, PROBABILITY
WRITE(LUWN,908)
908 FORMAT('0','BARTLETT''S APPROXIMATE CHI-SQUARE STATISTIC')
WRITE(LUWN,909)
909 FORMAT('0','J = DISCRIMINANT; F = DISCRIMINANTS REMOVED')
WRITE(LUWN,910)
910 FORMAT('0','(NATURAL LOGARITHMS ARE USED)')
WRITE(LUWN,911)
911 FORMAT('0',5X,'J',3X,'LN(1 + LAMBDA)',6X,'V',6X,'DF',7X,'F',7X,
1'VT',5X,'DFT',3X,'PROBABILITY')
WRITE(LUWN,912)
IDFT = NPA*NDFH
IDF = NPA + NDFH - 1
DO 75 J = 1,ND
JM = J - 1
DFT = IDFT
PROB = CHIPRA(VTOT,DFT)
WRITE(LUWN,912)JM,VTOT,IDFT,PROB
912 FORMAT(' ',44X,I2,F11.3,I6,F12.5)
WRITE(LUWN,913)J,CV(J,6),CV(J,7),IDF
913 FORMAT(' ',I6,F14.5,F13.3,I5)
IDFT = IDFT - IDF
IDF = IDF - 2
75 VTOT = VTOT - CV(J,7)
C COMPUTE WEIGHTS
CALL XMATML(NRMF,NRMTM1,NRMTM2,F,TM1,TM2,NAT,NPAAAA,ND)
DO 78 K = 1,ND
T1 = 0.0D0
DO 76 J = 1,NAT
76 T1 = T1 + GPM(NGP,J)*TM2(J,K)
IF(T1.GE.0.0D0) GO TO 78
DO 77 J = 1,NAT
77 TM2(J,K) = -TM2(J,K)
78 CONTINUE
IF(ISWP.EQ.0) GO TO 80
CALL TITLE('WEIGHT MATRIX A',16)
CALL XPRINT(NRMTM2,TM2,NAT,ND,'E')
80 CONTINUE
DO 85 K = 1,ND
TEMP = SQRT(CV(K,4))
DO 85 J = 1,NAT
TM1(J,K) = TM2(J,K)*RTNT*CV(K,4)
85 TM2(J,K) = TM1(J,K)*CV(J,8)
CALL TITLE('WEIGHT MATRIX FROM RAW ATTRIBUTE MEASURES',41)
CALL TITLE('TO COMPOSITE SCORES STANDARDIZED IN TERMS OF DENOMINAT
1OR EFFECT',63)
CALL XPRINT(NRMTM1,TM1,NAT,ND,'F')
CALL TITLE('WEIGHT MATRIX FROM ATTRIBUTES STANDARDIZED IN TERMS OF
1 DENOMINATOR EFFECT',73)
CALL TITLE('TO COMPOSITE SCORES STANDARDIZED IN TERMS OF DENOMINAT
1OR EFFECT',63)
CALL XPRINT(NRMTM2,TM2,NAT,ND,'F')
C GROUP AND TOTAL MEANS
CALL XMATML(NRMGPM,NRMTM1,NRMTM2,GPM,TM1,TM2,NGP,NAT,ND)
DO 87 I = 1,NG
DO 87 J = 1,ND
87 TM2(I,J) = TM2(I,J) - TM2(NGP,J)
CALL TITLE('GROUP MEANS ON DISCRIMINANT COMPOSITES STANDARDIZED IN
1 TERMS OF DENOMINATOR EFFECT',82)
CALL XPRINT(NRMTM2,TM2,NG,ND,'F')
DO 90 J = 1,NAT
90 TM2(1,J) = TM2(NGP,J)
CALL TITLE('TOTAL MEANS ON DISCRIMINANT COMPOSITES STANDARDIZED IN
1 TERMS OF DENOMINATOR EFFECT',82)
CALL XPRINT(NRMTM2,TM2,1,ND,'F')
95 CONTINUE
IZINCR = IZINCR - 1
RETURN
END
C SUBROUTINE MISNR (XDA,XBA,XVA,CXY,RXY,XMI,T,TV,XN,ZXY,
C 2 NS,NV,IOPT,IONE)
C TITLE: MISNR
C PROGRAMMER: DR. ALLEN I. FLEISHMAN
C DATE: DECEMBER, 1978
C LOCATION: LEDERLE LABS
C PURPOSE: TO COMPUTE CORRELATION MATRIX AND VECTOR OF MEANS
C AND VARIANCES, EVEN IF MISSING DATA ARE PRESENT.
C OPTIONS:
C
C 1 COMPUTE VECTOR OF MEANS AND SAMPLE SIZES ONLY
C (OUTPUT TO DIAGONAL OF XBA). ALL OTHER MATRICES ZEROED.
C 2 NO MISSING DATA, WILL NOT CHECK (MOST EFFICIENT
C IF THERE IS NO MISSING DATA; WILL PRODUCE AN
C ERROR IF ANY DATA IS MISSING).
C 3 REPLACE MISSING DATA BY MEAN.
C 4 DELETION OF ALL DATA FOR ANY SUBJECT WHO HAS
C ANY MISSING DATA PRESENT.
C 5 COMPUTE COVARIANCES BASED ON AVAILABLE PAIRS OF
C DATA. MEANS AND COVARIANCES WILL BE BASED ON
C AVAILABLE PAIRWISE DATA FOR THAT PAIR OF VARIABLES
C ONLY. CORRELATIONS WILL BE BASED ON THESE COVARIANCES
C AND VARIANCES FROM DIAGONAL OF COV. MATRIX.
C (DEFAULT)
C 6 COMPUTE PAIRWISE COVARIANCE AS IN OPTION 5 ABOVE.
C ESTIMATES FOR VARIANCES ARE BASED ON AVAILABLE
C PAIRWISE DATA FOR THAT PAIR OF VARIABLES ONLY.
C 7 COVARIANCES BASED ON AVAILABLE PAIRS, HOWEVER THE
C FULL DATA MEANS WILL BE UTILIZED. CORRELATIONS
C WILL BE BASED ON FULL DATA VARIANCES (DIAG OF COV.
C MATRIX).
C 8 COVARIANCE MATRIX COMPUTED AS IN OPTION 7 ABOVE.
C VARIANCE FOR CORRELATION MATRIX COMPUTED FROM AVAILABLE
C PAIRS OF DATA WHEN FULL DATA MEAN IS UTILIZED.
C
C VARIABLES USED
C XDA RECTANGULAR MATRIX CONTAINING RAW DATA (DIMENSIONS
C INPUT NS - NUMBER OF SUBJECTS X NV - NUMBER OF VARIABLES)
C XBA MATRIX OF MEANS OF VARIABLES. IF OPTIONS 5 OR 6 ARE
C OUTPUT USED THIS MATRIX WILL BE MEAN FOR ROW VARIABLE FOR
C ALL CASES WHEN THERE WAS DATA FOR THIS VARIABLE AND
C COLUMN VARIABLE.
C XVA MATRIX OF VARIANCES OF VARIABLES. IF OPTIONS 6 OR 8 ARE
C OUTPUT USED THIS MATRIX WILL BE VARIANCE FOR COLUMN VARIABLE FOR
C ALL CASES WHEN THERE WAS DATA FOR THIS VARIABLE AND
C ROW VARIABLE.
C CXY MATRIX OF COVARIANCES.
C OUTPUT
C RXY MATRIX OF CORRELATIONS.
C OUTPUT
C XMI VECTOR OF LENGTH NV CONTAINS INFORMATION AS TO HOW MISSING
C INPUT DATA WERE CODED (PER VARIABLE). MISSING DATA FOR ANY
C VARIABLE IS IDENTIFIED BY CODING IT AS IN TV(I).
C THIS SUBROUTINE CAN NOT DIFFERENTIATE BETWEEN -0.0 (OR
C BLANKS) AND 0.0 (OR ZERO).
C IF THIS PROCEDURE CANNOT DESCRIBE HOW MISSING
C DATA IS CODED IN THE USERS STUDY, WRITE
C A NEW SUBROUTINE SUBMIS (SEE THAT SUBROUTINE
C FOR ITS PARAMETERS).
C T A DUMMY DOUBLE PRECISION VECTOR OF LENGTH NV.
C TV A DUMMY DOUBLE PRECISION VECTOR OF LENGTH NV.
C XN A DUMMY MATRIX OF SIZE (NV X NV).
C ZXY MATRIX OF SIZE NV X NV CONTAINING SAMPLE SIZE FOR
C OUTPUT PAIRWISE COVARIANCES.
C NS NUMBER OF SUBJECTS
C INPUT
C NV NUMBER OF VARIABLES
C INPUT
C IOPT OPTION ON HOW TO HANDLE MISSING DATA (SEE ABOVE)
C INPUT
C IONE DENOMINATOR OF VARIANCES AND COVARIANCES.
C INPUT IF IONE: -1 THEN DENOMINATOR IS N - 1 (UNBIASED),
C 0 THEN DENOMINATOR IS N (MAXIMUM LIKELIHOOD), OR
C 1 THEN DENOMINATOR IS N + 1 (MINIMUM ERROR).
C
C
C KMISNR (SET IN COMMON BLOCK C) SWITCH FOR MEAN SUBTRACTION:
C INPUT IF KMISNR = 1, THE MEAN WILL NOT BE SUBTRACTED FROM
C ANY OF THE COVARIANCES (MAKING THEM
C CROSS PRODUCTS) NOR FROM THE CORREL-
C ATIONS [MAKING THEM STANDARDIZED OR
C NORMED LENGTHS (CROSS PRODUCTS)].
C = 0 (DEFAULT) THE MEAN WILL BE APPROPRIATELY
C SUBTRACTED AS USUAL.
C
SUBROUTINE MISNR (XDA,XBA,XVA,CXY,RXY,XMI,T,TV,XN,ZXY,
2 NS,NV,IOPT,IONE)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION XDA(IO,IO),XBA(IO,IO),XVA(IO,IO),CXY(IO,IO),RXY(IO,IO)
DIMENSION XMI(IO,IO),T(IO,IO)
DIMENSION TV(IO,IO),XN(IO,IO),ZXY(IO,IO)
COMMON /B/IO
COMMON /ZINCR/ IZINCR
CALL QINCR ('MISNR')
NRMXDA=IO
NRMXBA=IO
NRMXVA=IO
NRMCXY=IO
NRMRXY=IO
NRMXMI=IO
NRMT=IO
NRMTV=IO
NRMXN=IO
NRMZXY=IO
CALL XMISNR(NRMXDA,NRMXBA,NRMXVA,NRMCXY,NRMRXY,NRMXMI,NRMT,
2 NRMTV,NRMXN,NRMZXY,XDA,XBA,XVA,CXY,RXY,XMI,T,TV,XN,ZXY,NS,
3 NV,IOPT,IONE)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XMISNR(NRMXDA,NRMXBA,NRMXVA,NRMCXY,NRMRXY,NRMXM1,NRMT,
2NRMTV,NRMXN,NRMZXY,XDAT,XBAR,XVAR,CXY,RXY,XMISS,T,TV,XN,
3ZXY,NS,NV,IOPT,IONE)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION XDAT(NRMXDA,1),XBAR(NRMXBA,1),XVAR(NRMXVA,1)
DIMENSION CXY(NRMCXY,1),RXY(NRMRXY,1),XMISS(NRMXM1)
DIMENSION T(NRMT),TV(NRMTV),XN(NRMXN,1),ZXY(NRMZXY,1)
COMMON /B/IO
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
COMMON/KMISNR/KMISNR
CALL QINCR ('XMISNR')
CALL XDIMCH(NRMXDA,NS,'MISNR',0)
CALL XDIMCH(NRMXBA,NV,'MISNR',0)
CALL XDIMCH(NRMXVA,NV,'MISNR',0)
CALL XDIMCH(NRMCXY,NV,'MISNR',0)
CALL XDIMCH(NRMRXY,NV,'MISNR',0)
CALL XDIMCH(NRMXM1,NV,'MISNR',0)
CALL XDIMCH(NRMT,NV,'MISNR',0)
CALL XDIMCH(NRMTV,NV,'MISNR',0)
CALL XDIMCH(NRMXN,NV,'MISNR',0)
CALL XDIMCH(NRMZXY,NV,'MISNR',0)
ACC = 10.0D-5
KMISS = 1 - KMISNR
IF (KMISS.LT.0.OR.KMISS.GT.1) KMISS = 1
C ACC IS THE ACCURACY OF THE MACHINE
C CHECK IOPT AND RESET TO DEFAULT (IOPT=7) IF OUT OF BOUNDS.
IOP=IOPT
IF (IOP.LT.1) IOP=7
IF (IOP.GT.8) IOP=7
C CHECK IONE AND RESET TO DEFAULT (IONE=0) IF OUT OF BOUNDS
ION=IONE
IF (ION.LT.-1) ION=0
IF (ION.GT.+1) ION=0
C ZERO MATRICES
Z = 0.0D0
CALL XGENR8 (NRMCXY,Z,CXY,NV,NV)
CALL XGENR8 (NRMRXY,Z,RXY,NV,NV)
CALL XGENR8 (NRMXBA,Z,XBAR,NV,NV)
CALL XGENR8 (NRMXVA,Z,XVAR,NV,NV)
CALL XGENR8 (NRMZXY,Z,ZXY,NV,NV)
DO 30 I = 1,NS
IF (IOP.NE.4) GOTO 33
DO 32 J = 1,NV
32 IF (DABS(XDAT(I,J) - XMISS(J)).LT.ACC) GOTO 30
33 DO 30 J = 1,NV
IF (IOP.NE.2.AND.DABS(XDAT(I,J)-XMISS(J)).LT.ACC) GOTO 30
XBAR(J,1)=XBAR(J,1)+XDAT(I,J)
ZXY(J,1)=ZXY(J,1) + 1.0D0
30 CONTINUE
DO 31 I = 1,NV
XBAR(I,1) = XBAR(I,1)/ZXY(I,1)
31 TV(I) = XBAR(I,1)
DO 50 I = 2,NV
IF (IOP.EQ.1) ZXY(I,I) = ZXY(I,1)
XBAR(I,I) = XBAR(I,1)
XBAR(I,1)=Z
50 ZXY(I,1) = 0.0D0
IF (IOP.NE.1) ZXY(1,1) = 0.0D0
IF (IOP.EQ.1) GOTO 9999
IF (IOP.GT.4) GOTO 500
C
200 CONTINUE
NSS = 0
DO 201 I = 1,NS
DO 202 J = 1,NV
T(J) = XDAT(I,J)
IF (IOP.EQ.2) GOTO 202
IF (IOP.EQ.3.AND.(DABS(T(J)-XMISS(J)).LT.ACC)) T(J) = XBAR(J,J)
IF (IOP.EQ.4.AND.(DABS(T(J)-XMISS(J)).LT.ACC)) GOTO 201
202 T(J) = T(J) - XBAR(J,J)*KMISS
NSS = NSS + 1
DO 203 K = 1,NV
DO 203 J = K,NV
203 CXY(K,J) = CXY(K,J) + T(K)*T(J)
201 CONTINUE
DO 205 I = 1,NV
DO 204 J = I,NV
CXY(I,J) = CXY(I,J)/(NSS + ION)
204 CXY(J,I) = CXY(I,J)
ZXY(I,I) = NSS
205 XVAR(I,I)=CXY(I,I)
C
GOTO (211,250,250,250,211), IOP
211 CALL END (1)
WRITE (LUWN,221) IOP
221 FORMAT ('0 SOMETHING TERRIBLE IS WRONG AT LINE 211 IN PROGRAM',
1' MISNR. IOPT WAS ', I20///'0BRING THIS OUTPUT TO DR. ',
2'FLEISHMAN'////)
CALL END (2)
CALL END (3)
250 CONTINUE
260 CONTINUE
DO 262 I = 1,NV
DO 262 J = 1,NV
262 XVAR(I,J)=CXY(I,I)
263 CONTINUE
INOTE = 0
DO 261 I = 1,NV
IF (XVAR(I,I).LT.ACC) WRITE (LUWN,290) I
290 FORMAT (' WARNING: VARIABLE ', I3,' DOES NOT VARY')
DO 261 J = I,NV
IF (XVAR(I,J)*XVAR(J,I).GT.ACC) GOTO 265
RXY(I,J) = Z
RXY(J,I) = Z
GOTO 267
265 RXY(I,J) = CXY(I,J)/DSQRT(XVAR(I,J)*XVAR(J,I))
RXY(J,I) = RXY(I,J)
267 IF (I.EQ.J.AND.RXY(I,I).LT.0.999D0) RXY(I,I) = 1.0D0
IF (RXY(I,J).GT.1.0D0) INOTE = 1
261 CONTINUE
IF (INOTE.GT.0) WRITE (LUWN,280)
280 FORMAT (' WARNING: ONE OF YOUR CORRELATIONS IS GREATER THAN ONE')
GOTO 9999
C
500 CONTINUE
DO 504 I = 1,NV
504 XBAR(I,I)=Z
DO 503 I = 1,NS
DO 502 J = 1,NV
502 T(J) = XDAT(I,J)
CALL SUBMIS(NRMXN,T,ACC,XMISS,XN,NV)
DO 505 K = 1,NV
505 T(K) = T(K) - TV(K)*KMISS
DO 503 J = 1,NV
DO 503 K = 1,NV
IF (XN(J,K).LT.0.5D0) GOTO 503
IF (IOP.GE.7) GOTO 507
XBAR(J,K)=XBAR(J,K) + T(J)
507 ZXY(J,K) = ZXY(J,K) + 1.0D0
CXY(J,K) = CXY(J,K) + T(J)*T(K)
IF (IOP.EQ.5.OR.IOP.EQ.7) GOTO 503
XVAR(J,K) = XVAR(J,K) + (T(J))**2
503 CONTINUE
DO 520 I = 1,NV
DO 520 J = 1,NV
IF (IOP.GE.7) GOTO 521
IF (IOP.EQ.5) GOTO 524
XVAR(I,J)=XVAR(I,J) - ((XBAR(I,J)**2)*KMISS/ZXY(I,J))
524 CXY(I,J)=CXY(I,J)-((XBAR(I,J)*XBAR(J,I)*KMISS)/ZXY(I,J))
XBAR(I,J)=XBAR(I,J)/ZXY(I,J) + TV(I)*KMISS
GOTO 523
521 XBAR(I,J) = TV(I)
523 XVAR(I,J) = XVAR(I,J)/(ZXY(I,J)+ION)
CXY (I,J) = CXY (I,J)/(ZXY(I,J)+ION)
520 CONTINUE
IF (IOP.EQ.5.OR.IOP.EQ.7) GOTO 260
GOTO 263
9999 IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE SUBMIS (NRMXN,T,ACC,XMISS,XN,NV)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION T(1), XMISS(1), XN(NRMXN,1)
C
C TITLE: SUBJECT'S MISSING DATA (SUBROUTINE OF MISNR)
C PURPOSE: TO INDICATE WHICH DATA POINTS ARE MISSING
C THIS SUBROUTINE READS IN A SINGLE SUBJECT'S DATA - IN A
C VECTOR CALLED T, SEES WHETHER IT IS EQUAL TO XMISS - A
C VECTOR TELLING THE PROGRAM HOW MISSING DATA WAS CODED -
C THEN RETURNS A VALUE OF 1.0 IN MATRIX XN(I,J) IF THE DATA
C FOR THAT PAIR OF VARIABLES IS GOOD, OR 0.0 IF THAT PAIR
C HAD MISSING DATA.
C
C SYMBOLS:
C
C T = VECTOR OF DATA FOR A SINGLE SUBJECT
C NRMXN = NUMBER OF ROWS DIMENSIONED FOR XN IN MAIN
C ACC = HOW SMALL A DIFFERENCE BETWEEN OBSERVED AND MISSING
C VALUE IS BEFORE CONSIDERED MISSING
C (SET IN MAIN TO BE [0.00001]). RELATED TO
C ACCURACY OF COMPUTER.
C XMISS = VECTOR CONTAINING ELEMENTS INDICATING HOW MISSING DATA
C WERE CODED PER EACH VARIABLE.
C XN = MATRIX CONTAINING ELEMENTS OF EITHER ZERO (IGNORE THIS
C PAIR FOR THIS SUBJECT) OR ONE (RETAIN DATA FOR THIS
C PAIR FOR THIS SUBJECT).
C NV = NUMBER OF VARIABLES CORRELATED.
C
C NO CALL TO QINCR WILL BE MADE.
C
DO 1 I = 1,NV
XN(I,I) = 1.0D0
1 IF (DABS(T(I) - XMISS(I)).LT.ACC) XN(I,I) = 0.0D0
N=NV-1
DO 2 I = 1,N
K=I+1
DO 2 J = K,NV
XN(I,J)=0.0D0
IF (XN(I,I)*(XN(J,J)).GT.0.5D0) XN(I,J)=1.0D0
2 XN(J,I)=XN(I,J)
RETURN
END
SUBROUTINE RNDCR(CS,RS,CP,F,NIN,NAT,NCF,IRN,MF)
C GENERATES A COVARIANCE MATRIX AND CORRELATION MATRIX AMONG NAT ATTRIBUTES
C FOR A SAMPLE SIZE NIN DRAWN FROM A POPULATION DISTRIBUTED MULTIDIMENSIONAL
C NORMAL WITH GIVEN COVARIANCE MATRIX AND MEAN VECTOR
C PROGRAM WRITTEN BY LEDYARD R TUCKER, APRIL 1978
C MATRICES; ROW DIMENSION OF ARRAYS TO CONTAIN MATRICES IS I1 IN
C COMMON BLOCK B , MUST BE AT LEAST NAT+2
C COLUMN DIMENSION OF ARRAYS TO CONTAIN MATRICES MUST BE AT LEAST NAT
C CS CONTAINS SAMPLE COVARIANCE MATRIX WITH SAMPLE MEANS IN ROW NAT+1
C AND SAMPLE STANDARD DEVIATIONS IN ROW NAT+2; OUTPUT
C RS CONTAINS SAMPLE CORRELATION MATRIX; OUTPUT
C CP CONTAINS POPULATION COVARIANCE MATRIX IN FIRST NAT ROWS AND
C POPULATION MEAN VECTOR IN ROW NAT+1 ; INPUT/OUTPUT
C F CONTANINS FACTOR MATRIX OF POPULATION COVARIANCE MATRIX;
C MAY BE INPUT, OTHERWISE IS OBTAINED BY A CHOLESKY
C DECOMPOSITION OF CP ; OUTPUT
C SCALARS:
C NRMCS IS THE NUMBER OF ROWS OF CS IN MAIN = I1
C NRMRS IS THE NUMBER OF ROWS OF RS IN MAIN = I1
C NRMCP IS THE NUMBER OF ROWS OF CP IN MAIN = I1
C NRMF IS THE NUMBER OF ROWS OF F IN MAIN = I1
C NIN IS SAMPLE SIZE; INPUT/OUTPUT
C NAT IS NUMBER OF ATTRIBUTES; INPUT/OUTPUT
C NCF IS NUMBER OF FACTORS IN F ; INPUT WHEN F IS INPUT; OUTPUT
C IRN IS A 9 DIGIT INTEGER USED IN GENERATING RANDOM
C VALUES; INPUT, CHANGED BEFORE OUTPUT
C RAN IS THE SUBROUTINE WHICH FURNISHES THE RANDOM NUMBERS.
C SETRAN 'SETS' THE SEED FOR RAN, AND
C SAVRAN RETREIVES THAT SEED FOR FUTURE USE.
C MF IS A FLAG CONCERNING MATRIX F ; INPUT:
C = 0 WHEN F IS NOT INPUT
C = 1 WHEN F IS INPUT,
C ON OUTPUT MF = 1
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION CS(I1,1),RS(I1,1),CP(I1,1),F(I1,1)
COMMON/B/I1
COMMON /ZINCR/ IZINCR
CALL QINCR ('RNDCR')
NRMCS = I1
NRMRS = I1
NRMCP = I1
NRMF = I1
CALL XRNDCR(NRMCS,NRMRS,NRMCP,NRMF,CS,RS,CP,
2 F,NIN,NAT,NCF,IRN,MF)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XRNDCR(NRMCS,NRMRS,NRMCP,NRMF,CS,RS,CP,
2 F,NIN,NAT,NCF,IRN,MF)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION CS(NRMCS,1),RS(NRMRS,1),CP(NRMCP,1),F(NRMF,1)
COMMON/ZINCR/IZINCR
CALL QINCR ('XRNDCR')
C OBTAIN SAMPLE Q MATRIX
CALL XRANDQ (NRMCS,NRMCP,NRMF,NRMRS,CS,CP,F,RS,MIN,
2NAT,NCF,IRN,MF)
C CONVERT Q TO SAMPLE C
NM1 = NIN - 1
NATP2 = NAT + 2
NATP1 = NAT + 1
CALL XDIMCH (NRMCS,NATP2,'RNDCR',0)
CALL XDIMCH (NRMCP,NATP1,'RNDCR',0)
CALL XDIMCH (NRMRS,NAT,'RNDCR',0)
CALL XDIMCH(NRMF,NAT,'RNDCR',0)
DO 5 I = 1,NAT
DO 5 J = 1,I
CS(I,J) = CS(I,J)/NM1
5 CS(J,I) = CS(I,J)
C OBTAIN STANDARD DEVIATIONS AND CORRELATION MATRIX
DO 10 I = 1,NAT
CS(NATP2,I) = SQRT(CS(I,I))
DO 10 J = 1,I
T = CS(NATP2,I)*CS(NATP2,J)
RS(I,J) = CS(I,J)/T
10 RS(J,I) = RS(I,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE RANDQ(Q,CP,F,S,NIN,NAT,NCF,IRN,MF)
C GENERATES A SAMPLE SSCP MATRIX AND MEAN VECTOR
C THE SSCP MATRIX CONTAINS SUMS OF SQUARES AND CROSS PRODUCTS OF
C MEASURES DEVIATED FROM SAMPLE MEANS
C THE POPULATION IS MULTIDIMENSIONAL NORMAL WITH A GIVEN MEAN VECTOR
C AND COVARIANCE MATRIX
C PROGRAM WRITTEN BY LEDYARD R TUCKER, MARCH 1978
C MATRICES; ROW DIMENSION OF ARRAYS TO CONTAIN MATRICES IS I1 IN
C COMMON BLOCK B , MUST BE AT LEAST NAT+1
C COLUMN DIMENSION OF ARRAYS TO CONTAIN MATRICES MUST BE AT LEAST
C NAT
C Q CONTAINS SAMPLE SSCP MATRIX IN FIRST NAT ROWS AND
C SAMPLE MEAN VECTOR IN ROW NAT+1 ; OUTPUT
C CP CONTAINS POPULATION COVARIANCE MATRIX IN FIRST NAT ROWS AND
C POPULATION MEAN VECTOR IN ROW NAT+1 ; INPUT/OUTPUT
C F CONTANINS FACTOR MATRIX OF POPULATION COVARIANCE MATRIX;
C MAY BE INPUT, OTHERWISE IS OBTAINED BY A CHOLESKY
C DECOMPOSITION OF CP ; OUTPUT
C S IS A SCRATCH MATRIX
C SCALARS:
C NRMQ IS THE NUMBER OF ROWS OF Q IN MAIN = I1
C NRMCP IS THE NUMBER OF ROWS OF CP IN MAIN = I1
C NRMF IS THE NUMBER OF ROWS OF F IN MAIN = I1
C NRMS IS THE NUMBER OF ROWS OF S IN MAIN = I1
C NIN IS SAMPLE SIZE; INPUT/OUTPUT
C NAT IS NUMBER OF ATTRIBUTES; INPUT/OUTPUT
C NCF IS NUMBER OF FACTORS IN F ; INPUT WHEN F IS INPUT; OUTPUT
C IRN IS A 15 DIGIT INTEGER ENDING IN 1 USED IN GENERATING RANDOM
C VALUES; INPUT, CHANGED BEFORE OUTPUT
C MF IS A FLAG CONCERNING MATRIX F ; INPUT:
C = 0 WHEN F IS NOT INPUT
C = 1 WHEN F IS INPUT,
C ON OUTPUT MF = 1
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION Q(I1,1),CP(I1,1),F(I1,1),S(I1,1)
COMMON/B/I1
COMMON /ZINCR/IZINCR
NRMQ=I1
NRMCP=I1
NRMF=I1
NRMS=I1
CALL QINCR ('RANDQ')
CALL XRANDQ(NRMQ,NRMCP,NRMF,NRMS,Q,CP,F,S,NIN,NAT,NCF,
2 IRN,MF)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XRANDQ(NRMQ,NRMCP,NRMF,NRMS,Q,CP,F,S,NIN,NAT,NCF,
2 IRN,MF)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION Q(NRMQ,1),CP(NRMCP,1),F(NRMF,1),S(NRMS,1)
COMMON/ZINCR/IZINCR
CALL QINCR ('XRANDQ')
NATP1 = NAT + 1
CALL XDIMCH (NRMQ ,NATP1,'RANDQ',0)
CALL XDIMCH (NRMS ,NATP1,'RANDQ',0)
CALL XDIMCH (NRMCP,NATP1,'RANDQ',0)
CALL XDIMCH (NRMF ,NATP1,'RANDQ',0)
C COMPUTE F WHEN NECESSARY
IF(MF.NE.0) GO TO 10
CALL XCHLSK(NRMCP,NRMF,CP,F,NAT)
NCF = NAT
MF = 1
10 CONTINUE
C OBTAIN Q AND M FOR STANDARD POPULATION
CALL XRNDQM(NRMQ,Q,NIN,NCF,IRN)
C TRANSFORM Q AND M
NCFP1 = NCF + 1
DO 15 I = 1,NCFP1
DO 15 J = 1,NAT
S(I,J) = 0.0D0
DO 15 K = 1,NCF
15 S(I,J) = S(I,J) + Q(I,K)*F(J,K)
DO 30 I = 1,NAT
DO 25 J = 1,I
Q(I,J) = 0.0D0
DO 20 K = 1,NCF
20 Q(I,J) = Q(I,J) + F(I,K)*S(K,J)
25 Q(J,I) = Q(I,J)
30 Q(NATP1,I) = S(NCFP1,I) + CP(NATP1,I)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE RNDQM(Q,NIN,NAT,IRN)
C GENERATES A SAMPLE SSCP MATRIX AND MEAN VECTOR
C THE SSCP MATRIX CONTAINS SUMS OF SQUARES AND CROSS PRODUCTS OF
C MEASURES DEVIATED FROM SAMPLE MEANS
C POPULATION IS MULTIDIMENSIONAL NORMAL WITH MEAN VECTOR = 0 AND
C COVARIANCE MATRIX = I
C PROGRAM WRITTEN BY LEDYARD R TUCKER, MARCH 1978
C Q CONTAINS SAMPLE SSCP MATRIX IN FIRST NAT ROWS AND
C SAMPLE MEAN VECTOR IN ROW NAT+1 ; OUTPUT
C ROW DIMENSION OF ARRAY Q IS I1 IN COMMON BLOCK B
C ROW DIMENSION OF ARRAY Q MUST BE AT LEAST NAT+1
C COLUMN DIMENSION OF ARRAY Q MUST BE AT LEAST NAT
C NIN IS SAMPLE SIZE; INPUT/OUTPUT
C NAT IS NUMBER OF ATTRIBUTES; INPUT/OUTPUT
C IRN IS A 15 DIGIT INTEGER ENDING IN 1 USED IN GENERATING RANDOM
C VALUES; INPUT, CHANGED BEFORE OUTPUT
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION Q(I1,1)
COMMON/B/I1
COMMON /ZINCR/ IZINCR
CALL QINCR ('RNDQM')
NRMQ = I1
CALL XRNDQM(NRMQ,Q,NIN,NAT,IRN)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XRNDQM(NRMQ,Q,NIN,NAT,IRN)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION Q(NRMQ,1)
COMMON/ZINCR/IZINCR
CALL QINCR ('XRNDQM')
CALL SETRAN(IRN)
NATP1 = NAT + 1
CALL XDIMCH (NRMQ,NATP1,'RNDQM',0)
CALL XGENR8(NRMQ,0.0D0,Q,NRMQ,NRMQ)
C PROTECT AGAINST NIN.LT.1
IF(NIN.LT.1) IZINCR = IZINCR - 1
IF(NIN.LT.1) RETURN
C CASE WHEN NIN.EQ.1
IF(NIN.EQ.1) GO TO 47
C DETERMINE ND
ND = NIN - 1
IF(ND.GT.NAT) ND = NAT
C INSERT RANDOM NORMAL DEVIATES BELOW DIAGONAL IN ND COLUMNS OF H
DO 10 J = 1,ND
IF(J.GE.NAT) GO TO 12
JP = J + 1
DO 10 I = JP,NAT
10 Q(I,J) = QRAND1(T1)
12 CONTINUE
C INSERT RANDOM CHI IN DIAGONAL OF H , DF = NIN - J
DO 25 J = 1,ND
NDF = NIN - J
IF(NDF.GT.5) GO TO 22
T2 = 0.0D0
DO 21 I = 1,NDF
T1 = QRAND1(T3)
21 T2 = T2 + T1**2
GO TO 25
22 T1 = RAN(T2)
T2 = CHIVLT(T1,NDF)
25 Q(J,J) = SQRT(T2)
C SAMPLE SSCP MATRIX
DO 45 J = 1,NAT
JB = NATP1 - J
DO 40 I = 1,JB
T2 = 0.0D0
DO 35 K = 1,I
35 T2 = T2 + Q(I,K)*Q(JB,K)
40 Q(I,JB) = T2
DO 45 I = 1,JB
45 Q(JB,I) = Q(I,JB)
C SAMPLE MEANS
47 CONTINUE
T1 = NIN
T1 = SQRT(T1)
DO 50 J = 1,NAT
50 Q(NATP1,J) = QRAND1(T2)/T1
CALL SAVRAN(IRN)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE TSAMX(X,C,F,NIN,NAT,NCF,IRN,MF)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(IO,IO), C(IO,IO), F(IO,IO)
COMMON /B/IO
COMMON /ZINCR/IZINCR
CALL QINCR ('TSAMX')
NRMX = IO
NRMC = IO
NRMF = IO
CALL XTSAMX(NRMX,NRMC,NRMF,X,C,F,NIN,NAT,NCF,IRN,MF)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XTSAMX(NRMX,NRMC,NRMF,X,C,F,NIN,NAT,NCF,IRN,MF)
C OBTAINS SAMPLE DATA MATRIX FROM MULTIDIMENSIONAL NORMAL POPULATION
C HAVING GIVEN POPULATION MEAN VECTOR AND COVARIANCE MATRIX
C PROGRAM WRITTEN BY LEDYARD R TUCKER, MARCH 1978
C MATRICES:
C X IS SAMPLE DATA MATRIX ON OUTPUT
C ARRAY TO CONTAIN X MUST HAVE AT LEAST NIN ROWS AND NAT COLUMNS
C C IS POPULATION COVARIANCE MATRIX WITH MEAN VECTOR IN ROW NAT+1;
C INPUT/OUTPUT
C ARRAY TO CONTAIN C MUST HAVE AT LEAST NAT+1 ROWS AND COLUMNS
C F IS FACTOR MATRIX OF POPULATION COVARIANCE MATRIX, MAY BE INPUT,
C OTHERWISE IS OBTAINED BY A CHOLESKY DECOMPOSITION; OUTPUT
C ARRAY TO CONTAIN F MUST HAVE AT LEAST NAT+1 ROWS AND COLUMNS
C SCALARS:
C NIN IS NUMBER OF INDIVIDUALS IN SAMPLE; INPUT/OUTPUT
C NAT IS NUMBER OF ATTRIBUTES; INPUT/OUTPUT
C NCF IS NUMBER OF COLUMNS OF F ; INPUT IF F IS INPUT; OUTPUT
C NRMX, NRMC, NRMF ARE NUMBER OF ROWS IN MAIN OF ARRAYS TO CONTAIN
C MATRICES X , C , F ; INPUT/OUTPUT
C IRN IS A 15 DIGIT INTEGER ENDING IN 1; INPUT, ALTERED BEFORE OUTPUT
C MF IS A FLAG CONCERNING MATRIX F ; INPUT:
C = 0 IF F IS NOT INPUT,
C = 1 IF F IS INPUT
C IS SET TO 1 FOR OUTPUT
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION X(NRMX,1),C(NRMC,1),F(NRMF,1)
COMMON /ZINCR/ IZINCR
CALL QINCR ('XTSAMX')
NATP1 = NAT + 1
CALL XDIMCH(NRMX,NIN,'XTSAMX',0)
CALL XDIMCH(NRMC,NATP1,'XTSAMX',0)
CALL XDIMCH(NRMF,NATP1,'XTSAMX',0)
CALL SETRAN(IRN)
C COMPUTE F WHEN NECESSARY
IF(MF.NE.0) GO TO 10
CALL XCHLSK(NRMC,NRMF,C,F,NAT)
NCF=NAT
MF = 1
10 CONTINUE
C GENERATE DATA ONE ROW AT A TIME
DO 20 IN = 1,NIN
C OBTAIN VECTOR OF NORMAL RANDOM DEVIATES IN ROW NAT+1 OF F
DO 15 J = 1,NAT
T = RAN(T)
15 F(NATP1,J) = QDVNRM(T)
C TRANSFORM TO ROW OF X
DO 20 K = 1,NAT
X(IN,K) = C(NATP1,K)
DO 20 J = 1,NAT
20 X(IN,K) = X(IN,K) + F(NATP1,J)*F(K,J)
CALL SAVRAN(IRN)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE FAWMV(S1,S2,S3,S4,NAT,NCF)
C FACTOR ANALYSIS OF C WITH MV IN DIAGONALS
C DRAFT BY LRT, MARCH, 1977, REVISED APRIL 1978
C RUNS ON FORMAL
C COVARIANCE MATRIX C IS INPUT/OUTPUT IN S1
C DIMENSIONALITY MUST BE LESS THAN I1
C S2, S3, S4 ARE SCRATCH ON INPUT
C S2 ON OUTPUT CONTAINS C WITH MV IN DIAGONAL
C NOTE: IF C IS NOT DESIRED ON OUTPUT, S1 AND S2 CAN BE IDENTICAL
C S3 ON OUTPUT CONTAINS THE FACTOR MATRIX
C S4 ON OUTPUT CONTAINS COMMUNALITIES FOR EACH GIVEN NUMBER OF FACTORS
C NAT IS NUMBER OF ATTRIBUTES; INPUT/OUTPUT
C NCF IS MAXIMUM NUMBER OF FACTORS TO BE OBTAINED; IF NCF IS GREATER THAN
C THE NUMBER OF POSITIVE EIGENVALUES OF C WITH MULTIPLE VARIANCES IN THE
C DIAGONAL CELLS, NCF WILL BE REDUCED TO THE NUMBER OF POSITIVE
C EIGENVALUES. NCF MUST BE A SYMBOLIC REPRESENTATION OF AN INTEGER.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION S1(I1,I1), S2(I1,I1), S3(I1,I1), S4(I1,I1)
COMMON/B/I1
COMMON /ZINCR/ IZINCR
CALL QINCR ('FAWMV')
NRMS1 = I1
NCMS1 = I1
NRMS2 = I1
NRMS3 = I1
NCMS3 = I1
NRMS4 = I1
CALL XFAWMV(NRMS1,NCMS1,NRMS2,NRMS3,NCMS3,NRMS4,S1,S2,S3,S4,
2 NAT,NCF)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE FAHIR(S1,S2,S3,S4,NAT,NCF)
C FACTOR ANALYSIS OF C WITH HIGHEST COVARIANCE IN DIAGONALS
C BY ALLEN FLEISHMAN, 1979
C RUNS ON FORMAL
C COVARIANCE MATRIX C IS INPUT/OUTPUT IN S1
C DIMENSIONALITY MUST BE LESS THAN I1
C S2, S3, S4 ARE SCRATCH ON INPUT
C S2 ON OUTPUT CONTAINS C WITH HIGHEST COVARIANCE IN DIAGONAL
C NOTE: IF C IS NOT DESIRED ON OUTPUT, S1 AND S2 CAN BE IDENTICAL
C S3 ON OUTPUT CONTAINS THE FACTOR MATRIX
C S4 ON OUTPUT CONTAINS COMMUNALITIES FOR EACH GIVEN NUMBER OF FACTORS
C NAT IS NUMBER OF ATTRIBUTES; INPUT/OUTPUT
C NCF IS MAXIMUM NUMBER OF FACTORS TO BE OBTAINED; IF NCF IS GREATER THAN
C THE NUMBER OF POSITIVE EIGENVALUES OF C WITH MULTIPLE VARIANCES IN THE
C DIAGONAL CELLS, NCF WILL BE REDUCED TO THE NUMBER OF POSITIVE
C EIGENVALUES. NCF MUST BE A SYMBOLIC REPRESENTATION OF AN INTEGER.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION S1(I1,I1), S2(I1,I1), S3(I1,I1), S4(I1,I1)
COMMON/B/I1
COMMON /ZINCR/ IZINCR
CALL QINCR ('FAHIR')
NRMS1 = I1
NCMS1 = I1
NRMS2 = I1
NRMS3 = I1
NCMS3 = I1
NRMS4 = I1
CALL XFAHIR (NRMS1,NCMS1,NRMS2,NRMS3,NCMS3,NRMS4,S1,S2,S3,S4,
2 NAT,NCF)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XFAWMV(NRMS1,NCMS1,NRMS2,NRMS3,NCMS3,NRMS4,S1,S2,
2 S3,S4,NAT,NCF)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION S1(NRMS1,1),S2(NRMS2,1),S3(NRMS3,1),S4(NRMS4,1)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
CALL QINCR ('XFAWMV')
IF (NAT.LT.4) CALL END (1)
IF (NAT.LT.4) WRITE (LUWN, 77) NAT
77 FORMAT (40X,'A FACTOR ANALYSIS WAS DESIRED WITH ',I10,' VARIABLES'
3 /40X, 'NO FACTOR ANALYSIS WILL BE ATTEMPTED, ONE NEED AT LEAST'/
4 40X, ' 4 VARIABLES.')
IF (NAT.LT.4) CALL END (2)
IF (NAT.LT.4) CALL END (3)
CALL XDIMCH(NRMS1,NAT,'XFAWMV',0)
CALL XDIMCH(NRMS2,NAT,'XFAWMV',0)
CALL XDIMCH(NRMS3,NAT,'XFAWMV',0)
CALL XDIMCH(NRMS4,NAT,'XFAWMV',0)
CALL XEIGEN(NRMS1,NRMS3,NRMS4,S1,S3,S4,NAT,NAT)
IF (S4(NAT,NAT).LT.0.000001D0) GOTO 79
IZERO = 0
CALL XRECIP (NRMS4,NRMS2,S4,S2,NAT,NAT,IZERO)
GOTO 84
79 CALL TITLE ('MATRIX IS SINGULAR, HIGHEST COVARIANCE IS USED AS COM
2MUNALITY ESTIMATE',80)
GOTO 85
ENTRY XFAHIR (NRMS1,NCMS1,NRMS2,NRMS3,NCMS3,NRMS4,S1,S2,
2 S3,S4,NAT,NCF)
CALL QINCR ('XFAHIR')
IF (NAT.LT.4) CALL END (1)
IF (NAT.LT.4) WRITE (LUWN, 77) NAT
IF (NAT.LT.4) CALL END (2)
IF (NAT.LT.4) CALL END (3)
CALL XDIMCH(NRMS1,NAT,'XFAHIR',0)
CALL XDIMCH(NRMS2,NAT,'XFAHIR',0)
CALL XDIMCH(NRMS3,NAT,'XFAHIR',0)
CALL XDIMCH(NRMS4,NAT,'XFAHIR',0)
85 N1 = NAT - 1
S2(1,1) = DABS(S1(1,2))
DO 80 I = 2, NAT
80 S2(I,I) = DABS(S1(1,I))
DO 81 I = 1, N1
K = I + 1
DO 81 J = K, NAT
IF (DABS(S1(I,J)).GT.S2(I,I)) S2(I,I) = DABS(S1(I,J))
S2(I,J) = S1(I,J)
81 S2(J,I) = S1(I,J)
DO 82 I = 2, N1
82 IF (DABS(S1(NAT,I)).GT.S2(NAT,NAT)) S2(NAT,NAT) = DABS(S1(NAT,I))
CALL TITLE ('FACTOR ANALYSIS OF COVARIANCE MATRIX WITH HIGHEST COV
2ARIANCE IN DIAGONAL',72)
GOTO 83
84 CALL TITLE('FACTOR ANALYSIS OF COVARIANCE MATRIX WITH MULTIPLE VAR
1IANCES IN DIAGONAL',72)
C COMPUTE MV IN DIAGONAL OF C
DO 71 I = 1, NAT
DO 71 J = 1, NAT
S4(I,J) = 0.0D0
DO 71 K = 1, NAT
71 S4(I,J) = S4(I,J) + S3(I,K)*S2(K,K)*S3(J,K)
CALL XCOPY(NRMS1,NRMS2,S1,S2,NAT,NAT)
CALL TITLE('MULTIPLE VARIANCES',18)
DO 10 I = 1,NAT
T1 = 1.0D0/S4(I,I)
S2(I,I) = S2(I,I) - T1
10 WRITE(LUWN,900)I,S2(I,I)
900 FORMAT('0',I5,F15.5)
C EIGENSOLUTION OF CWMV
83 CALL XEIGEN(NRMS2,NRMS3,NRMS4,S2,S3,S4,NAT,NCF)
CALL TITLE('EIGENVALUES OF COVARIANCE MATRIX MINUS COMMUNALITY ES
2TIMATE',67)
CALL TITLE(' EIGENVALUES DIFFERENCES',36)
WRITE(LUWN,901)
DO 15 I = 1,NAT
IP = I + 1
WRITE(LUWN,901)I,S4(I,I)
901 FORMAT(' ',I5,F15.5)
IF(IP.GT.NAT) GO TO 15
T1 = S4(I,I) - S4(IP,IP)
WRITE(LUWN,902) T1
902 FORMAT(' ',20X,F15.5)
15 CONTINUE
C TEST POSITIVENESS OF EIGENVALUES AND REDUCE NCF WHEN NECESSARY
NCFT = 0
DO 20 I = 1,NCF
IF(S4(I,I).LE.0.0D0) GO TO 25
NCFT = NCFT + 1
20 CONTINUE
WRITE(LUWN,903)NCF
903 FORMAT('-','FACTOR SOLUTION FOR',I4,2X,'FACTORS')
GO TO 30
25 NCF = NCFT
WRITE(LUWN,904)NCF
904 FORMAT('-','NCF REDUCED TO',I4,2X,'FACTORS')
C FACTOR MATRIX
30 DO 45 J = 1,NCF
T1 = 0.0D0
T2 = DSQRT(S4(J,J))
DO 35 I = 1,NAT
S3(I,J) = S3(I,J)*T2
35 T1 = T1 + S3(I,J)
IF(T1.GE.0.0D0) GO TO 45
DO 40 I = 1,NAT
40 S3(I,J) = -S3(I,J)
45 CONTINUE
CALL TITLE('FACTOR MATRIX',13)
CALL XPRINT(NRMS3,S3,NAT,NCF,'F')
C COMMUNALITIES
CALL TITLE('COMMUNALITIES FOR GIVEN NUMBER OF FACTORS',41)
DO 55 I = 1,NAT
S4(I,1) = S3(I,1)**2
IF(NCF.LE.1) GO TO 55
DO 50 J = 2,NCF
JM = J - 1
50 S4(I,J) = S4(I,JM) + S3(I,J)**2
55 CONTINUE
CALL XPRINT(NRMS4,S4,NAT,NCF,'F')
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE ITPFA(C,A,HV,S1,S2,NAT,NF)
C ITERATIVE PRINCIPAL AXES FACTOR ANALYSIS
C NOTE: THIS PROGRAM DOES NOT HAVE A PROTECTION AGAINST HEYWOOD CASES
C PROGRAM WRITTEN BY LEDYARD R TUCKER, APRIL 1978
C RUNS ON FORMAL
C C IS COVARIANCE MATRIX; INPUT/OUTPUT
C A IS RESULTING FACTOR MATRIX;OUTPUT
C HV IS A FORMAL MATRIX WITH COLUMNS USED AS VECTORS:
C 1 FOR COMMUNALITIES OF PREVIOUS TRIAL,
C 2 FOR COMMUNALITIES COMPUTED FROM FACTOR MATRIX,
C 3 FOR DIFFERENCES BETWEEN INPUT COMMUNALITIES AND COMPUTED
C COMMUNALITIES FOR PREVIOUS TRIAL,
C 4 FOR DIFFERENCES FOR PRESENT TRIAL.
C S1 CONTAINS THE EIGENVALUES FOR THE PRESENT TRIAL IN THE DIAGONAL; OUTPUT
C S2 CONTAINS THE COVARIANCE MATRIX WITH THE OBTAINED COMMUNALITIES; OUTPUT
C NAT IS THE NUMBER OF ATTRIBUTES; INPUT/OUTPUT
C NF IS THE NUMBER OF FACTORS; INPUT/OUTPUT
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION C(I1,1), A(I1,1), HV(I1,1), S1(I1,1), S2(I1,1)
COMMON/B/I1
COMMON/ZINCR/IZINCR
CALL QINCR ('ITPFA')
NRMC = I1
NRMA = I1
NRMHV= I1
NRMS1= I1
NCMS1= I1
NRMS2= I1
NCMS2= I1
CALL XITPFA(NRMC,NRMA,NRMHV,NRMS1,NCMS1,NRMS2,NCMS2,
2 C,A,HV,S1,S2,NAT,NF)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XITPFA(NRMC,NRMA,NRMHV,NRMS1,NCMS1,NRMS2,NCMS2,
2 C,A,HV,S1,S2,NAT,NF)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION C(NRMC,1),A(NRMA,1),HV(NRMHV,1),S1(NRMS1,NCMS1)
DIMENSION S2(NRMS2,NCMS2)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
CALL QINCR ('XITPFA')
CALL XDIMCH(NRMC,NAT,'ITPFA',0)
CALL XDIMCH(NRMA,NAT,'ITPFA',0)
CALL XDIMCH(NRMHV,NAT,'ITPFA',0)
CALL XDIMCH(NRMS1,NAT,'ITPFA',0)
CALL XDIMCH(NRMS2,NAT,'ITPFA',0)
CALL TITLE('ITERATIVE PRINCIPAL FACTOR ANALYSIS',35)
NFS = NF
T1 = NAT
T1 = T1/2.0D0
ITEM = T1
IF(NF.LE.ITEM) GO TO 5
NF = ITEM
5 CONTINUE
EPS = 1.0D-5
EPSE = 1.0D-20
ST = 1.0D0
NCYLT = 20
NCY = 0
NFP1 = NF + 1
C INITIAL HSQ = MV
CALL XCOPY(NRMC,NRMS2,C,S2,NAT,NAT)
CALL XINVRS(NRMS2,NCMS2,NRMS1,NCMS1,S2,S1,NAT)
DO 10 I = 1,NAT
HV(I,1) = S2(I,I)
10 HV(I,3) = -1.0D0/S1(I,I)
GO TO 105
C START OF FULL CYCLE
100 NCY = NCY + 1
C COMPUTE NEW DIAGONALS
105 DO 110 I = 1,NAT
110 S2(I,I) = HV(I,1) + ST*HV(I,3)
C OBTAIN EIGENSOLUTION, FACTOR MATRIX, OUTPUT COMMUNALITIES
CALL XEIGEN(NRMS2,NRMA,NRMS1,S2,A,S1,NAT,NF)
C FOR NCY = 0 CHECK NF
IF(NCY.GT.0) GO TO 114
ITEM = 0
DO 112 I = 1,NF
IF(S1(I,I).GT.EPSE) ITEM = ITEM + 1
112 CONTINUE
NF = ITEM
114 DO 115 J = 1,NF
T1 = S1(J,J)
IF(T1.LT.EPSE) T1 = EPSE
T1 = SQRT(T1)
DO 115 I = 1,NAT
115 A(I,J) = A(I,J)*T1
DO 120 I = 1,NAT
HV(I,2) = 0.0D0
DO 120 J = 1,NF
120 HV(I,2) = HV(I,2) + A(I,J)**2
C OBTAIN NEW D AND CRITERION
SSD = 0.0D0
DO 125 I = 1,NAT
HV(I,4) = HV(I,2) - S2(I,I)
125 SSD = SSD + HV(I,4)**2
CRT = 0.0D0
DO 130 I = NFP1,NAT
130 CRT = CRT + S1(I,I)**2
CRT = CRT - SSD
C TEST NEW D'S AGAINST EPS
DO 135 I = 1,NAT
IF(ABS(HV(I,4)).GT.EPS) GO TO 140
135 CONTINUE
GO TO 165
140 IF(NCY.LE.1) GO TO 155
C TEST FOR REDUCED CRITERION; OTHERWISE REDUCE ST AND REPEAT
IF(CRT.LT.CRTP) GO TO 145
ST = ST/2.0D0
IF(ST.LT.EPS) GO TO 165
GO TO 105
145 CONTINUE
C COMPUTE COS DP WITH D AND ADJUST ST
COS = 0.0D0
DO 150 I = 1,NAT
150 COS = COS + HV(I,3)*HV(I,4)
T1 = SQRT(SSDP*SSD)
COS = COS/T1
ST = ST*(3.0D0 + COS)/(3.0D0 - COS)
C MOVE DIAGONAL CWHSQ TO HSQP, D TO DP, SSD TO SSDP AND CRT TO CRTP
155 DO 160 I = 1,NAT
HV(I,1) = S2(I,I)
160 HV(I,3) = HV(I,4)
SSDP = SSD
CRTP = CRT
C TEST NUMBER OF CYCLES, REPEAT WHEN NCY.LT.NCYLT
IF(NCY.LT.NCYLT) GO TO 100
C PRINT RESULTS
WRITE(LUWN,910)NCY
910 FORMAT('0','COMPUTATIONS STOPPED AFTER',I4,3X,'CYCLES')
GO TO 170
165 WRITE(LUWN,911)NCY
911 FORMAT('0','CONVERGENCE OBTAINED IN',I4,3X,'CYCLES')
170 IF(NF.LT.NFS) WRITE(LUWN,909)NF
909 FORMAT('0','NUMBER OF FACTORS REDUCED TO',I4)
CALL TITLE('OBTAINED COMMUNALITIES',22)
DO 175 I = 1,NAT
175 WRITE(LUWN,912)I,HV(I,2)
912 FORMAT('0',I6,F14.5)
CALL TITLE('EIGENVALUES',11)
CALL PRDIF(S1,NAT,NRMS1)
C REVERSE FACTORS WHEN NECCESSARY
DO 185 J = 1,NF
T1 = 0.0D0
DO 180 I = 1,NAT
180 T1 = T1 + A(I,J)
IF(T1.GE.0.0D0) GO TO 185
DO 182 I = 1,NAT
182 A(I,J) = -A(I,J)
185 CONTINUE
CALL TITLE('FACTOR MATRIX',13)
CALL XPRINT(NRMA,A,NAT,NF,'F')
ITEM = NAT*(NAT-1)
IF(CRT.LT.EPSE) CRT = EPSE
T1 = CRT/ITEM
T1 = SQRT(T1)
WRITE(LUWN,913)T1
913 FORMAT('0','ROOT MEAN SQUARE RESIDUAL',F14.5)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE ABS(X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO), Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('ABS')
NRMX=IO
NRMY=IO
CALL XABS(NRMX,NRMY,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE ADD(X,Y,Z,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO), Y(IO,IO), Z(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('ADD')
NRMX=IO
NRMY=IO
NRMZ=IO
CALL XADD(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE CHLSK(X,Y,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO), Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('CHLSK')
NRMX=IO
NRMY=IO
CALL XCHLSK(NRMX,NRMY,X,Y,N)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE CHSYM(X,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('CHSYM')
NRMX=IO
CALL XCHSYM(NRMX,X,N)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE COLSM(X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO), Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('COLSM')
NRMX=IO
NRMY=IO
CALL XCOLSM(NRMX,NRMY,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE COPY(X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO), Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('COPY')
NRMX=IO
NRMY=IO
CALL XCOPY(NRMX,NRMY,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE DET(X,D,N,Y)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO), Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('DET')
NRMX=IO
NCMX=IO
NRMY=IO
NCMY=IO
CALL XDET(NRMX,NCMX,NRMY,NCMY,X,D,N,Y)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE DIAG(X,Y,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO), Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('DIAG')
NRMX=IO
NRMY=IO
CALL XDIAG(NRMX,NRMY,X,Y,N)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE EIGEN(R,V,D,N,NVE)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION R(IO,IO),V(IO,IO),D(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('EIGEN')
NRMR=IO
NRMV=IO
NRMD=IO
CALL XEIGEN(NRMR,NRMV,NRMD,R,V,D,N,NVE)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE ELEDI(X,Y,Z,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO),Z(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('ELEDI')
NRMX=IO
NRMY=IO
NRMZ=IO
CALL XELEDI(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE ELEML(X,Y,Z,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO),Z(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('ELEML')
NRMX=IO
NRMY=IO
NRMZ=IO
CALL XELEML(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE FNFRT(A,NMP,TP,PHI,PJT,B,Q,NAT,NCF,ICN,ICT,ICP,
1ICB)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(IO,IO),TP(IO,IO),PHI(IO,IO),PJT(IO,IO),B(IO,IO)
DIMENSION Q(IO,IO)
DOUBLE PRECISION NMP(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('FNFRT')
NRMA=IO
NRMNMP=IO
NRMTP=IO
NRMPHI=IO
NRMPJT=IO
NRMB=IO
NRMQ=IO
CALL XFNFRT(NRMA,NRMNMP,NRMTP,NRMPHI,NCMPHI,NRMPJT,NRMB,NRMQ,
2 NCMQ,A,NMP,TP,PHI,PJT,B,Q,NAT,NCF,ICN,ICT,ICP,ICB)
IZINCR = IZINCR - 1
RETURN
END
C SUBROUTINE XFNFRT(NRMA,NRMNMP,NRMTP,NRMPHI,NCMPHI,NRMPJT,
C 2 NRMB,NRMQ,NCMQ,A,NMP,TP,PHI,PJT,B,Q,NAT,NCF,ICN,ICT,ICP,ICB)
C SUBROUTINE FNFRT
C PERFORMS FINAL COMPUTATIONS FOR FACTOR ANALYSIS ROTATION OF AXES
C RESULTS ARE PRINTED AND PASSED TO MAIN PROGRAM
C
C PROGRAM WRITTEN BY LEDYARD R TUCKER
C MARCH 1976
C
C USES FORMAL
C
C NUMBER OF COMMON FACTORS MUST BE 2 OR GREATER
C AND AT MOST 1 LESS THAN THE MATRIX DIMENSIONS IN FORMAL
C INPUT-OUTPUT MATRICES
C A = COORDINATE AXES FACTOR MATRIX
C NMP = A DOUBLE PRECISION MATRIX FOR NORMALS, N'
C INPUT FORM CONTROLLED BY ICN
C ICN = 0 FOR N' INPUT (NORMALS ARE COLUMNS)
C = 1 FOR N INPUT (NORMALS ARE ROWS)
C = 2 FOR N' TO BE COMPUTED FROM T'
C NORMALS ARE ASSUMED TO BE UNIT VECTORS WHEN INPUT
C OUTPUT N' (NORMALS ARE COLUMNS)
C TP = MATRIX FOR PRIMARY VECTORS
C INPUT FORM CONTROLLED BY ICT
C ICT = 0 FOR T' INPUT (PRIMARY VECTORS ARE COLUMNS)
C = 1 FOR T INPUT (PRIMARY VECTORS ARE ROWS)
C = 2 FOR T' TO BE COMPUTED FROM N'
C PRIMARY VECTORS ARE ASSUMED TO BE UNIT VECTORS WHEN INPUT
C OUTPUT T' (PRIMARY VECTORS ARE COLUMNS)
C NOTE EITHER NMP OR TP OR BOTH MUST BE INPUT
C PHI = FACTOR INTERCORRELATION MATRIX ON OUTPUT = TP'*TP
C ROW NCF + 1 CONTAINS DIAGONAL ENTRIES OF MATRIX D
C D = (DIAG(PHI-INVERSE))**-.5
C PJT = MATRIX FOR PROJECTIONS ON NORMALS = A*NMP
C INPUT FORM CONTROLLED BY ICP
C ICP = 0 FOR PJT INPUT
C = 1 FOR PJT TO BE CALCULATED
C PJT IS OUTPUT
C B = FACTOR WEIGHT MATRIX = PJT*D-INVERSE
C INPUT FORM CONTROLLED BY ICB
C ICB = 0 FOR B INPUT
C = 1 FOR B TO BE COMPUTED
C B IS OUTPUT
C Q = MATRIX OF COVARIANCES OF ATTRIBUTES WITH FACTORS= B*PHI
C Q IS OUTPUT
C INPUT SCALARS
C NAT = NUMBER OF ATTRIBUTES
C NCF = NUMBER OF COMMON FACTORS
C CONTROL FLAGS DESCRIBED ABOVE ICN, ICT, ICP, ICB
C XDIMCH & QINCR ARE FORMAL UTILITY SUBROUTINES
C COMMON /A/ CONTAINS THE PROGRAM COUNTER VARIABLE
SUBROUTINE XFNFRT(NRMA,NRMNMP,NRMTP,NRMPHI,NCMPHI,NRMPJT,
2 NRMB,NRMQ,NCMQ,A,NMP,TP,PHI,PJT,B,Q,NAT,NCF,ICN,ICT,ICP,ICB)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION A(NRMA,1),TP(NRMTP,1),PHI(NRMPHI,1),PJT(NRMPJT,1)
DIMENSION Q(NRMQ,1),B(NRMB,1)
DOUBLE PRECISION NMP(NRMNMP,1)
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
901 FORMAT('-',5X,'NEITHER N'' (OR N) NOR T'' (OR T) WAS INPUT')
902 FORMAT('0',5X,'COMPUTATIONS ARE NOT POSSIBLE BY THIS SUBROUTINE')
CALL QINCR ('XFNFRT')
IF(NCF.LT.2)GOTO101
NCFP=NCF+1
NCFM=NCF-1
CALL XDIMCH(NRMA,NAT,'FNFRT',0)
CALL XDIMCH(NRMNMP,NAT,'FNFRT',0)
CALL XDIMCH(NRMTPS,NAT,'FNFRT',0)
CALL XDIMCH(NRMPHI,NAT,'FNFRT',0)
CALL XDIMCH(NRMPJT,NAT,'FNFRT',0)
CALL XDIMCH(NRMB,NAT,'FNFRT',0)
CALL XDIMCH(NRMQ,NAT,'FNFRT',0)
CALL XDIMCH(NRMA,NCFP,'FNFRT',0)
CALL XDIMCH(NRMNMP,NCFP,'FNFRT',0)
CALL XDIMCH(NRMTPS,NCFP,'FNFRT',0)
CALL XDIMCH(NRMPHI,NCFP,'FNFRT',0)
CALL XDIMCH(NRMPJT,NCFP,'FNFRT',0)
CALL XDIMCH(NRMB,NCFP,'FNFRT',0)
CALL XDIMCH(NRMQ,NCFP,'FNFRT',0)
C BRANCH FOR N' (OR N) NOT INPUT
IF(ICN.GT.1)GOTO100
C IF N INPUT, TRANSPOSE TO N'
IF(ICN.EQ.0)GOTO15
DO 10I=1,NCFM
IP=I+1
DO 10J=IP,NCF
TEM1=NMP(I,J)
NMP(I,J)=NMP(J,I)
10 NMP(J,I)=TEM1
C COMPUTE NN' STORED IN Q
15 CALL XXTX(NRMNMP,NRMQ,NMP,Q,NCF,NCF)
C IF T' (OR T) NOT INPUT, COMPUTE T'
IF(ICT.LT.2)GOTO35
C COMPUTE (NN') INVERSE, STORE IN PHI
CALL XINVRS(NRMQ,NCMQ,NRMPHI,NCMPHI,Q,PHI,NCF)
C COMPUTE DIAGONAL ENTRIES OF D, STORE IN ROW NCFP OF PHI
DO 20I=1,NCF
20 PHI(NCFP,I)=1.0D0/DSQRT(PHI(I,I))
C COMPUTE T'
CALL XMATML(NRMNMP,NRMPHI,NRMTP,NMP,PHI,TP,NCF,NCF,NCF)
DO 25I=1,NCF
DO 25J=1,NCF
25 TP(I,J)=TP(I,J)*PHI(NCFP,J)
C COMPUTE PHI
DO 30I=1,NCF
DO 30J=I,NCF
PHI(I,J)=PHI(NCFP,I)*PHI(I,J)*PHI(NCFP,J)
30 PHI(J,I)=PHI(I,J)
C OUT TO PRINT TRANSFORMATION MATRICES
GOTO 200
C T' (OR T) INPUT
C IF T INPUT, TRNSPOSE TO T'
35 IF(ICT.EQ.0)GOTO45
DO 40I=1,NCFM
IP=I+1
DO 40J=IP,NCF
TEM1=TP(I,J)
TP(I,J)=TP(J,I)
40 TP(J,I)=TEM1
C COMPUTE PHI AND D
45 CALL XXTX(NRMTP,NRMPHI,TP,PHI,NCF,NCF)
DO 50I=1,NCF
PHI(NCFP,I)=0.D0
DO 50J=1,NCF
50 PHI(NCFP,I)=PHI(NCFP,I)+NMP(J,I)*TP(J,I)
C OUT TO PRINT TRANSFORMATION MATRICES
GOTO 200
C WHEN N' (OR N) WAS NOT INPUT
C TEST THAT T' (OR T) WAS INPUT, OTHERWISE NOTE AND RETURN
100 IF(ICT.LT.2)GOTO105
CALL END (1)
WRITE(LUWN,901)
WRITE(LUWN,902)
CALL END (2)
IZINCR = IZINCR - 1
RETURN
C IF T INPUT, TRNSPOSE TO T'
105 IF(ICT.EQ.0)GOTO115
DO 110I=1,NCFM
IP=I+1
DO 110J=IP,NCF
TEM1=TP(I,J)
TP(I,J)=TP(J,I)
110 TP(J,I)=TEM1
C COMPUTE PHI
115 CALL XXTX(NRMTP,NRMPHI,TP,PHI,NCF,NCF)
C COMPUTE PHI INVERSE, D, NN', N'
CALL XINVRS(NRMPHI,NCMPHI,NRMQ,NCMQ,PHI,Q,NCF)
DO 120I=1,NCF
120 PHI(NCFP,I)=1.0D0/DSQRT(Q(I,I))
CALL XMATML(NRMTP,NRMQ,NRMNMP,TP,Q,NMP,NCF,NCF,NCF)
DO 125I=1,NCF
DO 125J=1,NCF
125 NMP(I,J)=PHI(NCFP,J)*NMP(I,J)
DO 130I=1,NCF
DO 130J=I,NCF
Q(I,J)=PHI(NCFP,I)*Q(I,J)*PHI(NCFP,J)
130 Q(J,I)=Q(I,J)
C PRINT TRANSFORMATION MATRICES
200 CALL TITLE('MATRIX N'' OF NORMALS TO HYPERPLANES',35)
CALL TITLE('R0WS COORDINATE AXES, COLUMNS NORMALS',41)
CALL XPRINT(NRMNMP,NMP,NCF,NCF,'F')
CALL TITLE('MATRIX NN'' OF COSINES OF ANGLES BETWEEN NORMALS',47)
CALL XPRINT(NRMQ,Q,NCF,NCF,'F')
CALL TITLE('MATRIX T'' OF PRIMARY VECTORS',28)
CALL TITLE('ROWS COORDINATE AXES, COLUMNS PRIMARY VECTORS',
+49)
CALL XPRINT(NRMTP,TP,NCF,NCF,'F')
CALL TITLE('MATRIX PHI OF FACTOR INTERCORRELATIONS',38)
CALL XPRINT(NRMPHI,PHI,NCF,NCF,'F')
DO 205I=1,NCF
205 Q(1,I)=PHI(NCFP,I)
CALL TITLE('DIAGONAL ENTRIES OF MATRIX D',28)
CALL TITLE(
+'COSINES OF ANGLES BETWEEN CORRESPONDING NORMALS AND PRIMARY VECTO
+RS',67)
CALL XPRINT(NRMQ,Q,1,NCF,'F')
C COMPUTE AND PRINT ATTRIBUTE - FACTOR MATRICES
C PROJECTIONS ON NORMALS
IF(ICP.EQ.0)GOTO225
CALL XMATML(NRMA,NRMNMP,NRMPJT,A,NMP,PJT,NAT,NCF,NCF)
225 CALL TITLE('ATTRIBUTE - PART FACTOR RELATION MATRIX P',42)
CALL TITLE('(PROJECTIONS ON NORMALS)',24)
CALL TITLE('(STRUCTURE LOADINGS ON REFERENCE VECTORS)',41)
CALL XPRINT(NRMPJT,PJT,NAT,NCF,'F')
C FACTOR WEIGHTS
IF(ICB.EQ.0)GOTO235
DO 230I=1,NAT
DO 230J=1,NCF
230 B(I,J)=PJT(I,J)/PHI(NCFP,J)
235 CALL TITLE('FACTOR WEIGHT MATRIX B',23)
CALL TITLE('(PATTERN LOADINGS ON PRIMARY VECTORS)',37)
CALL XPRINT(NRMB,B,NAT,NCF,'F')
C COVARIANCES OF ATTRIBUTES WITH FACTORS
CALL XMATML(NRMB,NRMPHI,NRMQ,B,PHI,Q,NAT,NCF,NCF)
CALL TITLE('ATTRIBUTE - FACTOR RELATION MATRIX Q',37)
CALL TITLE('(COVARIANCES OF ATTRIBUTES WITH FACTORS)',40)
CALL TITLE('(STRUCTURE LOADINGS ON PRIMARY VECTORS)',39)
CALL XPRINT(NRMQ,Q,NAT,NCF,'F')
IZINCR = IZINCR - 1
RETURN
101 CALL END (1)
WRITE(LUWN,102)NCF,NAT
102 FORMAT('0* * * * * FATAL ERROR. THE NUMBER OF COMMON FACTORS SPEC
1IFIED (',I5,') IS SUPPOSED TO BE .GE.2 AND .LT.',I4)
CALL END (2)
CALL END (3)
END
SUBROUTINE GENR8(C,X,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('GENR8')
CALL XGENR8(IO,C,X,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE GINVR(X,Y,NR,NC,A,D,C)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION X(IO,IO),Y(IO,IO),A(IO,IO),D(IO,IO),C(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('GINVR')
NRMX=IO
NRMY=IO
NRMA=IO
NRMD=IO
NRMC=IO
CALL XGINVR(NRMX,NRMY,NRMA,NRMD,NRMC,X,Y,NR,NC,A,D,C)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE HORIZ(X,Y,Z,NR,NCX,NCY)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO),Z(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('HORIZ')
NRMX=IO
NRMY=IO
NRMZ=IO
CALL XHORIZ(NRMX,NRMY,NRMZ,X,Y,Z,NR,NCX,NCY)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE IDENT(X,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('IDENT')
CALL XIDENT(IO,X,N)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE INPUT(X)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('INPUT')
CALL XINPUT(IO,X)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE TINPT (XINPUT, X)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('TINPT')
CALL XTINPT(IO,XINPUT,X)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE INVRS(X,Y,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('INVRS')
NRMX=IO
NCMX=IO
NRMY=IO
NCMY=IO
CALL XINVRS(NRMX,NCMX,NRMY,NCMY,X,Y,N)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE KRONK(X,Y,Z,NRX,NCX,NRY,NCY)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO),Z(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('KRONK')
NRMX=IO
NRMY=IO
NRMZ=IO
CALL XKRONK(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NCX,NRY,NCY)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE LSQHY(A,S,N,W,D,C,NR,NC,IFL)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION A(IO,IO),N(IO,IO),S(IO,IO),
2 W(IO,IO),D(IO,IO),C(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('LSQHY')
NRMA=IO
NRMS=IO
NRMN=IO
NRMW=IO
NRMD=IO
NRMC=IO
CALL XLSQHY(NRMA,NRMS,NRMN,NRMW,NRMD,NRMC,A,S,N,W,D,C,
2 NR,NC,IFL)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE MATML(X,Y,Z,NRX,NXY,NCY)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO),Z(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('MATML')
NRMX=IO
NRMY=IO
NRMZ=IO
CALL XMATML(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NXY,NCY)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE PAGE
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
CALL QINCR ('PAGE')
WRITE (LUWN,1)
1 FORMAT ('1')
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE PARTT(X,Y,IBR,IER,IBC,IEC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('PARTT')
NRMX=IO
NRMY=IO
CALL XPARTT(NRMX,NRMY,X,Y,IBR,IER,IBC,IEC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE PRDDI(X,N,PD)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('PRDDI')
CALL XPRDDI(IO,X,N,PD)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE PRINT(X,NR,NC,IFLAG)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('PRINT')
CALL XPRINT(IO,X,NR,NC,IFLAG)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE PUNCH(X,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('PUNCH')
CALL XPUNCH(IO,X,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE RANK(X,Y,NR,NC,RNK)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('RANK')
NRMX=IO
NCMX=IO
NRMY=IO
CALL XRANK(NRMX,NCMX,NRMY,X,Y,NR,NC,RNK)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE RECIP(X,Y,NR,NC,IZERO)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('RECIP')
NRMX=IO
NRMY=IO
CALL XRECIP(NRMX,NRMY,X,Y,NR,NC,IZERO)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE RO2DI(X,Y,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('RO2DI')
NRMX=IO
NRMY=IO
CALL XRO2DI(NRMX,NRMY,X,Y,N)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE RO2MT(X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('RO2MT')
NRMX=IO
NRMY=IO
CALL XRO2MT(NRMX,NRMY,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE SCAML(C,X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('SCAML')
NRMX=IO
NRMY=IO
CALL XSCAML(NRMX,NRMY,C,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE SIMLN(X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('SIMLN')
NRMX=IO
NRMY=IO
CALL XSIMLN(NRMX,NRMY,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XSIMLN(NRMX,NRMY,X,Y,NR,NC)
C TITLE: SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS
C PROGRAMMER: CARL FINKBEINER
C DATE: AUGUST, 1975
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: FIND THE SOLUTION SET OF SEVERAL LINEAR EQUATIONS IN SEVERAL
C UNKNOWNS.
C SYMBOLS:
C X = INPUT MATRIX OF THE FORM (A:B) WHERE A IS NR X NR, C IS NR X N1
C Y = RESULTING SOLUTION FOR Y IN THE EQUATION A*Y = C
C NR = # OF ROWS OF X; # OF ROWS & COLS OF A
C NC = # OF COLS OF X
C QINCR = PROGRAM COUNTER ROUTINE
C XDIMCH = DIMENSIONALITY CHECK ROUTINE
C QWRGDM = ERROR MESSAGE ROUTINE
C IER = ERROR FLAG USED BY QGELGZ
C N = BEGINNING COL OF C PORTION OF X
C N1 = # OF COLS OF C PORTION OF X
C QGELGZ = NUMBER CRUNCHER ROUTINE
C REORD = A REORDERING ROUTINE FOR THE RESULT MATRIX
C I,J,K = SCRATCH VARIABLES
C PROCEDURE: THE ORDER OF X IS CHECKED FOR APPROPRIATENESS, NAMELY THE
C # OF ROWS CANNOT BE 1 OR GREATER THAN OR EQUAL TO THE # OF COLS.
C THE # OF COLS IN THE C PORTION OF X IS COMPUTED, THE SOLUTION FOR Y
C IN THE EQUATION A*Y = C IS FOUND AND STORED IN THE UPPER LEFT CORNER
C OF Y AS A MATRIX THE SAME ORDER AS C. A AND C ARE CONCATENATED (A:
C AND INPUT TO THIS ROUTINE AS X. X AND Y MAY BE THE SAME MATRIX IN
C THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XSIMLN')
CALL XDIMCH(NRMX,NR,'SIMLN',0)
CALL XDIMCH(NRMY,NC,'SIMLN',0)
IF(NR.EQ.1)CALL QWRGDM(NR,-9999)
IF(NR.GE.NC)CALL QWRGDM(NR,NC)
C FIND C PORTION OF X AND FIND SOLUTION
IER=1
N=NR+1
N1=NC-NR
CALL QGELGZ(NRMX,NRMX,X(1,N),X,NR,N1,Y(1,N),Y,1.D-7,IER)
IF(IER.NE.0)CALL QWRGDM(-9999,-9999)
C GET Y INTO EXPECTED FORM
CALL XREORD(NRMY,Y(1,N),NR,N1)
J=NR
DO 1I=1,N1
J=J+1
DO 1K=1,NR
1 Y(K,I)=Y(K,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE SQRT(X,Y,NR,NC,EPS)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('SQRT')
NRMX=IO
NRMY=IO
CALL XSQRT(NRMX,NRMY,X,Y,NR,NC,EPS)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE SUB(X,Y,Z,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO),Z(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('SUB')
NRMX=IO
NRMY=IO
NRMZ=IO
CALL XSUB(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE TRACE(X,N,T)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('TRACE')
CALL XTRACE(IO,X,N,T)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE TRNSP(X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('TRNSP')
NRMX=IO
NRMY=IO
CALL XTRNSP(NRMX,NRMY,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE VRTCL(X,Y,Z,NRX,NRY,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO),Z(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('VRTCL')
NRMX=IO
NRMY=IO
NRMZ=IO
CALL XVRTCL(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NRY,NC)
IZINCR = IZINCR - 1
RETURN
END
C SUBROUTINE XRANK(NRMX,NCMX,NRMY,X,Y,NR,NC,RNK)
C TITLE:RANK
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO FIND THE RANK OF A REAL MATRIX
C SYMBOLS:
C X = INPUT MATRIX
C Y = SCRATCH ARRAY
C NR = # OF ROWS OF X
C NC = # OF COLS OF X
C RNK = SMALLEST NON-ZERO EIGENVALUE TO BE COUNTED: INPUT BY USER =
C RANK
C OF X ON OUTPUT.
C A = LABELLED COMMON AREAS CONTAINING IPGM.
C IPGM = PROGRAM COUNTER
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C XXTX = CROSS PRODUCTS SUBROUTINE
C NX = SCRATCH VARIABLE
C QTRED1 = EISPACK TRIDIAGONALIZATION SUBROUTINE
C QIMTQL = EISPACK EIGENVALUES OF A TRIDIAGONAL MATRIX SUBROUTINE
C IER = ERROR FLAG FOR QIMTQL
C PROCEDURE: THE PROGRAM CALLS XXTX TO FIND THE CROSS-PRODUCTS OF X (NOTE
C IPGM MUST BE SET BACK TO COMPENSATE HERE). THE RANK OF X'X = RANK
C OF NON-ZERO EIGENVALUES OF X'X. HENCE THE EISPACK SUBROUTINES ARE
C ISSUED FOR FINDING THE EIGENVALUES OF X'X. X AND Y MAY BE THE SAME
C ADDRESS IF THE USER DOESN'T MIND THAT X WILL BE DESTROYED.
C FUNCTIONAL 0 IS USED BY THE USER THRU RNK ON INPUT. A NUMBER IN THE
C RANGE 1.D-5 TO 1.D-10 USUALLY ADEQUATE. (MATHEMATICAL NOTE: THE
C THING THAT MAKES THIS ALSO WORK IS THE FACT THAT X'X (OR XX') IS
C ALWAYS POSITIVE DEFINITE WHEN COMPOSED OF REAL NUMBERS.
SUBROUTINE XRANK(NRMX,NCMX,NRMY,X,Y,NR,NC,RNK)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),RNK,EPS
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XRANK')
NCP1=NC+1
NCP2=NC+2
CALL XDIMCH(NRMX,NR,'RANK',0)
CALL XDIMCH(NRMY,NC,'RANK',0)
CALL XDIMCH(NCMX,NCP2,'RANK',0)
IF(NCP1.GE.NCMX)GOTO 5
IF(NCP1.GE.NRMX)GOTO 5
IF(NCP1.GE.NRMY)GOTO 5
C FIND CROSS-PRODUCTS OF X AND STORE IN Y
CALL XXTX(NRMX,NRMY,X,Y,NR,NC)
NX=NCP1
NX1=NCP2
C USE EISPACK SUBROUTINES FOR TRIDIAGONALIZING THE CROSS-PRODUCTS MATRIX
C THEN FINDING THE EIGENVALUES.
CALL QTRED1(NX,NC,Y,X(1,NX),Y(1,NX),Y(1,NX1))
CALL QIMTQL(NC,X(1,NX),Y(1,NX),IER)
IF(IER.NE.0)GOTO3
EPS=RNK
RANK=0.0D0
C COPY THE EIGENVALUES OF THE CROSS-PRODUCTS MATRIX INTO THE 1ST COL OF
C AND CHECK FOR # OF EIGENVALUES .GT. EPS
J=NC+1
DO 1 I=1,NC
J=J-1
IF(X(J,NX).GT.EPS)RNK=RNK+1.D0
1 CONTINUE
C ROUND RNK TO AN EXACT INTEGER
I=RNK+.499999D0
RNK=I
IZINCR = IZINCR - 1
RETURN
3 CALL END (1)
WRITE(LUWN,4)
4 FORMAT(13X,' * * * THE ALGORITHM FOR FINDING THE RANK OF A MATRIX
1WILL NOT WORK FOR THE MATRIX PASSED'/13X' THE SUBROUTINE FAILED W
2HILE CALCULATING EIGENVALUE NUMBER',I4,'.')
CALL END (2)
CALL END (3)
5 CALL END (1)
WRITE(LUWN,6)NC,NCMX
6 FORMAT(12X,' * * * THE NUMBER OF COLUMNS (=',I4,') OF THE MATRIX O
1F WHICH YOU ARE TRYING TO FIND THE RANK '/
312X,'MUST BE LESS THAN THE DIMENSION OF THE MATRIX (=',I4,' MINUS
42).')
CALL END (2)
CALL END (3)
END
SUBROUTINE XTX(X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('XTX')
NRMX=IO
NRMY=IO
CALL XXTX(NRMX,NRMY,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XXT(X,Y,NR,NC)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(IO,IO),Y(IO,IO)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('XXT')
NRMX=IO
NRMY=IO
CALL XXXT(NRMX,NRMY,X,Y,NR,NC)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE PSTCR(P,Q,C,R,NIN,NAT,MPRN)
C COMPUTES MEANS, STANDARD DEVIATIONS, COVARIANCE MATRIX, CORRELATION
C MATRIX FROM PRODUCT MATRIX AND ATTRIBUTE SUMS
C ROW DIMENSION OF ARRAYS FOR P, Q, C, AND R IS IO IN COMMON BLOCK B
C AS IN FORMAL, MUST BE AT LEAST NAT+2
C COLUMN DIMENSION OF ARRAYS FOR P, Q, C, AND R MUST BE AT LEAST NAT+2
C P CONTAINS PRODUCT MATRIX, INPUT
C ROW NAT + 1 CONTAINS ATTRIBUTE SUMS
C Q CONTAINS WITHIN SAMPLE SSCP MATRIX, ROW NAT+1 CONTAINS MEANS; OUTPUT
C C TO CONTAIN COVARIANCE MATRIX, OUTPUT
C ROW NAT+1 CONTAINS MEANS
C ROW NAT+2 CONTAINS STANDARD DEVIATIONS
C MAY USE SAME ARRAY AS P
C R TO CONTAIN CORRELATION MATRIX, OUTPUT
C MAY USE SAME ARRAY AS C (C NOT OUTPUT) OR P
C NIN IS NUMBER OF INDIVIDUALS IN SAMPLE
C NAT IS NUMBER OF ATTRIBUTES
C MPRN IS AN INTEGER INDICATOR FOR PRINTING OUTPUT
C = 0 FOR NO PRINTING; = 1 FOR PRINTED OUTPUT
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION P(IO,1),Q(IO,1),C(IO,1),R(IO,1)
COMMON/B/IO
COMMON/ZINCR/IZINCR
CALL QINCR ('PSTCR')
NRMP= IO
NRMQ= IO
NRMC= IO
NRMR= IO
CALL XPSTCR(NRMP,NRMQ,NRMC,NRMR,P,Q,C,R,NIN,NAT,MPRN)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XPSTCR(NRMP,NRMQ,NRMC,NRMR,P,Q,C,R,NIN,NAT,MPRN)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DIMENSION P(NRMP,1),Q(NRMQ,1),C(NRMC,1),R(NRMR,1)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
CALL QINCR ('XPSTCR')
NATP1 = NAT + 1
NATP2 = NAT + 2
CALL XDIMCH(NRMP,NATP1,'PSTCR',0)
CALL XDIMCH(NRMQ,NATP1,'PSTCR',0)
CALL XDIMCH(NRMC,NATP2,'PSTCR',0)
CALL XDIMCH(NRMR,NAT ,'PSTCR',0)
C MEANS, STANDARD DEVIATIONS, COVARIANCE MATRIX
NM1 = NIN - 1
DO 10 I = 1,NAT
Q(NATP1,I) = P(NATP1,I)/NIN
C(NATP1,I) = Q(NATP1,I)
DO 5 J = 1,I
Q(I,J) = P(I,J) - P(NATP1,I)*Q(NATP1,J)
Q(J,I) = Q(I,J)
C(I,J) = Q(I,J)/NM1
5 C(J,I) = C(I,J)
10 C(NATP2,I) = SQRT(C(I,I))
C PRINT MEANS, STANDARD DEVIATIONS, COVARIANCE MATRIX
IF(MPRN.EQ.0) GO TO 20
WRITE(LUWN,900)
900 FORMAT('-',17X,'MEAN',7X,'STANDARD DEVIATION')
901 FORMAT('0',I5,2F18.5)
DO 15 I = 1,NAT
15 WRITE(LUWN,901)I,C(NATP1,I),C(NATP2,I)
WRITE(LUWN,902)
902 FORMAT('-','COVARIANCE MATRIX')
CALL XPRINT(NRM,C,NAT,NAT,'F')
20 CONTINUE
C CORRELATION MATRIX
DO 30 I = 1,NAT
T = 1.0D0/C(NATP2,I)
DO 30 J = 1,I
R(I,J) = T*C(I,J)/C(NATP2,J)
30 R(J,I) = R(I,J)
C PRINT CORRELATION MATRIX
IF(MPRN.EQ.0) GO TO 40
WRITE(LUWN,903)
903 FORMAT('-','CORRELATION MATRIX')
CALL XPRINT(NRM,R,NAT,NAT,'F')
40 IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XABS(NRMX,NRMY,X,Y,NR,NC)
C TITLE: ABSOLUTE VALUE OF A MATRIX
C PROGRAMMER: ALLEN FLEISHMAN
C DATE:JULY, 1979
C LOCATION: LEDERLE LABS
C PURPOSE: ABSOLUTE VALUE OF MATRIX
C SYMBOLS:
C NRMX = NUMBER OF ROWS OF X IN CALLING PROGRAM
C NRMY = NUMBER OF ROWS OF Y IN CALLING PROGRAM
C X = INPUT MATRIX.
C Y = OUTPUT MATRIX.
C NR,NC = # OF ROWS AND COLUMNS, RESPECTIVELY, OF X AND Y.
C QINCR = PROGRAM COUNTER SUBROUTINE.
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE.
C I,J = SCRATCH VARIABLES.
C PROCEDURE: ELEMENTS OF X ARE MADE INTO ABSOLUTE VALUE
C AND THE RESULT IS STORED IN THE CORRESPONDING ELEMENTS OF Y.
C X OR Y MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON /ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XABS')
CALL XDIMCH(NRMX,NR,'ABS',0)
CALL XDIMCH(NRMY,NR,'ABS',0)
DO 1I=1,NR
DO 1J=1,NC
1 Y(I,J)=DABS(X(I,J))
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XADD(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
C TITLE: ADDITION OF 2 MATRICES
C PROGRAMMER: CARL FINKBEINER
C DATE: SEPTEMBER, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: ADDITION OF 2 MATRICES
C SYMBOLS:
C NRMX = NUMBER OF ROWS OF X IN CALLING PROGRAM
C NRMY = NUMBER OF ROWS OF Y IN CALLING PROGRAM
C NRMZ = NUMBER OF ROWS OF Z IN CALLING PROGRAM
C X = 1ST INPUT MATRIX.
C Y = 2ND INPUT MATRIX.
C Z = RESULTANT MATRIX.
C NR,NC = # OF ROWS AND COLUMNS, RESPECTIVELY, OF X, Y, AND Z.
C QINCR = PROGRAM COUNTER SUBROUTINE.
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE.
C I,J = SCRATCH VARIABLES.
C PROCEDURE: CORRESPONDING ELEMENTS OF X AND Y ARE ADDED
C AND THE RESLUT IS STORED IN THE CORRESPONDING ELEMENTS OF Z. ANY OF
C X, Y, OR Z MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),Z(NRMZ,1)
COMMON /ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XADD')
CALL XDIMCH(NRMX,NR,'ADD',0)
CALL XDIMCH(NRMY,NR,'ADD',0)
CALL XDIMCH(NRMZ,NR,'ADD',0)
C ADD X+Y AND STORE RESULT IN Z.
DO 1I=1,NR
DO 1J=1,NC
1 Z(I,J)=X(I,J)+Y(I,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XCHLSK(NRMX,NRMY,X,Y,N)
C TITLE: CHOLESKY DECOMPOSITION
C PROGRAMMER: CARL FINKBEINER
C DATE: AUGUST, 1975
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO OBTAIN A LOWER TRIANGULAR MATRIX Y, SUCH THAT YY' = X,
C A SYMMETRIC MATRIX.
C SYMBOLS:
C X = INPUT SYMMETRIC MATRIX
C Y = RESULTING LOWER TRIANGULAR MATRIX
C N = ORDER OF X AND Y
C A = LABELLED COMMON AREAS CONTAINING IPGM AND ISUBRU
C IPGM = PROGRAM COUNTER
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK ROUTINE
C I,J,K,I1,K1 = SCRATCH VARIABLES
C PROCEDURE: X IS CONVERTED TO A LOWER TRIANGULAR MATRIX BY TRANSFOR-
C MATION OF SUCCESSIVE ROWS OF ITS LOWER TRIANGLE. THE RESULT IS
C STORED IN Y. THE INPUT MATRIX MUST BE SYMMETRIC AND POSITIVE
C DEFINITE. X AND Y MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
C ALGORITHM REFERENCE: G.W.STEWART, INTRODUCTION TO MATRIX COMPUTATION
C ACADEMIC PRESS, 1973, PP.139-142.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION X(NRMX,1),Y(NRMY,1)
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON /ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XCHLSK')
CALL XDIMCH(NRMX,N,'CHLSK',0)
CALL XDIMCH(NRMY,N,'CHLSK',0)
C CHECK (1,1) ELEMENT; IF 0 OR NEGATIVE, ERROR
IF(X(1,1).LT.1.D-14)GOTO5
C FIND (1,1) ELEMENT
Y(1,1)=DSQRT(X(1,1))
IF(N.EQ.1)IZINCR=IZINCR - 1
IF(N.EQ.1)RETURN
C FIND (2,1) AND (2,2) ELEMENTS.
C ZERO OUT THE (1,2) ELEMENT.
Y(2,1)=X(2,1)/Y(1,1)
Y(2,2)=X(2,2)-Y(2,1)**2
IF(Y(2,2).LT.1.D-14)GOTO5
Y(2,2)=DSQRT(Y(2,2))
Y(1,2)=0.D0
IF(N.EQ.2)IZINCR=IZINCR - 1
IF(N.EQ.2)RETURN
C BEGIN LOOP FOR REMAINING ROWS OF LOWER TRIANGLE OF X
DO 1K=3,N
K1=K-1
C FIND THE FIRST ELEMENT OF THE ROW AND ZERO OUT THE TRANSPOSE ELEMENT
Y(K,1)=X(K,1)/Y(1,1)
Y(1,K)=0.D0
C LOOP FOR THE REST OF THE OFF-DIAGONAL ELEMENTS OF THE ROW; ZERO OUT TH
C TRANSPOSE ELEMENTS
DO 2I=2,K1
Y(K,I)=X(K,I)
I1=I-1
DO 3J=1,I1
3 Y(K,I)=Y(K,I)-Y(I,J)*Y(K,J)
Y(K,I)=Y(K,I)/Y(I,I)
2 Y(I,K)=0.D0
C FIND THE DIAGONAL ELEMENT OF THE ROW
Y(K,K)=X(K,K)
DO 4J=1,K1
4 Y(K,K)=Y(K,K)-Y(K,J)**2
C CHECK THE DIAGONAL ELEMENT; IF 0 OR NEGATIVE, ERROR
IF(Y(K,K).LT.1.D-14)GOTO5
1 Y(K,K)=DSQRT(Y(K,K))
IZINCR = IZINCR - 1
RETURN
5 CALL END (1)
WRITE(LUWN,6)
6 FORMAT ('0* * * * * THE SYMMETRIC MATRIX YOU PASSED TO THE CHLSK P
1ROGRAM ',/12X,'WAS NOT POSITIVE DEFINITE.'/12X,
2 'IF SOME, BUT NOT ALL, OF THE NUMBERS IN THIS MATRIX WERE VERY LA
3RGE, THIS COULD BE DUE TO MACHINE PRECISION.')
CALL END (2)
CALL END (3)
END
SUBROUTINE XCOLSM(NRMX,NRMY,X,Y,NR,NC)
C TITLE: COLUMN SUMMATION
C PROGRAMMER: CARL FINKBEINER
C DATE: AUGUST,1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO ADD THE ELEMENTS DOWN COLUMNS OF THE INPUT MATRIX AND
C STORE AS A ROW VECTOR.
C SYMBOLS:
C X = INPUT MATRIX WHOSE COLUMNS ARE TO BE SUMMED
C Y = OUTPUT MATRIX CONTAINING SUMS IN FIRST ROW
C NR = # OF ROWS OF X
C NC = # OF COLS OF X & Y
C NRMX = DIMENSIONS OF X IN CALLING PROGRAM
C NRMY = DIMENSIONS OF Y IN CALLING PROGRAM
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C PROCEDURE: THE COLUMNS OF X ARE SUMMED SIMULTANEOUSLY AND STORED IN Y.
C X & Y MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON /ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XCOLSM')
CALL XDIMCH(NRMX,NR,'COLSM',0)
CALL XDIMCH(NRMY, 1,'COLSM',0)
C STORE FIRST ROW OF X INTO FIRST ROW OF Y
DO 1I=1,NC
1 Y(1,I)=X(1,I)
IF(NR.EQ.1)IZINCR=IZINCR - 1
IF(NR.EQ.1)RETURN
C ADD ELEMENTS OF X DOWN THE COLUMNS SIMULTANEOUSLY AND STORE IN FIRST
C ROW OF Y.
DO 2I=2,NR
DO 2J=1,NC
2 Y(1,J)=Y(1,J)+X(I,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XCOPY(NRMX,NRMY,X,Y,NR,NC)
C TITLE: COPY
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE,1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO COPY ONE MATRIX INTO ANOTHER.
C SYMBOLS:
C X = MATRIX TO COPIED.
C Y = MATRIX INTO WHICH X IS COPIED.
C NR = # OF ROWS OF X AND Y.
C NC = # OF COLS OF X AND Y.
C NRMX = DIMENSIONALITY OF X.
C NRMY = DIMENSIONALITY OF Y.
C QINCR = PROGRAM COUNTING SUBROUTINE.
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE.
C I,J = SCRATCH VARIABLE.
C PROCEDURE: MATRIX X IS ASSIGNED TO MATRIX Y.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XCOPY')
CALL XDIMCH(NRMX,NR,'COPY',0)
CALL XDIMCH(NRMY,NR,'COPY',0)
C COPY X INTO Y
DO 1I=1,NR
DO 1J=1,NC
1 Y(I,J)=X(I,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XDET(NRMX,NCMX,NRMY,NCMY,X,D,N,Y)
C TITLE: DETERMINANT
C PROGRAMMER: CARL FINKBEINER
C DATE: AUGUST, 1975
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: FIND THE DETERMINANT OF A SQUARE , REAL MATRIX
C SYMBOLS:
C X = INPUT SQUARE MATRIX
C D = DETERMINANT OF X
C N = ORDER OF X
C Y = SCRATCH ARRAY
C A = LABELLED COMMON AREAS CONTAINING IPGM AND ISUBRU
C IPGM = PROGRAM COUNTER
C NRMX = ROW DIMENSION OF X IN THE CALLING PROGRAM
C NRMY = ROW DIMENSION OF Y IN THE CALLING PROGRAM
C NCMX = COLUMN DIMENSION OF X IN THE CALLING PROGRAM
C NCMY = COLUMN DIMENSION OF Y IN THE CALLING PROGRAM
C QINCR = PROGRAM COUNTER ROUTINE
C XDIMCH = DIMENSIONALITY CHECK ROUTINE
C QWRGDM = ERROR MESSAGE ROUTINE
C IER = ERROR FLAG USED BY QMINVZ
C IDET = EXPONENT OF DETERMINANT RETURNED FROM QMINVZ
C QMINVZ = MATRIX INVERSION ROUTINE WHICH RETURNS DETERMINANT
C PROCEDURE: THE DETERMINANT IS FOUND BY QMINVZ DURING THE PROCESS OF
C INVERSION. THUS, IF D IS NOT 0, Y CONTAINS THE INVERSE OF X ON
C OUTPUT. X AND Y MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
C NOTICE THAT X WILL BE DESTROYED IF THIS IS SO. N MUST BE STRICTLY
C LESS THAN NRMX.
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),D
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON /ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XDET')
CALL XDIMCH(NRMX,N,'DET',0)
CALL XDIMCH(NRMY,N,'DET',0)
IF(N.EQ.1)CALL QWRGDM(N,-9999)
I01=NCMY-1
IF(N.GE.I01)GOTO1
C FIND INVERSE OF X AND ITS DETERMINANT
IER=1
IDET=1
CALL QMINVZ(NRMX,NRMY,NRMY,X,N,DET,IDET,Y,Y(1,I01),IER,X(1,NCMX))
IF(IER.GT.0)D=0.D0
D=DET*(10.D0**IDET)
IZINCR = IZINCR - 1
RETURN
1 CALL END (1)
WRITE(LUWN,2)N,NCMY
2 FORMAT('0* * * * * TO FIND THE DETERMINANT OF A MATRIX THE ORDER O
1F THE MATRIX (=',I4,') MUST BE LESS THAN THE DIMENSIONS (=',I4,
2').')
CALL END (2)
CALL END (3)
END
SUBROUTINE QGELGZ(NRMB,NRMAA,B,AA,N,NR,X,A,EPS,IER)
C-----------------------------------------------------------------------
C SOLVES SIMULTANEOUS LINEAR EQUATIONS WITH >1 RIGHT HAND SIDE BY
C GAUSS ELIMINATION WITH COMPLETE PIVOTING.
C PARAMETERS ARE:
C B - N BY NR MATRIX OF RIGHT HAND SIDES
C AA - N BY N COEFFICIENT MATRIX
C N - NUMBER OF EQUATIONS (N.LE.NRMB.AND.N.LE.NRMAA)
C NR - NUMBER OF RIGHT HAND SIDES
C X - SOLUTION IS STORED IN THE 1'ST N*NR LOCATIONS OF X. IT IS
C PERMISSIBLE FOR X AND B TO SHARE THE SAME STORAGE IN WHICH
C CASE B IS OVERWRITTEN BY THE SOLUTION. OTHERWISE, B IS RE-
C TAINED IN THE COMPUTATION.
C A - SCRATCH ARRAY USED BY QGELGZ. IT IS PERMISSIBLE FOR A AND AA
C TO SHARE THE SAME STORAGE IN WHICH CASE AA IS DESTROYED BY
C THE COMPUTATION. OTHERWISE AA IS RETAINED IN THE COMPUTATION
C EPS - INPUT CONSTANT USED IN TEST ON LOSS OF SIGNIFICANT DIGITS
C SUGGESTED VALUE =1.D-15
C IER - ON INPUT IF IER:
C = 0 NO MESSAGES WILL BE PRINTED.
C = 1 MESSAGES WILL BE PRINTED.
C ON OUTPUT IF IER:
C = 0 NO ERROR DETECTED.
C =-1 N.LE.1, NR.LT.1, OR A PIVOT ELEMENT = 0.D0 .
C NO SOLUTION IS COMPUTED BY QGELGZ.
C = K QGELGZ DETECTED A LOSS OF SIGNIFICANT DIGITS IN THE
C COMPUTED SOLUTION. THE K'TH PIVOT ELEMENT (K=1,N-1) WAS
C LESS THAN EPS TIMES THE LARGEST ELEMENT OF THE INPUT
C MATRIX. THE INPUT MATRIX IS PROBABLY SINGULAR.
C-----------------------------------------------------------------------
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION B(NRMB,1),AA(NRMAA,1),X(1),A(1)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
CALL QINCR ('QGELGZ')
IER1=IER
IER=0
NB=N*2
C-----TEST FOR ERROR IN INPUT FOR N AND NR.
IF(N.LE.1)GOTO150
IF(NR.LT.1)GOTO150
C-----MAKES THE B MATRIX 1-DIMENSIONAL.
K=0
DO 10J=1,NR
DO 10I=1,N
K=K+1
10 X(K)=B(I,J)
C-----MAKES THE AA MATRIX 1-DIMENSIONAL.
K=0
DO 20J=1,N
DO 20I=1,N
K=K+1
20 A(K)=AA(I,J)
C-----SEARCH FOR GREATEST ELEMENT IN MATRIX A.
PIV=0.D0
MM=N*N
NM=NR*N
DO 30L=1,MM
TB=DABS(A(L))
IF(TB.LE.PIV)GOTO30
PIV=TB
I=L
30 CONTINUE
TOL=EPS*PIV
C-----A(I) IS PIVOT ELEMENT. PIV CONTAINS THE ABSOLUTE VALUE OF A(I).
C-----START ELIMINATION.
LST=1
DO 110K=1,N
C-----TEST ON SINGULARITY.
IF(PIV.EQ.0.D0)GOTO150
IF(IER.NE.0)GOTO40
IF(PIV.GT.TOL)GOTO40
IF(IER1.EQ.1) CALL END (1)
IF(IER1.EQ.1)WRITE(LUWN,2001)K
2001 FORMAT(1X,'**** WARNING ****'/1X,'LOSS OF SIGNIFICANT DIGITS IN TH
*E SOLUTION PRODUCED BY QGELGZ.'/1X,'THE PIVOT ELEMENT AT STEP ',I3
*,' WAS LESS THAN EPS TIMES THE LARGEST ELEMENT OF THE INPUT MATRIX
*.'/1X,'THE INPUT MATRIX IS PROBABLY SINGULAR.')
CALL END (2)
IER=K
40 PIVI=1.D0/A(I)
J=(I-1)/N
I=I-J*N-K
J=J+1-K
C-----ROW REDUCTION AND ROW INTERCHANGE IN RIGHT HAND SIDE MATRIX.
C-----DIVIDE PIVOT ROW ELEMENTS BY PIVOT ELEMENT AND MAKE PIVOT ROW
C-----THE K'TH ROW.
DO 50L=K,NM,N
LL=L+I
TB=PIVI*X(LL)
X(LL)=X(L)
50 X(L)=TB
C-----COLUMN INTERCHANGE IN MATRIX A.
IF(K.EQ.N)GOTO120
LEND=LST+N-K
C-----IF J.LE.0 THEN THE PIVOT ELEMENT IS IN THE K'TH COLUMN SO DON'T
C-----NEED TO INTERCHANGE COLUMNS.
IF(J.LE.0)GOTO70
II=J*N
DO 60L=LST,LEND
TB=A(L)
LL=L+II
A(L)=A(LL)
60 A(LL)=TB
C-----ROW INTERCHANGE AND PIVOT ROW REDUCTION IN MATRIX A.
70 DO 80L=LST,MM,N
LL=L+I
TB=PIVI*A(LL)
A(LL)=A(L)
80 A(L)=TB
C-----SAVE COLUMN INTERCHANGE INFORMATION ALONG THE MAIN DIAGONAL.
A(LST)=J
C-----ELEMENT REDUCTION IN A AND NEXT PIVOT SEARCH.
PIV=0.D0
C-----LST IS INDEX OF (K,K) DIAGONAL ELEMENT. THE PIVOT ELEMENT IS ALSO
C-----ON MAIN DIAGONAL SO NEED TO ADD 1 TO LST TO GET AT 1ST NON-PIVOT
C-----ELEMENT OF ROW.
LST=LST+1
J=0
C-----II IS THE ROW NUMBER WE ARE WORKING ON.
DO 100II=LST,LEND
PIVI=-A(II)
IST=II+N
J=J+1
DO 90L=IST,MM,N
LL=L-J
A(L)=A(L)+PIVI*A(LL)
C-----SEARCH FOR NEXT PIVOT IN 'A' MINUS THE FIRST K COLUMNS AND ROWS.
TB=DABS(A(L))
IF(TB.LE.PIV)GOTO90
PIV=TB
I=L
90 CONTINUE
C-----ELEMENT REDUCTION IN RIGHT HAND SIDE MATRIX.
DO 100L=K,NM,N
LL=L+J
100 X(LL)=X(LL)+PIVI*X(L)
C-----MAKE LST POINT AT THE NEXT DIAGONAL ELEMENT.
110 LST=LST+N
C-----END OF ELIMINATION LOOP.
C-----BACK SUBSTITUTION AND BACK INTERCHANGE.
120 IST=MM+N
LST=N+1
DO 140I=2,N
II=LST-I
IST=IST-LST
L=IST-N
C-----MAKE SURE YOU GET THE SAME INTEGER BACK WHICH YOU STORED.
L=A(L)+.5D0
DO 140J=II,NM,N
TB=X(J)
LL=J
DO 130K=IST,MM,N
LL=LL+1
130 TB=TB-A(K)*X(LL)
K=J+L
X(J)=X(K)
140 X(K)=TB
IZINCR = IZINCR - 1
RETURN
C-----ERROR RETURN.
150 IF(IER1.EQ.1) CALL END (1)
IF(IER1.EQ.1)WRITE(LUWN,2002)
2002 FORMAT(1X,'NO SOLUTION COMPUTED BY QGELGZ BECAUSE N.LE.1, NR.LT.1,
* OR A PIVOT ELEMENT WAS EQUAL TO 0.D0 .')
CALL END (2)
IER=-1
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XDIAG(NRMX,NRMY,X,Y,N)
C TITLE: DIAGONAL MATRIX
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO TAKE THE DIAGONAL OF A MATRIX AND CREATE A DIAGONAL MATRIX
C SYMBOLS:
C X = INPUT MATRIX
C Y = OUTPUT DIAGONAL MATRIX = DIAG.(X)
C N = # OF ROWS AND COLS OF X AND Y
C NRMX = DIMENSIONS OF X IN THE CALLING PROGRAM
C NRMY = DIMENSIONS OF Y IN THE CALLING PROGRAM
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C N1,I,J,K = SCRATCH VARIABLES
C PROCEDURE: FOR A GIVEN ROW AND COLUMN OF X, THE CORRESPONDING OFF-DIAG
C ARE SET TO 0 IN Y, AND THE CORRESPONDING ELEMENT OF Y IS SET EQUAL
C TO THE ELEMENT OF X. X AND Y CAN BE THE SAME MATRIX IN THE CALLING
C PROGRAM. NOTICE THAT X AND Y ARE ASSUMED TO BE SQUARE AND LOCATED
C IN THE UPPER LEFT OF THE CORRESPONDING MATRICES IN THE CALLING
C PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XDIAG')
CALL XDIMCH(NRMX,N,'DIAG',0)
CALL XDIMCH(NRMY,N,'DIAG',0)
C ZERO OUT OFF DIAGONALS OF Y AND RETAIN DIAG(X) IN DIAG(Y)
N1=N-1
DO 1I=1,N1
J=I+1
DO 2K=J,N
Y(I,K)=0.D0
2 Y(K,I)=0.D0
1 Y(I,I)=X(I,I)
Y(N,N)=X(N,N)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XEIGEN(NRMR,NRMV,NRMD,R,V,D,N,NVE)
C TITLE: EIGEN SOLUTION
C PROGRAMMER: CARL FINKBEINER
C DATE: AUGUST, 1975
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: FIND THE EIGEN SOLUTION FOR A SYMMETRIC MATRIX
C SYMBOLS:
C R = INPUT SYMMETRIC MATRIX
C V = EIGENVECTORS
C D = DIAGONAL MATRIX OF EIGENVALUES
C N = ORDER OF R
C NVE = # OF EIGENVECTORS REQUESTED BY USER
C A = LABELLED COMMON AREAS CONTAINING IPGM AND ISUBRU
C NRMR = ORDER OF R IN CALLING PROGRAM
C NRMV = ORDER OF V IN CALLING PROGRAM
C NRMD = ORDER OF D IN CALLING PROGRAM
C IPGM = PROGRAM COUNTER
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK ROUTINE
C IERR = ERROR FLAG USED BY QLRSVM
C I,J,K = SCRATCH VARIABLES
C PROCEDURE:
C SUBROUTINE QLRSVM IS INVOKED TO DO THE REAL NUMBER CRUNCHING
C AFTER SUCCESSFUL EXECUTION OF THAT ROUTINE, ALL OF THE EIGENVALUES
C ARE STORED IN THE DIAGONAL ELEMENTS OF D AND THE OFF-DIAGONALS ARE
C ZEROED OUT. IF NVE IS .GT. N, IT IS CHANGED TO N; IF NVE IS .LT. 1
C IT IS CHANGED TO 1. V CONTAINS ONLY NVE COLUMNS ON OUT, ONE FOR
C EACH EIGENVECTOR. THE ORDER OF THE INPUT MATRICES (I.E., NRMR,ETC.)
C MUST BE .GE.N
C R, V, AND D MUST BE UNIQUE MATRICES IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION R(NRMR,1),V(NRMV,1),D(NRMD,1)
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XEIGEN')
CALL XDIMCH(NRMR,N,'EIGEN',0)
CALL XDIMCH(NRMV,N,'EIGEN',0)
CALL XDIMCH(NRMD,N,'EIGEN',0)
IF (NRMD.LT.10) GOTO 1
NVE1=NVE
IF(NVE1.GT.N)NVE1=N
IF(NVE1.LT.1)NVE1=1
C EIGEN SOLUTION
CALL QLRSVM(NRMR,NRMV,NRMD,N,NVE1,R,V,D,D(1,2),D(1,3),D(1,4),
+ D(1,5),D(1,6),IERR)
IF(IERR.NE.0)GOTO2
C CREATE D AS A DIAGONAL MATRIX
DO 3I=2,N
D(I,I)=D(I,1)
J=I-1
DO 3K=1,J
D(I,K)=0.D0
3 D(K,I)=0.D0
IZINCR = IZINCR - 1
RETURN
C ERROR MESSAGES
1 CALL END (1)
WRITE(LUWN,4)
4 FORMAT('0* * * * * YOUR CALL TO THE EIGEN PROGRAM HAS FAILED.')
WRITE(LUWN,5)
5 FORMAT(12X,'THE ORDER OF THE INPUT MATRICES AS SPECIFIED IN THE DI
1MENSION STATEMENT CANNOT BE LESS THAN 10.')
CALL END (2)
CALL END (3)
2 CALL END (1)
WRITE(LUWN,4)
WRITE(LUWN,7)
7 FORMAT(12X,'SOMETHING IS PROBABLY WRONG WITH THE INPUT SYMMETRIC M
1ATRIX.')
CALL END (2)
CALL END (3)
END
SUBROUTINE XELEDI(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC,ICHECK)
C TITLE: ELEMENTWISE DIVISION
C PROGRAMMER: ALLEN FLEISHMAN
C DATE: JULY, 1979
C LOCATION: LEDERLE LABS
C PURPOSE: TO FIND THE ELEMENTWISE DIVISION OF 2 MATRICES.
C SYMBOLS:
C X = 1ST INPUT MATRIX
C Y = 2ND INPUT MATRIX
C Z = RESULT MATRIX
C NR = # OF ROWS OF X,Y, AND Z
C NC = # OF COLS OF X,Y, AND Z
C ICHECK = FLAG TO DETERMINE THE EFFECT OF A ZERO ENTRY IN Y.
C = 0 (ON INPUT) STOP PROGRAM WHEN ZERO IS FOUND. (DEFAULT)
C = 1 (ON INPUT) CONTINUE BUT SET Z(I,J) TO 0.0D0.
C NRMX = NUMBER OF ROWS OF X IN CALLING PROGRAM
C NRMY = NUMBER OF ROWS OF Y IN CALLING PROGRAM
C NRMZ = NUMBER OF ROWS OF Z IN CALLING PROGRAM
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C I,J = SCRATCH VARIABLES
C PROCEDURE: DIVIDE AN ELEMENT OF X BY THE CORRESPONDING ELEMENT OF Y
C STORE IN THE CORRESPONDING ELEMENT OF Z. X, Y, AND Z MUST HAVE THE
C SAME DIMENSIONS. X, Y, AND Z, OR ANY PAIR COMBINATION, MAY BE THE
C SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),Z(NRMZ,1)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XELEDI')
CALL XDIMCH(NRMX,NR,'ELEDI',0)
CALL XDIMCH(NRMY,NR,'ELEDI',0)
CALL XDIMCH(NRMZ,NR,'ELEDI',0)
C FIND ELEMENTWISE DIVISION
DO 1I=1,NR
DO 1J=1,NC
IF(DABS(Y(I,J)).LT.1.D-10)GOTO 2
IF (ICHECK.GT.1.OR.ICHECK.LT.0) ICHECK = 0
IF(IFLAG.EQ.0) CALL END (1)
IF(IFLAG.EQ.0) WRITE(LUWN,100) Y(I,J), I, J
100 FORMAT (25X,' YOU TRIED TO DIVIDE AN ELEMENT BY ',D20.15,
3 ' IT WAS ELEMENT ',I5,',',I5)
IF(IFLAG.EQ.0)CALL END (2)
IF(IFLAG.EQ.0)CALL END (3)
IF(IFLAG.EQ.1)Z(I,J)=0.0D0
IF(IFLAG.EQ.1)GOTO 1
2 Z(I,J)=X(I,J)/Y(I,J)
1 CONTINUE
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XELEML(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
C TITLE: ELEMENTWISE MULTIPLICATION
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO FIND THE ELEMENTWISE PRODUCT OF 2 MATRICES.
C SYMBOLS:
C X = 1ST INPUT MATRIX
C Y = 2ND INPUT MATRIX
C Z = RESULT MATRIX
C NR = # OF ROWS OF X,Y, AND Z
C NC = # OF COLS OF X,Y, AND Z
C NRMX = NUMBER OF ROWS OF X IN CALLING PROGRAM
C NRMY = NUMBER OF ROWS OF Y IN CALLING PROGRAM
C NRMZ = NUMBER OF ROWS OF Z IN CALLING PROGRAM
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C I,J = SCRATCH VARIABLES
C PROCEDURE: MULTIPLY AN ELEMENT OF X BY THE CORRESPONDING ELEMENT OF Y
C STORE IN THE CORRESPONDING ELEMENT OF Z. X, Y, AND Z MUST HAVE THE
C SAME DIMENSIONS. X, Y, AND Z, OR ANY PAIR COMBINATION, MAY BE THE
C SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),Z(NRMZ,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XELEML')
CALL XDIMCH(NRMX,NR,'ELEML',0)
CALL XDIMCH(NRMY,NR,'ELEML',0)
CALL XDIMCH(NRMZ,NR,'ELEML',0)
C FIND ELEMENTWISE PRODUCT
DO 1I=1,NR
DO 1J=1,NC
1 Z(I,J)=X(I,J)*Y(I,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XGENR8(NRMX,C,X,NR,NC)
C TITLE: GENERATE
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE8 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO GENERATE A MATRIX CONTAINING CONSTANTS
C SYMBOLS:
C C = THE CONSTANT (REAL)
C X = THE MATRIX TO BE GENERATED CONTAINING ALL C VALUES
C NR,NC = # OF ROWS AND COLS OF X
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C I,J = SCRATCH VARIABLES
C PROCEDURE: EVERY ELEMENT OF X IS SET EQUAL TO C. X MAY BE ANY SIZED
C MATRIX. C MAY BE A REAL*4 OR DOUBLE PRECISION NUMBER.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XGENR8')
CALL XDIMCH(NRMX,NR,'GENR8',0)
C FILL EVERY ELEMENT OF X WITH C.
DO 1I=1,NC
DO 1J=1,NR
1 X(J,I)=C
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XGINVR(NRMX,NRMY,NRMA,NRMD,NRMC,X,Y,NR,NC,A,D,C)
C TITLE: GENERALIZED INVERSE
C PROGRAMMER: CARL FINKBEINER
C DATE: MAY, 1975
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: COMPUTE THE GENERALIZED INVERSE (Y) OF A MATRIX X,
C SUCH THAT XYX = X
C SYMBOLS:
C X = INPUT MATRIX
C Y = GENERALIZED INVERSE OF X
C NR,NC = # OF ROWS & COLS OF X = # OF COLS & ROWS OF Y
C A,D,C = TEMPORARY WORK MATRICES
C NRMX,NRMY,NRMA,NRMD,NRMC = NUMBER OF ROWS OF X, Y, A, D, AND C
C RESPECTIVELY IN CALLING PROGRAM. MUST BE .GE. 10.
C QINCR = PROGRAM COUNTER ROUTINE
C XDIMCH = DIMENSIONALITY CHECK ROUTINE
C QLRSVM = IMPLICIT QL EIGEN ROUTINE
C PROCEDURE: THE MATRIX X'X IS FOUND AND ITS EIGENVALUES (D) AND
C EIGENVECTORS (A) ARE COMPUTED. Y IS THEN SIMPLY A(D**-1)A'X'.
C Y MAY BE THE SAME AS EITHER A OR D IN THE CALLING PROGRAM,
C BUT NO OTHER COMBINATION (INCLUDING A AND D) MAY BE THE SAME.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION X(NRMX,1),Y(NRMY,1),A(NRMA,1),D(NRMD,1),C(NRMC,1)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XGINVR')
CALL XDIMCH(NRMX,NR,'GINVR',0)
CALL XDIMCH(NRMC,NC,'GINVR',0)
CALL XDIMCH(NRMD,NC,'GINVR',0)
CALL XDIMCH(NRMA,NC,'GINVR',0)
CALL XDIMCH(NRMY,NC,'GINVR',0)
I0=NRMX
IF(NRMX.LT.10)GOTO1
I0=NRMY
IF(NRMY.LT.10)GOTO1
I0=NRMA
IF(NRMA.LT.10)GOTO1
I0=NRMD
IF(NRMD.LT.10)GOTO1
I0=NRMC
IF(NRMC.LT.10)GOTO1
C FIND X'X
DO 4I=1,NC
DO 4K=1,I
C(I,K)=0.D0
DO 5J=1,NR
5 C(I,K)=C(I,K)+X(J,I)*X(J,K)
4 C(K,I)=C(I,K)
C FIND EIGENVALUES AND EIGENVECTORS OF X'X
CALL QLRSVM(NRMC,NRMA,NRMD,NC,NC,C,A,D,D(1,2),D(1,3),D(1,4),
+ D(1,5),D(1,6),IERR)
IF(IERR.NE.0)GOTO2
C FIND # OF NON-ZERO EIGENVALUES OF X'X
DO 6J=1,NC
IF(D(J,1).LE.1.D-7)GOTO7
6 CONTINUE
J=NC
C COMPUTE A(D**-1)A' AND STORE IN C
7 DO 8I=1,NC
DO 8K=1,I
C(I,K)=0.D0
DO 9L=1,J
9 C(I,K)=C(I,K)+A(I,L)*A(K,L)/D(L,1)
8 C(K,I)=C(I,K)
C COMPUTE Y = CX'
DO 10I=1,NC
DO 10J=1,NR
Y(I,J)=0.D0
DO 10K=1,NC
10 Y(I,J)=Y(I,J)+C(I,K)*X(J,K)
IZINCR = IZINCR - 1
RETURN
1 CALL END (1)
WRITE(LUWN,3)I0
3 FORMAT ('0* * * * * YOUR CALL TO THE GINVR PROGRAM HAS FAILED BECA
1USE THIS PROGRAM REQUIRES NO LESS THAN 10 COLUMNS FOR THE DIMENSIO
2NS OF'/12X,'THE MATRICES. YOU SPECIFIED',I4,'.')
CALL END (2)
CALL END (3)
2 CALL END (1)
WRITE(LUWN,12)
12 FORMAT('0* * * * * SOMETHING WAS VERY WRONG WITH THE NUMBERS PASSE
1D TO THE GINVR PROGRAM IN THE INPUT MATRIX.')
CALL END (2)
CALL END (3)
END
SUBROUTINE XHORIZ(NRMX,NRMY,NRMZ,X,Y,Z,NR,NCX,NCY)
C TITLE: HORIZONTAL CONCATENATE
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE,1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO CONCATENATE THE ROWS OF 2 MATRICES TO CREATE A SUPER
C MATRIX.
C SYMBOLS:
C X = 1ST INPUT MATRIX
C Y = 2ND INPUT MATRIX
C Z = RESULTING SUPER MATRIX = (X:Y)
C NR = # OF ROWS IN BOTH X AND Y
C NCX = # OF COLS IN X
C NCY = # OF COLS IN Y
C NRMX, NRMY, AND NRMZ = DIMENSIONS OF X, Y, AND Z IN CALLING PROGRAM
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C NC = # OF COLS IN Z = NCX + NCY
C I,J,K = SCRATCH VARIABLES
C PROCEDURE: A MATRIX Z = (X:Y) IS CREATED BY STORING X IN COLS 1 THRU
C NCX, ROWS 1 THRU NR OF Z AND THEN BY STORING Y IN COLS (NCX+1) THRU
C NC AND ROW 1 THRU NR OF Z. THUS Z WILL BE NR X NC. ALL THREE
C MATRICES,X,Y, AND Z MAY BE THE SAME MATRIX IN THE CALLING PROGRAM;
C MATRICES X AND Y MAY BE THE SAME BUT DIFFERENT FROM Z; MATRICES X
C AND Z MAY BE THE SAME BUT DIFFERENT FROM Y; HOWEVER, MATRICES Y AND
C Z MAY NOT BE THE SAME IF X IS DIFFERENT IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),Z(NRMZ,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XHORIZ')
CALL XDIMCH(NRMX,NR,'HORIZ',0)
CALL XDIMCH(NRMY,NR,'HORIZ',0)
NC=NCX+NCY
CALL XDIMCH(NRMZ,NC,'HORIZ',1)
C LOOP FOR SETTING LEFT MOST PORTION OF Z EQUAL TO X.
DO 2J=1,NCX
DO 2K=1,NR
2 Z(K,J)=X(K,J)
C LOOP FOR SETTING NEXT PORTION TO THE RIGHT IN Z EQUAL TO Y.
I=NCX
DO 3J=1,NCY
I=I+1
DO 3K=1,NR
3 Z(K,I)=Y(K,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XIDENT(NRMX,X,N)
C TITLE: IDENT
C PROGRAMMER: ALLEN I. FLEISHMAN
C DATE: DECEMBER, 1978
C LOCATION: LEDERLE LABS
C PURPOSE: CREATE IDENTITY MATRIX OF ORDER N
C SYMBOLS:
C N = ORDER OF IDENTITY MATRIX
C X = IDENTITY MATRIX
C A = LABELED COMMON AREAS CONTAINING IPGM AND ISUBRU
C IPGM = PROGRAM COUNTER
C NRMX = DIMENSION OF X IN MAIN PROGRAM
C QINCR = PROGRAM COUNTER ROUTINE
C XDIMCH = DIMENSIONALITY CHECK ROUTINE
C PROCEDURE: CREATES A SQUARE IDENTITY MATRIX OF ORDER N
C IF N > 0
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1)
COMMON /A/IPGM,ISUBRU
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XIDENT')
CALL XDIMCH (NRMX,N,'IDENT',0)
DO 1 I=1,N
DO 2 J=1,N
X(I,J)=0.0D0
2 CONTINUE
X(I,I)=1.0D0
1 CONTINUE
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XINPUT(NRMX,X)
C TITLE: DATA INPUT FROM CARDS
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO INPUT DATA FROM CARDS.
C SYMBOLS:
C X = MATRIX IN WHICH TO STORE INPUT DATA
C FMT = ARRAY FOR CONTAINING INPUT FORMAT (READ FROM THE FIRST CARD OF
C DATA DECK).
C C = LABELED COMMON AREAS CONTAINING NMAT, AND LURN.
C NRMX = DIMENSIONALITY OF X.
C NMAT = MATRIX COUNTER.
C LURN = DATA SET REFERENCE NUMBER
C QINCR = PROGRAM NUMBER COUNTER SUBROUTINE.
C NR = NUMBER OF ROWS INPUT TO X.
C NC = NUMBER OF COLS INPUT TO X.
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE.
C PROCEDURE: THIS SUBROUTINE READS THE HEADER RECORD FROM THE DATA FILE,
C INCREMENTS THE MATRIX COUNTER, PRINTS OUT NMAT,LURN,NR,NC,& FMT FOR
C USER, AND, IF ALL IS OK, READS THE ACTUAL DATA FILE. IF, FOR ANY
C REASON THE END OF THE FILE IS ENCOUNTERED BEFORE EXPECTED,
C PROCESSING TERMINATE AFTER PRINTING AN APPROPRIATE MESSAGE.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),FMT(9),XIN,OUTPUT
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XIN(5),OUTPUT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XINPUT')
C READ DATA DECK HEADER RECORD
READ(LURN,1,END=3)NR,NC,FMT
1 FORMAT(2I4,9A8)
C INCREMENT MATRIX COUNTER AND PRINT SPECIFICATIONS
NMAT=NMAT+1
WRITE(LUWN,4)NMAT,LURN,NR,NC,FMT,XIN(1)
4 FORMAT('0SPECIFICATIONS FOR INPUT MATRIX #',I3,'(LURN =',I3,'):'/
1 3X,' NUMBER OF ROWS =',I4,'; NUMBER OF COLS =',I4,'; FORMAT =
2',9A8/'0DATA READ FROM DATA FILE ', A10)
C MORE BOOKKEEPING
CALL XDIMCH(NRMX,NR,'INPUT',0)
C READ ACTUAL DATA FILE
DO 2I=1,NR
READ(LURN,FMT,END=3)(X(I,J),J=1,NC)
2 CONTINUE
C NORMAL RETURN
IZINCR = IZINCR - 1
RETURN
C IF END OF FILE ON LURN OCCURS BEFORE EXPECTED, PRINT ERROR MESSAGE
3 CALL END (1)
WRITE(LUWN,5)
5 FORMAT('0* * * * * YOU HAVE RUN OUT OF INPUT DATA RECORDS.THIS COU
1LD HAVE OCCURRED BECAUSE OF INCORRECT FORMAT OR DIMENSION SPECIFIC
2ATIONS.')
CALL END (2)
CALL END (3)
STOP
END
SUBROUTINE XTINPT(NRMX,XIN,X)
C TITLE: TEMPORARY DATA INPUT FROM DATA UNIT 25
C PROGRAMMER: ALLEN I. FLEISHMAN
C DATE: AUGUST 1979
C LOCATION: LEDERLE LABS
C PURPOSE: TO INPUT DATA FROM TEMPORARY DATA UNIT 25.
C SYMBOLS:
C XIN = NAME OF DATA FILE WHICH CONTAINS THE DATA TO BE READ.
C X = MATRIX IN WHICH TO STORE INPUT DATA.
C FMT = ARRAY FOR CONTAINING INPUT FORMAT (READ FROM THE FIRST CARD OF
C DATA DECK).
C C = LABELED COMMON AREAS CONTAINING LURN, LUWN, NMAT, XINPUT AND OUTPUT.
C NRMX = DIMENSIONALITY OF X.
C NMAT = MATRIX COUNTER.
C QINCR = PROGRAM NUMBER COUNTER SUBROUTINE.
C NR = NUMBER OF ROWS INPUT TO X.
C NC = NUMBER OF COLS INPUT TO X.
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE.
C PROCEDURE: THIS SUBROUTINE OPENS UP UNIT 4 A FILE WHICH IS CALLED
C 'XIN', THEN READS THE HEADER RECORD FROM THE DATA FILE,
C INCREMENTS THE MATRIX COUNTER, PRINTS OUT NMAT,LURN,NR,NC,& FMT FOR
C USER, AND, IF ALL IS OK, READS THE ACTUAL DATA FILE. IF, FOR ANY
C REASON THE END OF THE FILE IS ENCOUNTERED BEFORE EXPECTED,
C PROCESSING TERMINATE AFTER PRINTING AN APPROPRIATE MESSAGE.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),FMT(9),XIN(5)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),OUTPUT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XTINPT')
C READ DATA DECK HEADER RECORD
C OPEN UNIT 4 WHICH WILL CONTAIN THE FILE 'XIN'
OPEN (UNIT=LUEX,DIALOG=XIN,ACCESS='SEQIN')
READ(LUEX,1,END=3)NR,NC,FMT
1 FORMAT(2I4,9A8)
C INCREMENT MATRIX COUNTER AND PRINT SPECIFICATIONS
NMAT=NMAT+1
WRITE(LUWN,4)NMAT,LUEX,NR,NC,FMT,XIN(1)
4 FORMAT('0SPECIFICATIONS FOR INPUT MATRIX #',I3,'(LUEX =',I3,'):'/
1 3X,' NUMBER OF ROWS =',I4,'; NUMBER OF COLS =',I4,'; FORMAT =
2',9A8/'0DATA READ FROM DATA FILE ',A10)
C MORE BOOKKEEPING
CALL XDIMCH(NRMX,NR,'TINPT',0)
C READ ACTUAL DATA FILE
DO 2I=1,NR
READ(LUEX,FMT,END=3)(X(I,J),J=1,NC)
2 CONTINUE
C NORMAL RETURN
CLOSE (UNIT=LUEX,DIALOG=XIN)
IZINCR = IZINCR - 1
RETURN
C IF END OF FILE ON LUEX OCCURS BEFORE EXPECTED, PRINT ERROR MESSAGE
3 CALL END (1)
WRITE(LUWN,5)
5 FORMAT('0* * * * * YOU HAVE RUN OUT OF INPUT DATA RECORDS.THIS COU
1LD HAVE OCCURRED BECAUSE OF INCORRECT FORMAT OR DIMENSION SPECIFIC
2ATIONS.')
CALL END (2)
CALL END (3)
STOP
END
SUBROUTINE XINVRS(NRMX,NCMX,NRMY,NCMY,X,Y,N)
C TITLE: INVERSE
C PROGRAMMER: CARL FINKBEINER
C DATE: AUGUST, 1975
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: INVERT A SQUARE , REAL MATRIX
C SYMBOLS:
C X = INPUT SQUARE MATRIX
C Y = RESULTING INVERSE OF X
C N = ORDER OF X
C A = LABELLED COMMON AREAS CONTAINING IPGM AND ISUBRU
C IPGM = PROGRAM COUNTER
C NRMX & NRMY = DIMENSIONS OF X AND Y IN MAIN PROGRAM
C QINCR = PROGRAM COUNTER ROUTINE
C XDIMCH = DIMENSIONALITY CHECK ROUTINE
C QWRGDM = ERROR MESSAGE ROUTINE
C IER = ERROR FLAG USED BY QMINVZ
C QMINVZ = MATRIX INVERSION NUMBER CRUNCHER
C DET,IDET = DETERMINANT CONSTANTS
C PROCEDURE: THE ORDER OF X IS CHECKED TO MAKE SURE IT IS CORRECT AND
C QMINVZ IS CALLED. X AND Y MAY BE THE SAME MATRIX IN THE CALLING
C PROGRAM. N MUST BE 2 LESS THAN NCMX & NCMY & LESS THAN EQUAL
C TO NRMX & NRMY. IF X IS SINGULAR, A MESSAGE IS PRINTED AND
C EXECUTION CEASED.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),DET
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XINVRS')
NP2=N+2
CALL XDIMCH(NCMX,NP2,'INVRS',0)
CALL XDIMCH(NCMY,NP2,'INVRS',0)
CALL XDIMCH(NRMX,N,'INVRS',0)
CALL XDIMCH(NRMY,N,'INVRS',0)
IF(N.EQ.1)Y(1,1)=1.0D0/X(1,1)
IF(N.EQ.1)IZINCR = IZINCR - 1
IF(N.EQ.1)RETURN
I01=NCMY-1
IF(N.GE.I01)GOTO1
C CALL QMINVZ TO GET INVERSE
IER=1
CALL QMINVZ(NRMX,NRMY,NRMY,X,N,DET,IDET,Y,Y(1,I01),IER,X(1,NCMX))
IF(IER.NE.0)CALL QWRGDM(-9999,-9999)
IZINCR = IZINCR - 1
RETURN
C ERROR MESSAGE
1 CALL END (1)
WRITE(LUWN,2) N,NCMY
2 FORMAT('0* * * * * TO INVERT A SQUARE MATRIX ',
2 'THE ORDER OF THE SQUARE MATRIX (=',I4,') MUST BE 2 LESS
3 THAN THE DIMENSIONS (=',I4,').')
CALL END (2)
CALL END (3)
STOP
END
SUBROUTINE XKRONK(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NCX,NRY,NCY)
C TITLE: KRONECKER PRODUCT
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO FIND THE KRONECKER PRODUCT OF 2 MATRICES.
C SYMBOLS:
C X = 1ST INPUT MATRIX
C Y = 2ND INPUT MATRIX
C Z = RESULT MATRIX
C NRX,NCX = # OF ROWS & COLS OF X
C NRY,NCY = # OF ROWS & COLS OF Y
C NRZ,NCZ = # OF ROWS & COLS OF Z
C NRMX, NRMY & NRMZ = DIMENSIONS OF X,Y, AND Z IN CALLING PROGRAM
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C I,J,K,L,M,N,IK,JL = SCRATCH VARIABLES
C PROCEDURE: A SINGLE ELEMENT OF X IS SCALAR MULTIPLIED BY Y AND THE
C RESULT MATRIX IS PLACED IN A PARTITION OF Z CORRESPONDING TO THE
C RELATIVE POSITION OF THE SINGLE ELEMENT OF X. THIS PRODUCT IS
C CALLED A KRONECKER PRODUCT. X AND Y MAY BE THE SAME ADDRESS,
C BUT Z MUST BE DIFFERENT.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),Z(NRMZ,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XKRONK')
CALL XDIMCH(NRMX,NRX,'KRONK',0)
CALL XDIMCH(NRMY,NRY,'KRONK',0)
NRZ=NRX*NRY
CALL XDIMCH(NRMZ,NRZ,'KRONK',0)
C LOOP AROUND ROWS OF X
N=-NRY
DO 1I=1,NRX
N=N+NRY
M=-NCY
C LOOP AROUND COLS OF X
DO 1J=1,NCX
M=M+NCY
C LOOP AROUND ROWS OF Y
DO 1K=1,NRY
IK=N+K
C LOOP AROUND COLS OF Y
DO 1L=1,NCY
JL=M+L
C ACTUAL PRODUCT
1 Z(IK,JL)=X(I,J)*Y(K,L)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE QLRSM(NRMA,N,A,D,E,IERR)
C*TITLE: LATENT ROOTS OF A SYMMETRIC MATRIX
C*PROGRAMMER: CARL FINKBEINER
C*DATE: JANUARY,1975
C*LOCATION: UNIVERSITY OF ILLINOIS (SOUPAC)
C*PURPOSE: TO FIND ALL LATENT ROOTS (EIGENVALUES) OF A SYMMETRIC MATRIX
C* *BY THE IMPLICIT QL METHOD.
C*SYMBOLS:
C* *NRMA = ROW DIMENSION OF THE ARRAY A IN CALLING PROGRAM
C* *N = ORDER OF THE SYMMETRIC MATRIX STORED IN A
C* *A = 2 DIMENSIONAL ARRAY CONTAINING THE SQUARE SYMMETRIC MATRIX
C* * * *FOR WHICH THE LATENT ROOTS ARE TO BE FOUND.
C* *D = CONTAINS THE N LATENT ROOTS ON OUTPUT
C* *E = INTERMEDIATE STORAGE
C* *IERR = ERROR FLAG. IF IERR = 0, EXECUTION SUCCEEDED; OTHERWISE, IER
C* * * * * CONTAINS THE INDEX NUMBER OF THE EIGENVALUE ON WHICH EXECUTIO
C* * * * * FAILED.
C*PROCEDURE: A IS FIRST TRIDIAGONALIZED AND THE LATENT ROOTS OF THIS
C* *MATRIX ARE FOUND (THESE ROOTS ARE THE SAME AS THE ROOTS OF A).
C* *THE SUBROUTINES QTRED1 AND QIMTQL WHICH ACCOMPLISH THIS ARE EISPACK
C* *ROUTINE SEE THE LISTINGS OF THESE ROUTINES FOR A DESCRIPTION OF THE
C* *METHODS AND REFERENCES. THE A MATRIX IS PARTIALLY DESTROYED BY
C* *QTRED1, SO IT RESTORED IN THIS ROUTINE. THE ROOTS OUTPUT BY QIMTQL
C* *ARE IN ASCENDING ORDER AND SO ARE REORDERED PRIOR TO OUTPUT.
C* *THIS SUBROUTINE IS FASTER THAN EITHER LRVSM OR QLRSVM WHEN ONLY THE
C* *ROOTS OF A ARE DESIRED.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION A(NRMA,1),D(1),E(1)
COMMON/ZINCR/IZINCR
CALL QINCR ('QLRSM')
C*DIMENSIONALITY CHECK
IERR=-1
IF(NRMA.LT.N.OR.N.LE.0)IZINCR = IZINCR - 1
IF(NRMA.LT.N.OR.N.LE.0)RETURN
C*TRIDIAGONALIZE A
CALL QTRED1(NRMA,N,A,D,E,E)
C*COMPUTE ROOTS
CALL QIMTQL(N,D,E,IERR)
C*RECONSTRUCT INPUT MATRIX
DO 1I=1,N
J=I-1
DO 1K=1,J
1 A(I,K)=A(K,I)
C*ERROR RETURN
IF(IERR.NE.0)IZINCR = IZINCR - 1
IF(IERR.NE.0)RETURN
C*REORDER EIGENVALUES
M=N/2
J=N+1
DO 2I=1,M
J=J-1
TEM=D(I)
D(I)=D(J)
2 D(J)=TEM
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XLSQHY(NRMA,NRMS,NRMN,NRMW,NRMD,NRMC,A,S,
2 N,W,D,C,NR,NC,IFL)
C TITLE: LEAST-SQUARES HYPERPLANE FITTING
C PROGRAMMER: CARL FINKBEINER,ADAPTED FROM A PROGRAM BY L. TUCKER
C DATE: MAY, 1975
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO ROTATE A FACTOR MATRIX TO THE BEST (LEAST-SQUARES)
C HYPERPLANE FIT USING MARKERS SPECIFIED BY THE USER.
C SYMBOLS:
C A = INPUT FACTOR MATRIX
C S = SELECTION MATRIX (ALL 1'S AND 0'S)
C N = OUTPUT FACTOR NORMALS MATRIX - TRANSPOSED
C W,D,C = TEMPORARY WORK FILES
C NR = NUMBER OF ROWS OF A AND S
C NC = NUMBER OF COLS OF A AND S, AND DIMENSIONS OF B
C IFL = O TO SKIP INTERMEDIATE PRINTING,
C 1 TO PRINT INTERMEDIATE OUTPUT
C NRMA, NRMS, NRMN, NRMW, NRMD, NRMC = DIMENSIONS OF INPUT MATRICES
C IN CALLING PROGRAM,MUST BE .GE. 10
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C QLRSVM = IMPLICIT QL EIGEN ROUTINE
C PROCEDURE: ACCORDING TO LEDYARD TUCKER. NONE OF THE INPUT MATRICES
C MAY BE THE SAME IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION A(NRMA,1),N(NRMN,1),S(NRMS,1),
2 W(NRMW,1),D(NRMD,1),X,C(NRMC,1),Y
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XLSQHY')
CALL XDIMCH(NRMA,NR,'LSQHY',0)
CALL XDIMCH(NRMS,NR,'LSQHY',0)
CALL XDIMCH(NRMN,NR,'LSQHY',0)
CALL XDIMCH(NRMW,NR,'LSQHY',0)
CALL XDIMCH(NRMD,NR,'LSQHY',0)
CALL XDIMCH(NRMC,NR,'LSQHY',0)
I0 = NRMA
IF(NRMA.LT.10)GOTO12
I0 = NRMS
IF(NRMS.LT.10)GOTO12
I0 = NRMN
IF(NRMN.LT.10)GOTO12
I0 = NRMW
IF(NRMW.LT.10)GOTO12
I0 = NRMD
IF(NRMD.LT.10)GOTO12
I0 = NRMC
IF(NRMC.LT.10)GOTO12
C SOLUTION FOR EACH FITTED FACTOR
DO 1IFF=1,NC
C PRODUCT MATRIX FROM INPUT FACTORS FOR SELECTED VARIABLES
NVP=0
DO 2I=1,NC
DO 2J=1,I
2 C(I,J)=0.D0
DO 3I=1,NR
IF(S(I,IFF).EQ.0.D0)GOTO3
NVP=NVP+1
DO 4J=1,NC
DO 4K=1,J
4 C(J,K)=C(J,K)+A(I,J)*A(I,K)
3 CONTINUE
DO 5I=2,NC
DO 5K=1,I
5 C(K,I)=C(I,K)
C EIGENSOLUTION OF PRODUCT MATRIX
CALL QLRSVM(NRMC,NRMW,NRMD,NC,NC,C,W,D,D(1,2),D(1,3),D(1,4),
+ D(1,5),D(1,6),IERR)
IF(IERR.NE.0)GOTO13
IF(IFL.EQ.0)GOTO7
WRITE(LUWN,6)IFF,NVP
6 FORMAT('0',5X,'EIGENVALUES FOR ',
1'FITTING PLANE',I4,5X,'NUMBER OF ATTRIBUTES IN THE PLANE =',I4)
DO 8I=1,NC
8 WRITE(LUWN,9)I,D(I,1)
9 FORMAT(6X,I2,F12.5)
C MOVE LAST EIGENVECTOR TO N, REFLECT IF NECESSARY
7 X=0.D0
DO 10I=1,NC
Y=0.D0
DO 11J=1,NR
11 Y=Y+A(J,I)
10 X=X+Y*W(I,NC)
Y=-1.D0
IF(X.GE.0.D0)Y=1.D0
DO 1I=1,NC
1 N(I,IFF)=Y*W(I,NC)
IZINCR = IZINCR - 1
RETURN
12 CALL END (1)
WRITE(LUWN,14)I0
14 FORMAT ('0* * * * * YOUR CALL TO THE LSQHY PROGRAM HAS FAILED BECA
1USE THIS PROGRAM REQUIRES NO LESS THAN 10 COLUMNS FOR THE DIMENSIO
2NS OF'/12X,'THE MATRICES YOU HAVE SPECIFIED',I4,'.')
CALL END (2)
CALL END (3)
13 CALL END (1)
WRITE(LUWN,16)
16 FORMAT('0* * * * * SOMETHING WAS VERY WRONG WITH THE NUMBERS PASSE
1D TO THE LSQHY PROGRAM IN THE INPUT MATRIX.')
CALL END (2)
CALL END (3)
STOP
END
C SUBROUTINE QLRSVM(NRMA,NRMZ,NRMRV,N,NRV,A,Z,W,D,E,E2,IND,RV,IERR)
C*TITLE: LATENT ROOTS AND SOME VECTORS OF A SYMMETRIC MATRIX.
C*PROGRAMMER: CARL FINKBEINER
C*DATE: JANUARY, 1975
C*LOCATION: UNIVERSITY OF ILLINOIS (SOUPAC)
C*PURPOSE: FIND ALL LATENT ROOTS (EIGENVALUES) AND A SMALLER NUMBER OF
C* *LATENT VECTORS (EIGENVECTORS) OF A SYMMETRIC MATRIX.
C*SYMBOLS:
C* *NRMA, NRMZ, NRMRV = ROW DIMENSION OF A, Z, AND RV, RESPECTIVELY,
C* * * * * * * * * * * * IN CALLING PROGRAM.
C* *N = ORDER OF INPUT SYMMETRIC MATRIX (MUST BE .LE. NRMA OR NRMZ OR NRMRV)
C* *NRV = NUMBER OF LATENT VECTORS TO BE COMPUTED (CORRESPONDING TO
C* * * * *THE NRV LARGEST LATENT ROOTS)
C* *A = CONTAINS THE INPUT SYMMETRIC MATRIX IN ITS FIRST N ROWS
C* * * *AND N COLUMNS
C* *Z = CONTAINS THE OUTPUT LATENT VECTORS OF A IN ITS FIRST NRV COLUMNS
C* *W = ALL N LATENT ROOTS OF A
C* *D,E,E2 = TEMPORARY STORAGE ARRAYS OF LENGTH AT LEAST N. CONTAIN THE
C* * * * * * DIAGONAL, SUBDIAGONAL, AND SQUARED SUBDIAGONAL ELEMENTS
C* * * * * * OF THE TRIDIAGONAL FORM OF A.
C* *RV = TEMPORARY STORAGE ARRAY WITH AT LEAST N*5 ELEMENTS IN IT.
C* *IERR = ERROR FLAG. IF IERR = 0 ON OUTPUT, THE COMPUTATIONS SUCCEED
C* * * * * IF IERR = -(N+1), THE INPUT DIMENSIONS WERE WRONG. IF THE
C* * * * * LATENT ROOT COMPUTATION FAILED, IERR IS SET TO THE INDEX OF
C* * * * * THE LATENT ROOT AT WHICH IT FAILED. IF THE LATENT VECTOR
C* * * * * COMPUTATION FAILED, IERR IS SET TO -1 TIMES THE INDEX OF THE
C* * * * * VECTOR AT WHICH THE ERROR OCCURRED.
C* *IND = AN INTEGER*4 ARRAY OF AT LEAST N ELEMENTS USED BY THIS
C* * * * *SUBROUTINE TO REFER LATENT ROOTS TO SUBMATRICES OF THE
C* * * * *TRIDIAGONAL FORM.
C*PROCEDURE: THE MATRIX A IS TRIDIAGONALIZED AND THE ROOTS FOR EACH
C* *SUBMATRIX OF THIS TRIDIAGONAL FORM (DEFINED BY 0'S IN THE SUBDIAGO-
C* *NAL ARE FOUND BY THE IMPLICIT QL METHOD. THE ARRAY IND CONTAINS
C* *POINTERS REFERRING THE LATENT ROOTS TO THE SUBMATRIX FROM WHICH
C* *THEY WERE COMPUTED. THIS IS DONE BECAUSE THE ROOTS ARE SORTED
C* *PRIOR TO FINDING THE LATENT VECTORS AND THE SUBROUTINE WHICH
C* *DOES THIS (BY INVERSE ITERATION) NEEDS THIS INFORMATION. THE
C* *VECTORS ARE BACK TRANSFORMED SO THAT THEY ARE VECTORS OF A. THE
C* *MAIN NUMBER CRUNCHERS ARE EISPACK SUBROUTINES (QTRED1, QIMTQL,
C* *QTINVT,AND QTRBAK) AND CONTAIN FULL INFORMATION ON AND REFERENCES
C* *FOR THE ALGORITHMS USED. THIS SUBROUTINE COMPUTES ALL LATENT ROOTS
C* *BUT THE USER MAY REQUEST A SMALLER NUMBER OF THE LATENT VECTORS
C* *WHICH CORRESPOND TO THE NRV LARGEST LATENT ROOTS. THIS METHOD
C* *GUARANTEES THE LATENT VECTORS TO BE ORTHONORMAL (SEE QTINVT) EVEN
C* *FOR IDENTICAL OR NEAR IDENTICAL LATENT ROOTS. THIS SUBROUTINE IS
C* *FASTER THAN LRVSM FOR FINDING NRV (.LT. N) LATENT VECTORS OF A;
C* *HOWEVER, IF NRV .GT. N-8, THEN LRVSM USES LESS CORE THAN THIS
C* *SUBROUTINE AND SO IS TO BE PREFERRED IN THESE CIRCUMSTANCES. THE
C* *FASTEST WAY TO OBTAIN ALL N LATENT ROOTS IS WITH SUBROUTINE QLRSM.
SUBROUTINE QLRSVM(NRMA,NRMZ,NRMRV,N,NRV,A,Z,W,D,E,E2,IND,RV,IERR)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION A(NRMA,1),D(1),E(1),E2(1),W(1),RV(NRMRV,1)
DIMENSION IND(1),Z(NRMZ,1)
DOUBLE PRECISION MACHEP
COMMON/ZINCR/IZINCR
CALL QINCR ('QLRSVM')
C*MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING THE RELATIVE
C* *PRECISION OF FLOATING POINT ARITHMETIC.
MACHEP=.5D-14
IERR=-(N+1)
C*RESET NUMBER OF VECTORS REQUESTED IF INCORRECTLY
C*SPECIFIED AND IF IT IS POSSIBLE TO DEDUCE A CORRECTION
IF(N.LT.NRV)NRV=N
IF(NRV.LT.0)NRV=0
C*DIMENSIONALITY CHECK
CALL XDIMCH(NRMA,N,'QLRSVM',0)
CALL XDIMCH(NRMZ,N,'QLRSVM',0)
C*TRIDIAGONALIZE A
CALL QTRED1(NRMA,N,A,D,E,E2)
C*TEST FOR SMALL SUBDIAGONAL ELEMENTS AND FIND LATENT ROOTS FOR
C*SUBMATRICES AND BUILD IND TO HAVE SUBMATRIX NUMBERS.
INDI=1
IND(1)=1
E2(1)=2.D0
W(1)=D(1)
L=1
DO 20IJ=2,N
IF(DABS(E(IJ)).LE.MACHEP*(DABS(D(IJ))+DABS(D(IJ-1))))GOTO30
W(IJ)=D(IJ)
RV(IJ,2)=E(IJ)
IND(IJ)=INDI
20 CONTINUE
IJ=N+1
30 LK=IJ-L
C*FIND LATENT ROOTS OF SUBMATRIX
CALL QIMTQL(LK,W(L),RV(L,2),IERR)
IF(IERR.NE.0)GOTO10
IF(IJ.EQ.N+1)GOTO40
C*SET NEW IND VALUE AND SET E2 ELEMENT TO TRUE 0.
INDI=INDI+1
IND(IJ)=INDI
E2(IJ)=0.D0
W(IJ)=D(IJ)
L=IJ
GOTO 20
C*SORT LATENT ROOTS AND REORDER IND CORRESPONDINGLY (BUBBLE SORT)
40 NM1=N-1
I=1
1 K=N-I
J=1
2 J1=J+1
IF(W(J).GE.W(J1))GOTO3
TEM=W(J)
W(J)=W(J1)
W(J1)=TEM
L=IND(J)
IND(J)=IND(J1)
IND(J1)=L
3 J=J1
IF(J.LE.K)GOTO2
I=I+1
IF(I.LE.NM1)GOTO1
C*OBTAIN LATENT VECTORS (IF REQUESTED) OF TRIDIAGONAL MATRIX
C* *BY INVERSE ITERATION.
70 IF(NRV.EQ.0)GOTO10
CALL QTINVT(NRMZ,N,D,E,E2,NRV,W,IND,Z,IERR,RV,
2 RV(1,2),RV(1,3),RV(1,4),RV(1,5))
IF(IERR.NE.0)GOTO10
C*BACK TRANSFORM THE LATENT VECTORS SO THEY ARE LATENT VECTORS OF A.
CALL QTRBAK(NRMA,NRMZ,N,A,E,NRV,Z)
C*RESTORE A
10 DO 150I=1,N
J=I-1
DO 150K=1,J
150 A(I,K)=A(K,I)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XMATML(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NXY,NCY)
C TITLE: MATRIX MULTIPLY
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO MULTIPLY 2 INPUT MATRICES AND RETURN THE RESULT IN A THIRD
C MATRIX
C SYMBOLS:
C X,Y = THE 2 INPUT MATRICES
C Z = THE RESULT MATRIX = X*Y
C NRX = # OF ROWS OF X = # OF ROWS OF Z
C NXY = # OF COLUMNS OF X = # OF ROWS OF Y
C NCY = # OF COLUMNS OF Y = # OF COLUMNS OF Z
C NRMX, NRMY, & NRMZ = DIMENSIONALITY OF X, Y, AND Z IN CALLING PROGRAM
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C I,J,K = SCRATCH VARIABLES
C PROCEDURE: EACH ROW OF X IS VECTOR MULTIPLIED BY EACH COLUMN OF Y TO
C PRODUCE THE SUCCESSIVE ELEMENTS OF Z. NOTE THAT X AND Y MAY BE THE
C MATRIX IN THE CALLING PROGRAM, HOWEVER, Z MUST BE A DIFFERENT MATRIX
C Z WILL HAVE NRX ROWS AND NCY COLUMNS.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DIMENSION X(NRMX,1),Y(NRMY,1),Z(NRMZ,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XMATML')
CALL XDIMCH(NRMX,NRX,'MATML',0)
CALL XDIMCH(NRMZ,NRX,'MATML',0)
CALL XDIMCH(NRMY,NXY,'MATML',0)
C MULTIPLY EACH ROW OF X BY EACH COLUMN OF Y AND STORE IN Z
DO 4I=1,NRX
DO 4J=1,NCY
Z(I,J)=0.0D0
DO 4K=1,NXY
4 Z(I,J)=Z(I,J)+X(I,K)*Y(K,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XPARTT(NRMX,NRMY,X,Y,IBR,IER,IBC,IEC)
C TITLE: PARTITION
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO SELECT OFF A PARTITION OF A MATRIX.
C SYMBOLS:
C X = INPUT MATRIX TO BE PARTITIONED.
C Y = PARTITION SELECTED OUT OF X.
C IBR,IER = BEGINNING & ENDING ROWS, RESPECTIVELY, OF PARTITION OF X.
C IBC,IEC = BEGINNING & ENDING COLS, RESPECTIVELY, OF PARTITION OF X.
C A = LABELLED COMMON AREA CONTAINING IPGM.
C IPGM = PRGRAM COUNTER
C NRMX & NRMY = DIMENSIONS OF X AND Y IN CALLING PROGRAM.
C QINCR = PROGRAM COUNTER SUBROUTINE.
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE.
C I,J,K,L = SCRATCH VARIABLES.
C PROCEDURE: A PARTITION OF X IS DESIGNATED BY IBR,IER,IBC,IEC. THE
C SUBMATRIX DESIGNATED BY THIS PARTITION IS COPIED INTO THE UPPER LEFT
C OF Y. A TEST IS MADE TO ASCERTAIN THAT IER(IEC).GE.IBR(IBC) AND THAT
C IER(IEC) IS NOT GREATER THAN NRMX. IF THIS TEST IS NOT SATISFIED, A
C DIAGNOSTIC ERROR MESSAGE IS PRINTED AND EXECUTION IS TERMINATED. X
C MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XPARTT')
C CHECK IBR,IER,IBC,IEC FOR CONFORMABILITY AMONG THEMSELVES AND COMPARED
C IF SOMETHING IS WRONG, WRITE ERROR MESSAGE AND TERMINATE.
IF(IBR.GT.IER.OR.IER.GT.NRMX.OR.IBC.GT.IEC.OR.IER.GT.NRMY)GOTO1
C EXTRACT A SUBMATRIX FROM X AND STORE IN Y.
K=0
DO 2I=IBR,IER
K=K+1
L=0
DO 2J=IBC,IEC
L=L+1
2 Y(K,L)=X(I,J)
IZINCR = IZINCR - 1
RETURN
C ERROR MESSAGE
1 CALL END (1)
WRITE(LUWN,3) IBR,IER,IBC,IEC
3 FORMAT('0* * * * * THE VALUES SPECIFYING THE PARTITION TO BE CREAT
1ED IN SUBROUTINE PARTT ARE INCORRECT.'/
412X,'THESE VALUES ARE: BEGINNING ROW # =',I4/33X,'ENDING ROW # =',
5I4/30X,'BEGINNING COL # =',I4/33X,'ENDING COL # =',I4)
CALL END (2)
CALL END (3)
STOP
END
SUBROUTINE XPRDDI(NRMX,X,N,PD)
C TITLE: PRODUCT OF THE DIAGONAL
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO FIND THE PRODUCT OF THE ELEMENTS IN THE MAJOR DIAGONAL OF
C INPUT MATRIX
C SYMBOLS:
C X = INPUT MATRIX
C N = # OF ELEMENTS IN MAJOR DIAGONAL
C PD = THE RESULTING PRODUCT
C NRMX = DIMENSIONS OF X IN CALLING PROGRAM
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C I = SCRATCH VARIABLE
C PROCEDURE:
C THE ELEMENTS OF THE MAJOR DIAGONAL OF X ARE SUCCESSIVELY MULTIPLIED
C BY THE PRODUCT OF THE PRECEDING ELEMENTS. X MAY BE NON-SQUARE IN
C WHERE THE N, THE LENGTH OF THE MAJOR DIAGONAL, IS EQUAL TO MIN(# OF
C ROWS OF X, MAY BE THE SAME AS COLS OF X).
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DOUBLE PRECISION X(NRMX,1),PD
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XPRDDI')
CALL XDIMCH(NRMX,N,'PRDDI',0)
C MULTIPLY DIAGONAL ELEMENTS AND STORE IN PD.
PD=1.D0
DO 1I=1,N
1 PD=PD*X(I,I)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XPRINT(NRMX,X,NR,NC,IFLAG)
C TITLE: MATRIX PRINT
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO PRINT A MATRIX IN F OR D FORMAT
C SYMBOLS:
C X = MATRIX TO BE PRINTED
C NR = # OF ROWS OF X
C NC = # OF COLS OF X
C IFLAG = FLAG SPECIFYING FORMAT TYPE (EITHER 'F' OR 'D' ON INPUT TO
C PRINT
C FMT = PRINTING FORMAT
C IF = A CONSTANT USED IN TESTING IFLAG
C NRMX = DIMENSIONALITY OF X
C A,C = ALTERNATIVE SECOND HALVES OF PRINT FORMAT
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C L = # OF 9 COLUMN BLOCKS OF X, TO BE PRINTED ONE AT A TIME
C N,I = BEGINNING AND ENDING COLUMN NUMBERS FOR A BLOCK
C JK,K = COLUMN AND ROW HEADINGS
C J,M,IF1,IF2,IFLAG1,IFLAG2 = SCRATCH VARIABLES
C PROCEDURE:
C OUTPUT FORMAT TYPE IS DECIDED ON BY TESTING IFLAG AGAINST IF
C (IF IS INITIALIZED TO AN INTEGER CORRESPONDING TO THE CHARACTER 'F')
C OBTAINED; N AND I ARE OBTAINED FOR EACH BLOCK AND USED TO SELECT OUT
C PARTITIONS OF X FOR PRINTING IN A BLOCK. ROW AND COLUMN HEADINGS
C ARE PRINTED; A DOUBLE SPACE OCCURS BEFORE RETURNING.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),FMT(2),A,C
C INTEGER*2 IFLAG,IF/'F '/,IFLAG2/' '/,IFLAG1
INTEGER IFLAG,IF,IFLAG2,IFLAG1
C LOGICAL*1 IF1,IF2
LOGICAL IF1,IF2
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
EQUIVALENCE(IF1,IFLAG1),(IF2,IFLAG2)
DATA IF/'F '/,IFLAG2/' '/
DATA FMT/'(''0'',I6,',' '/,A/'9F14.5)'/,C/'9E14.5)'/
C BOOKKEEPING
CALL QINCR ('XPRINT')
CALL XDIMCH(NRMX,NR,'PRINT',0)
C CHOOSE 2ND HALF OF FMT
FMT(2)=C
IFLAG1=IFLAG
IF2=IF1
IF(IFLAG2.EQ.IF)FMT(2)=A
C FIND # OF BLOCKS TO BE PRINTED
L=NC/9
IF(L*9.NE.NC)L=L+1
I=0
C DO LOOP FOR PRINTING BLOCKS
DO 3J=1,L
C FIND BEGINNING & ENDING COLUMN #'S FOR CURRENT BLOCK
N=I+1
I=MIN0(NC,J*9)
C PRINT COLUMN HEADINGS FOR BLOCK
WRITE(LUWN,4)(JK,JK=N,I)
4 FORMAT('0',2X,9I14/)
C PRINT ALL ROWS OF THIS BLOCK WITH ROW HEADINGS
DO 3K=1,NR
3 WRITE(LUWN,FMT)K,(X(K,M),M=N,I)
C DOUBLE SPACE AFTER MATRIX HAS BEEN PRINTED
WRITE(LUWN,5)
5 FORMAT('0')
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XPUNCH(NRMX,X,NR,NC)
C TITLE: MATRIX PUNCH
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE,1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO PUNCH A MATRIX ONTO CARDS
C SYMBOLS:
C X = MATRIX TO BE PUNCHED
C NR = # OF ROWS OF X TO BE PUNCHED
C NC = # OF COLS OF X TO BE PUNCHED
C NRMX = DIMENSIONALITY OF X
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C L = # OF CARDS PER ROW OF X
C M,J = BEGINNING AND ENDING ELEMENTS OF THE PART OF THE ROW OF X BEIN
C PUNCHED ON 1 CARD
C I,K = ROW # AND CARD #
C N = SCRATCH VARIABLE
C PROCEDURE: A HEADER CARD OF THE FORM RECOGNIZED BY THE INPUT SUBROUTIN
C FORMAL IS PUNCHED FIRST. THE NUMBER OF CARDS PER ROW IS CALCULATED
C IN THE DOUBLY NESTED DO LOOP WHEREIN EACH ROW OF X IS PUNCHED ONTO
C SUCCESSIVE CARDS WITH SEQUENCE NUMBERS.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XPUNCH')
CALL XDIMCH(NRMX,NR,'PUNCH',0)
C PUNCH FORMAL HEADER CARD
WRITE (LUWN, 5)
5 FORMAT (' SORRY BUT LEDERLY DOES NOT HAVE A CARD PUNCH.'/
2 ' PROCESSING SHALL TERMINATE')
STOP
9996 WRITE(LUPN,1)NR,NC
1 FORMAT(2I4,'(10X,5D14.7)')
C FIND # OF CARDS PER ROW OF X
L=NC/5
IF(L*5.NE.NC)L=L+1
C DO LOOP FOR ROWS
DO 2I=1,NR
J=0
C DO LOOP FOR CARDS/ROW
DO 2K = 1,L
C FIND BEGINNING AND ENDING COLUMNS OF PART OF ROW TO BE PUNCHED
M=J+1
J=MIN0(NC,K*5)
C PUNCH A CARD WITH SEQUENCE #'S
2 WRITE(LUPN,3)I,K,(X(I,N),N=M,J)
3 FORMAT(2I5,5D14.7)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XRECIP(NRMX,NRMY,X,Y,NR,NC,IZERO)
C TITLE: RECIPROCAL
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE,1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO FIND THE RECIPROCAL OF EVERY ELEMENT OF A MATRIX.
C SYMBOLS:
C X = INPUT MATRIX
C Y = OUTPUT MATRIX CONTAINING THE RECIPROCALS OF X
C NR,NC = # OF ROWS AND COLS OF X AND Y
C IZERO = 0 TO LEAVE 0 ELEMENTS AS 0'S;
C 1 TO TERMINATE PROCESSING IF A 0 ELEMENT IS FOUND
C DX = SCRATCH VARIABLE
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C DABS = FORTRAN SUPPLIED DOUBLE PRECISION DABSOLUTE VALUE FUNCTION
C I,J = SCRATCH VARIABLES
C PROCEDURE:
C EVERY ELEMENT OF X IS TESTED AGAINST 10**-70 AS THE CRITERIA OF
C FUNCTIONAL ZERO. IF THE ELEMENT IS NOT FUNCTIONALY ZERO, THE
C RECIPROCAL OF THE ELEMENT IS STORED IN THE CORRESPONDING
C ELEMENT OF Y. IF THE ELEMENT IS FUNCTIONALY ZERO,IT IS SET
C TO 0 OR PROCESSING IS TERMINATED, DEPENDINMG ON THE INPUT VALUE
C OF IZERO. X AND Y MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),DX
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XRECIP')
CALL XDIMCH(NRMX,NR,'RECIP',0)
CALL XDIMCH(NRMY,NR,'RECIP',0)
C INITIATE LOOP THRU X AND Y
DO 1I=1,NR
DO 1J=1,NC
C TEST ELEMENT OF X AND SEE IF IT IS A FUNCTIONAL ZERO.
DX=DABS(X(I,J))
IF(DX.GT.1.D-30)GOTO4
IF(IZERO.EQ.1)GOTO2
Y(I,J)=0.D0
GOTO 1
C STORE RECIPROCAL OF X IN Y
4 Y(I,J)=1.D0/X(I,J)
1 CONTINUE
IZINCR = IZINCR - 1
RETURN
C ERROR MESSAGE.
2 CALL END (1)
WRITE(LUWN,3)I,J
3 FORMAT('0* * * * * THE ERROR OPTION FOR THE RECIPROCAL SUBROUTINE
1IS ON.'/11X,
2 'ELEMENT',I4,',',I4,' OF THE INPUT MATRIX IS ZERO')
CALL END (2)
CALL END (3)
STOP
END
SUBROUTINE XREORD(NRMX,X,NR,NC)
C TITLE: REORDER
C PROGRAMMER: T. WANG
C DATE: SEPTEMBER, 1971
C LOCATION: UNIVERSITY OF ILLINOIS
C DESCRIPTION: THIS SUBROUTINE WAS ADAPTED SLIGHTLY FROM THE CODE PROVID
C FORTUOI TO ACCOMPANY QGELGZ FOR PURPOSES OF REARRANGING G SO IT IS
C EXPECTED MATRIX FORM (I.E. SUCH THAT AG=B). THIS SUBROUTINE IS CALL
C BY INVRS AFTER QGELGZ AND IS PASSED G ALONG WITH ITS DIMENSIONALITY
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1)
COMMON/ZINCR/IZINCR
CALL QINCR ('XREORD')
K=1
L=NR*NC
10 IF(L.LT.NRMX)GOTO20
L=L-NRMX
K=K+1
GOTO 10
20 DO 40JJ=1,NC
J=NC+1-JJ
DO 40II=1,NR
I=NR+1-II
X(I,J)=X(L,K)
IF(L.EQ.1)GOTO30
L=L-1
GOTO 40
30 IF(K.EQ.1)GOTO40
K=K-1
L=NRMX
40 CONTINUE
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XRO2DI(NRMX,NRMY,X,Y,N)
C TITLE: ROW-TO-DIAGONAL MATRIX
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE:TO RAISE THE FIRST ROW OF AN INPUT MATRIX TO A DIAGONAL MATRIX
C SYMBOLS:
C X = INPUT MATRIX CONTAINING THE ROW IN ITS FIRST ROW.
C Y = RESULTING DIAGONAL MATRIX.
C N = # OF COLS OF X
C NRMX & NRMY = DIMENSIONALITY OF X AND Y IN CALLING PROGRAM.
C QINCR = PROGRAM COUNTER SUBROUTINE.
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE.
C I,J,K = SCRATCH VARIABLES.
C PROCEDURE: THE ELEMENTS OF THE FIRST ROW OF X ARE COPIED INTO THE MAIN
C DIAGONAL OF Y AND THE REMAINING ELEMENTS OF Y ARE ZEROED OUT. IF X IS
C 1 X 1, THE COPY IS MADE AND NO ZEROING OCCURS. NOTE THAT X MUST BE
C A 1 X 1 MATRIX, A SCALAR, OR A 2 DIMENSIONAL ARRAY IN THE CALLING
C PROGRAM X CANNOT BE A 1 DIMENSIONAL ARRAY. X AND Y MAY BE THE SAME
C MATRIX IN CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XRO2DI')
CALL XDIMCH(NRMY,N,'RO2DI',0)
CALL XDIMCH(NRMX,N,'RO2DI',0)
C SET UPPER LEFT ELEMENT
Y(1,1)=X(1,1)
C IF X IS A SCALAR, RETURN.
IF(N.EQ.1)IZINCR = IZINCR - 1
IF(N.EQ.1)RETURN
C LOOP FOR MOVING THE FIRST ROW OF X INTO THE DIAGONAL OF Y AND ZEROING
C THE REST OF Y.
DO 1I=2,N
Y(I,I)=X(1,I)
J=I-1
DO 1K=1,J
Y(I,K)=0.D0
1 Y(K,I)=0.D0
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XRO2MT(NRMX,NRMY,X,Y,NR,NC)
C TITLE: ROW TO MATRIX
C PROGRAMMER: CARL FINKBEINER
C DATE: AUGUST,1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: REPEAT A ROW TO CREATE A AMTRIX.
C SYMBOLS:
C X = INPUT MATRIX TO BE REPEATED
C Y = OUTPUT MATRIX WITH EACH ROW EQUAL TO X
C NR = # OF ROWS TO BE GENERATED IN Y
C NC = # OF COLS OF X & Y
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C PROCEDURE: THE FIRST ROW OF X IS COPIED INTO Y UNTIL Y HAS NR ROWS.
C X & Y MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XRO2MT')
CALL XDIMCH(NRMY,NR,'RO2MT',0)
C COPY THE FIRST ROW OF X INTO Y UNTIL Y HAS NR ROWS
DO 1I=1,NR
DO 1J=1,NC
1 Y(I,J)=X(1,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XSCAML(NRMX,NRMY,C,X,Y,NR,NC)
C TITLE: SCALAR MULTIPLICATION
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: TO MULTIPLY EVERY ELEMENT OF A MATRIX BY A SCALAR
C SYMBOLS:
C C = A REAL TYPE SCALAR. NOTE - C MUST BE A REAL NUMBER, NOT AN
C INTEGER
C X = INPUT MATRIX.
C Y = OUTPUT MATRIX CONTAINING C*X.
C NR = # OF ROWS OF X
C NC = # OF COLS OF X
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C PROCEDURE:
C EVERY ELEMENT OF X IS MULTIPLIED BY C AND THE RESULT IS STORED
C IN Y X AND Y MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
C C NUST BE A REAR NUMBER; IF C HAS AN INTEGER VALUE, A DECIMAL
C POINT MUST BE PRESENT IN ORDER TO INDICATE THAT C IS A REAL TYPE
C NUMBER. EXAMPLE: USE '1.' RATHER THAN JUST '1'. ALSO NOTE
C THAT C MAY BE A MATRIX, IN WHICH CASE THE FINAL ELEMENT OF
C THE MATRIX IS USED AS C.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XSCAML')
CALL XDIMCH(NRMX,NR,'SCAML',0)
CALL XDIMCH(NRMY,NR,'SCAML',0)
C MULTIPLY EVERY ELEMENT OF X BY THE REAL # C AND STORE IN Y
DO 1I=1,NR
DO 1J=1,NC
1 Y(I,J)=C*X(I,J)
IZINCR = IZINCR - 1
RETURN
END
SUBROUTINE XSQRT(NRMX,NRMY,X,Y,NR,NC,EPS)
C TITLE: SQUARE ROOT
C PROGRAMMER: CARL FINKBEINER
C DATE: JUNE, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: FIND THE SQUARE ROOT OF ALL ELEMENTS OF A MATRIX
C SYMBOLS:
C X = INPUT MATRIX
C Y = OUTPUT MATRIX CONTAINING SQUARE ROOTS OF X
C NR,NC = # OF ROWS AND COLS OF X AND Y
C EPS = LOWER BOUND TOLERANCE FOR 0. VALUES (A NEGATIVE NUMBER)
C A = LABELLED COMMON AREA CONTAINING IPGM
C IPGM = CURRENT PROGRAM #
C QINCR = PROGRAM COUNTER SUBROUTINE
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE
C I,J = SCRATCH VARIABLES
C DSQRT = FORTRAN SUPPLIED DOUBLE PRECISION SQUARE ROOT FUNCTION
C PROCEDURE:
C IF ANY ELEMENT OF X IS NEGATIVE, AN APPROPRIATE ERROR MESSAGE
C WILL WRITTEN OUT AND EXECUTION IS TERMINATED. OTHERWISE, THE
C SQUARE ROOT OF EVERY ELEMENT OF X IS STORED IN THE CORRESPON-
C DING ELEMENT OF Y. X AND Y MAY BE THE SAME MATRIX IN THE
C CALLING PROGRAM. ANY ELEMENT BETWEEN 0.0 AND EPS IS SET TO 0.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),X1
COMMON/A/IPGM,ISUBRU
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XSQRT')
CALL XDIMCH(NRMX,NR,'SQRT',0)
CALL XDIMCH(NRMY,NR,'SQRT',0)
EPS1=EPS
IF(EPS1.GT.0.)EPS1=0.D0
C INITIATE LOOP THRU X AND Y
DO 1I=1,NR
DO 1J=1,NC
C IF AN ELEMENT OF X IS NEGATIVE, WRITE ERROR MESSAGE AND TERMINATE
X1=X(I,J)
IF(X1.LT.0..AND.X1.GE.EPS1)X1=0.D0
IF(X1.LT.0.)GOTO2
C STORE SQUARE ROOT OF X IN Y
1 Y(I,J)=DSQRT(X1)
IZINCR = IZINCR - 1
RETURN
C ERROR MESSAGE
2 CALL END (1)
WRITE(LUWN,3)
3 FORMAT('0* * * * * O O P S !!! YOU CAN''T TAKE THE SQUARE ROOT O
1F A NEGATIVE NUMBER.'//25X,'CHECK THE ELEMENTS OF THE MATRIX WHICH
2 YOU PASSED TO THE SQRT PROGRAM')
CALL END (2)
CALL END (3)
STOP
END
SUBROUTINE XSUB(NRMX,NRMY,NRMZ,X,Y,Z,NR,NC)
C TITLE: SUBTRACTION OF 2 MATRICES
C PROGRAMMER: CARL FINKBEINER
C DATE: SEPTEMBER, 1974
C LOCATION: UNIVERSITY OF ILLINOIS
C PURPOSE: SUBTRACTION OF 2 MATRICES
C SYMBOLS:
C X = 1ST INPUT MATRIX.
C Y = 2ND INPUT MATRIX.
C Z = RESULTANT MATRIX.
C NR,NC = # OF ROWS AND COLUMNS, RESPECTIVELY, OF X, Y, AND Z.
C QINCR = PROGRAM COUNTER SUBROUTINE.
C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE.
C I,J = SCRATCH VARIABLES.
C PROCEDURE: CORRESPONDING ELEMENTS OF X AND Y ARE SUBTRACTED
C AND THE RESULT IS STORED IN THE CORRESPONDING ELEMENTS OF Z. ANY OF
C X, Y, OR Z MAY BE THE SAME MATRIX IN THE CALLING PROGRAM.
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
DOUBLE PRECISION X(NRMX,1),Y(NRMY,1),Z(NRMZ,1)
COMMON/ZINCR/IZINCR
C BOOKKEEPING
CALL QINCR ('XSUB')
CALL XDIMCH(NRMX,NR,'SUB',0)
CALL XDIMCH(NRMY,NR,'SUB',0)
CALL XDIMCH(NRMZ,NR,'SUB',0)
C SUBTRACT X-Y AND STORE THE RESULT IN Z
DO 1I=1,NR
DO 1J=1,NC
1 Z(I,J)=X(I,J)-Y(I,J)
IZINCR = IZINCR - 1
RETURN
END
FUNCTION CHIPRA (X, Y)
IMPLICIT DOUBLE PRECISION (A-H, O-Z)
COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG
C
C TO COMPUTE THE PROBABILITY (CHIPRA) GIVEN A CHI SQUARE
C VARIATE (X) WITH ITS ASSOCIATED DEGREES OF FREEDOM (Y).
C THIS APPROX. IS NOT GOOD FOR D. F. LESS THAN 30.
C THIS APPROXIMATION IS IN ABRAMOWITZ AND STEGUN, HANDBOOK
C OF MATHEMATICAL FUNCTIONS.
C
Z = ((X/Y)**(.3333333333333333D0) - (1.0D0 - (1.0D0/(4.5D0*Y))))
2 /DSQRT(1.0D0/(4.5D0*Y))
WRITE (LUWN, 1)
CHIPRA = 2.0D0
C IF (Y.LT.30.0D0) WRITE (LUWN, 1)
1 FORMAT (' * * * * THERE IS AN ERROR IN THE CHI SQUARE ',
2'PROBABILITY * * * *')
RETURN
END