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