SUBROUTINE CPLOT(A,NROW,NCOL,XMIN,XMAX,YMIN,YMAX,ID) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(I1,I1) COMMON/B/I1 COMMON/ZINCR/IZINCR CALL QINCR ('XCPLOT') NRMA = I1 NCMA = I1 CALL XCPLOT(NRMA,NCMA,A,NROW,NCOL,XMIN,XMAX,YMIN,YMAX,ID) IZINCR=IZINCR-1 RETURN END SUBROUTINE PLOT(A,NROW,NCOL,ID) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(I1,I1) COMMON/B/I1 COMMON/ZINCR/IZINCR CALL QINCR ('PLOT') NRMA = I1 NCMA = I1 CALL XPLOT(NRMA,NCMA,A,NROW,NCOL,ID) IZINCR=IZINCR-1 RETURN END SUBROUTINE XCPLOT(NRMA,NCMA,A,NROW,NCOL,XMIN,XMAX,YMIN,YMAX,IDD) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(NRMA,1) COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XIN(5),OUTPUT(5),FRSTPG COMMON/ZINCR/IZINCR DATA IONE,ITWO/1,2/ CALL QINCR ('XCPLOT') GOTO 673 ENTRY XPLOT (NRMA,NCMA,A,NROW,NCOL,IDD) CALL QINCR ('XPLOT') XMIN = 0.0D0 YMIN = XMIN XMAX = XMIN YMAX = XMIN 673 ID = IDD IF (IDD.LT.-2.OR.IDD.GT.2) ID = -1 IF (IDD.LT.-2.OR.IDD.GT.2) CALL END (1) IF (IDD.LT.-2.OR.IDD.GT.2) WRITE (LUWN, 65) IDD 65 FORMAT ( '0',32X,'YOU TOLD THE PLOT', 3 ' PROGRAM TO USE AN INVALID TYPE OF PLOT STYLE'/'0',32X, 4 'YOU SAID USE IFLAG = ',I10,' IFLAG CAN ONLY TAKE VALUES FROM', 5 /'0',32X,'-2 TO +2. PROCESSING SHALL CONTINUE WITH ID SET', 6 ' TO -1!!') CALL XDIMCH (NRMA, IONE, 'PLOT', 0) CALL XDIMCH (NCMA, IONE, 'PLOT', 0) CALL XDIMCH (NCOL, IONE, 'PLOT', 0) NCOLP1 = NCOL + 1 C IF (IDD.NE.0) GOTO 60 CALL XDIMCH (NCMA, NCOLP1, 'PLOT', 0) ID = -2 GOTO 35 C 60 IF (NCOL.GT.1) GO TO 500 35 CALL XDIMCH (NCMA, ITWO, 'PLOT', 0) DO 70 I=1,NROW 70 A(I, NCMA) = I C DO 2 K = 1, NCOL WRITE (LUWN,75) K,NROW 75 FORMAT ('1 SUCCESSIVE VALUE PLOT. VARIABLE NUMBER:', I4, 2' ON Y AXIS',40X,I3,9H SUBJECTS) C 2 CALL ZPLOT (A(IONE,NCMA),A(IONE,K),XMIN,XMAX,YMIN,YMAX,NROW,ID) GO TO 510 C 500 NPLM = NCOL-1 DO 506 I=1,NPLM IPL = I+1 DO 505 J=IPL,NCOL C WRITE (LUWN,3003) J,I,NROW 3003 FORMAT(1H1,5X,16H PLOT VARIABLE,I3,22H (X-AXIS) VS. VARIABLE,I3, 19H (Y-AXIS),49X,I3,9H SUBJECTS) C CALL ZPLOT (A(IONE,J),A(IONE,I),XMIN,XMAX,YMIN,YMAX,NROW,ID) C 505 CONTINUE 506 CONTINUE 510 IZINCR=IZINCR - 1 RETURN END SUBROUTINE ZPLOT (X,Y,ZXMIN,ZXMAX,ZYMIN,ZYMAX,NPOIZ,ID) C C ROUTINE TO GENERATE A ONE PAGE PLOT OF ARRAY-X- VS -Y- C C THE PARAMETERS -XMAX- -XMIN- -YMAX- -YMIN- INDICATE THE UPPER C AND LOWER BOUNDS FOR EACH AXIS. IF XMAX=XMIN THE ROUTINE GENERATES C ITS OWN BOUNDS FOR THE X AXIS, AND SIMILARLY IF YMAX=YMIN. C C IT IS ASSUMED THAT THE ARRAYS HAVE -NPOI- ENTRIES. C C IF -ID- IS POSITIVE, AXES WILL BE INCLUDED ON THE PLOT, C IF -ID- IS NEGATIVE, NO AXES WILL APPEAR. C IF ID IS PLUS OR MINUS 2, THE POINTS WILL BE IDENTIFIED, C OTHERWISE THEY WILL BE COUNTED. C C WRITTEN BY FORREST W. YOUNG, NOVEMBER, 1965 C VERSION 3, APRIL 1967 C MINOR UPDATE, JANUARY, 1968 C THIS PROGRAM WAS TAKEN FROM SOUPAC AT THE UNIVERSITY OF ILLINOIS C IT WAS MODIFIED FOR LEDERLE LABS C BY C DR. ALLEN FLEISHMAN C OCTOBER, 1979 C C IMPLICIT DOUBLE PRECISION (A-G,O-Z) REAL ITEM,HOLL,PTID,AIE,AMINUS,DD,PLUS,BLANK DOUBLE PRECISION X(1), Y(1), SMALL(21) DIMENSION ITEM(53,101),HOLL(11),PTID(124) COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XIN(5),OUTPUT(5),FRSTPG COMMON/ZINCR/IZINCR CALL QINCR ('ZPLOT') C DATA HOLL/1H ,1HX,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9,1HM/, XAIE,AMINUS,DD,PLUS,BLANK/1HI,1H-,1H.,1H*,1H /, XPTID/1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO, X1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H0,1H1,1H2,1H3,1H4, X1H5,1H6,1H7,1H8,1H9,1Ha,1Hb,1Hc,1Hd,1He,1Hf,1Hg,1Hh,1Hi,1Hj,1Hk, X1Hl,1Hm,1Hn,1Ho,1Hp,1Hq,1Hr,1Hs,1Ht,1Hu,1Hv,1Hw,1Hx,1Hy,1Hz, X1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM,1HN,1HO, X1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ,1H0,1H1,1H2,1H3,1H4, X1H5,1H6,1H7,1H8,1H9,1Ha,1Hb,1Hc,1Hd,1He,1Hf,1Hg,1Hh,1Hi,1Hj,1Hk, X1Hl,1Hm,1Hn,1Ho,1Hp,1Hq,1Hr,1Hs,1Ht,1Hu,1Hv,1Hw,1Hx,1Hy,1Hz/ 3001 FORMAT(14X,103H*.****.****.****.****.****.****.****.****.****.**** 1.****.****.****.****.****.****.****.****.****.****.*) 3300 FORMAT(1H ,F11.2,1X,105A1,F11.2) 3301 FORMAT(15X,A1,10(F9.4,A1)) 3302 FORMAT(10X,11F10.4) XMAX = ZXMAX XMIN = ZXMIN YMAX = ZYMAX YMIN = ZYMIN NPOI = NPOIZ IF (ID*ID.NE.4.AND.NPOIZ.GT.124) NPOI = 124 DO 115 I=1,53 DO 115 J=1,101 115 ITEM(I,J)=HOLL(1) IF(ID.LE.0)GO TO 119 DO 117 I=1,53 117 ITEM(I,51)=AIE DO 118 I=1,101 118 ITEM(27,I)=AMINUS 119 CONTINUE IF((XMAX-XMIN).GT.1.0E-4)GO TO 122 XMAX=-1.0E30 XMIN=+1.0E30 DO 121 I=1,NPOI IF(X(I).GT.XMAX)XMAX=X(I) IF(X(I).LT.XMIN)XMIN=X(I) 121 CONTINUE 122 IF((YMAX-YMIN).GT.1.0E0)GO TO 124 YMAX=-1.0E30 YMIN=+1.0E30 DO 123 I=1,NPOI IF(Y(I).GT.YMAX)YMAX=Y(I) IF(Y(I).LT.YMIN)YMIN=Y(I) 123 CONTINUE 124 CONTINUE SPANX = XMAX-XMIN SPANY = YMAX-YMIN DELX=SPANX/100.0D0 DELY=SPANY/50.0D0 VALUE= YMAX+DELY SMALL(1) = XMIN DO 120 I=1,20 120 SMALL(I+1)=SMALL(I)+5.0D0*DELX DO 135 II=1,NPOI I=(YMAX-Y(II))/DELY+1.5 J = (X(II)-XMIN)/DELX+1.5 IF(I.GT.53.OR.I.LT.1.OR.J.GT.101.OR.J.LT.1)GO TO 135 IF(IABS(ID).NE.2)GO TO 133 ITEM(I,J)=PTID(II) GO TO 135 133 DO 134 JJ=1,10 IF(ITEM(I,J).EQ.HOLL(JJ))GO TO 136 134 CONTINUE IF(ITEM(I,J).EQ.AIE.OR.ITEM(I,J).EQ.AMINUS)ITEM(I,J)=HOLL(2) GO TO 135 136 ITEM(I,J)=HOLL(JJ+1) 135 CONTINUE WRITE(LUWN,3301)DD,(SMALL(I),DD,I=2,20,2) WRITE(LUWN,3302)(SMALL(I),I=1,21,2) WRITE (LUWN,3001) DO 330 I=1,53 VALUE=VALUE-DELY A=BLANK B = PLUS L=I+2 IF(L/3*3-L) 330,310,330 310 B=DD IF(L/2*2-L)330,320,330 320 A=DD 330 WRITE(LUWN,3300)VALUE,A,B,(ITEM(I,J),J=1,101),B,A,VALUE WRITE (LUWN,3001) WRITE(LUWN,3301)DD,(SMALL(I),DD,I=2,20,2) WRITE(LUWN,3302)(SMALL(I),I=1,21,2) IF(ID.LE.0)GO TO 360 DO 340 I=1,53 340 ITEM(I,51)=BLANK DO 350 I=1,101 350 ITEM(27,I)=BLANK 360 DO 370 II=1,NPOI I=(YMAX-Y(II))/DELY+1.5 J=(X(II)-XMIN)/DELX+1.5 IF(I.GT.54.OR.I.LT.1.OR.J.GT.101.OR.J.LT.0)GO TO 370 ITEM(I,J)=BLANK 370 CONTINUE 200 CONTINUE IZINCR=IZINCR-1 RETURN END C SUBROUTINE QTINVT(NRMZ,N,D,E,E2,M,W,IND,Z,IERR, C 1RV1,RV2,RV3,RV4,RV6) C*TITLE: EIGENVECTORS OF A TRIDIAGONAL MATRIX BY INVERSE ITERATION C*PROGRAMMER: EISPACK C*DATE: MAY, 1972 C*LOCATION: ARGONNE NATIONAL LABORATORY C*PURPOSE: DETERMINE THOSE EIGENVECTORS OF A TRIDIAGONAL MATRIX C* *CORRESPONDING TO A SET OF ORDERED EIGENVALUES, USING INVERSE C* *ITERATION. C*SYMBOLS: C* *NRMZ = ROW DIMENSION OF Z IN THE CALLING PROGRAM C* *N = ORDER OF THE TRIDIAGONAL MATRIX C* *D = DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX C* *E = SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX IN THE LAST N-1 C* * * *POSITIONS OF E. E(1) IS ARBITRARY. C* *E2 = SQUARES OF E WITH 0 CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E. C* * * * E2(1) CONTAINS 0 IF EIGENVALUES ARE IN ASCENDING ORDER AND 2 IF C* * * * THEY ARE IN DESCENDING ORDER. C* *M = # OF THE FIRST EIGENVALUES FOR WHICH CORRESPONDING EIGENVECTORS C* * * *ARE OBTAINED. C* *W = CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER C* * * *DEPENDING ON THE VALUE OF E2(1). C* *IND = SUBMATRIX INDICES ASSOCIATED WITH THE M EIGENVALUES IN W. C* * * * *SUBMATRICES ARE NUMBERED FROM 1..., AND ARE DETERMINED BY THE C* * * * *ZERO VALUES IN E2. C* *Z = CONTAINS THE OUTPUT ORTHONORMAL EIGENVECTORS IN AN ORDER C* * * *CORRESPONDING TO THE ORDER OF THE EIGENVALUES. Z IS DIMENSIONED C* * * *NRMZ BY AT LEAST M. C* *IERR = 0 IF EXECUTION TERMINATES NORMALLY. C* * * * * -R IF MORE THAN 5 ITERATIONS ARE REQUIRED FOR AN EIGENVECTOR, C* * * * * WHERE R IS THE INDEX OF THE LAST EIGENVECTOR FOR WHICH THIS C* * * * * OCCURRED. C* *RV1,RV2,RV3,RV4,RV6 = TEMPORARY STORAGE ARRAYS. C*PROCEDURE: FIRST THE E2 ARRAY IS INSPECTED FOR THE PRESENCE OF 0 C* *ELEMENTS DEFINING SUBMATRICES. EIGENVALUES BELONGING TO A GIVEN C* *SUBMATRIX ARE IDENTIFIED BY THEIR COMMON SUBMATRIX INDICES IN IND. C* *EIGENVECTORS OF A SUBMATRIX ARE THEN COMPUTED BY INVERSE ITERATION. C* *FIRST THE LU DECOMPOSITION OF THE SUBMATRIX WITH AN EIGENVALUE C* *SUBTRACTED FROM ITS DIAGONAL ELEMENTS IS ACHIEVED BY GAUSSIAN C* *ELIMINATION WITH PARTIAL PIVOTING. THE MULTIPLIERS DEFINING THE C* *LOWER TRIANGULAR MATRIX L ARE STORED IN RV4 AND THE UPPER C* *TRIANGULAR MATRIX U IS STORED IN RV1,RV2,AND RV3. THUS IF FURTHER C* *ITERATIONS ARE REQUIRED, THE LU DECOMPOSITION NEED NOT BE REPEATED. C* *AN APPROXIMATE VECTOR, STORED IN RV6, IS COMPUTED STARTING FROM AN C* *INITIAL VECTOR, AND THE NORM OF THE APPROXIMATE VECTOR IS COMPARED C* *WITH A NORM OF THE SUBMATRIX TO DETERMINE WHETHER THE GROWTH IS C* *SUFFICIENT TO ACCEPT IT AS AN EIGENVECTOR. IF ACCEPTED, ITS C* *EUCLIDEAN NORM IS MADE 1, IF NOT, THIS VECTOR IS USED AS AN INITIAL C* *VECTOR IN COMPUTING THE NEXT APPROXIMATE VECTOR. THIS ITERATION C* *PROCESS IS REPEATED AT MOST 5 TIMES. EIGENVECTORS COMPUTED IN THIS C* *WAY ARE ORTHOGONAL IF THE CORRESPONDING EIGENVALUES ARE WELL C* *SEPARATED. IF THE EIGENVALUES ARE CLOSE, BUT NOT IDENTICAL, C* *ORTHOGONALITY IS INSURED BY ORTHOGONALIZING EACH APPROXIMATE VECTOR C* *WITH RESPECT TO THE PREVIOUSLY COMPUTED 'CLOSE' VECTORS. IF THIS C* *ORTHOGONALIZATION RESULTS IN A ZERO VECTOR, A COLUMN OF THE IDENTITY C* *MATRIX IS USED AS AN INITIAL VECTOR FOR THE NEXT ITERATION. IF THE C* *EIGENVALUES ARE IDENTICAL, THEY ARE PERTURBED SLIGHTLY AND THESE C* *PERTURBATIONS ARE NOT RECORDED IN W. THE ABOVE STEPS ARE REPEATED C* *ON EACH SUBMATRIX UNTIL ALL REQUESTED EIGENVECTORS ARE COMPUTED. C*REFERENCE: SAME AS QTRBAK. SUBROUTINE QTINVT(NRMZ,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4, +RV6) IMPLICIT DOUBLE PRECISION(A-H,O-Z) INTEGER P,Q,R,S,TAG,GROUP,IND(1) DOUBLE PRECISION MACHEP,NORM,D(1),E(1),E2(1),W(1),Z(NRMZ,1), 2 RV1(1),RV2(1),RV3(1),RV4(1),RV6(1) COMMON/ZINCR/IZINCR CALL XDIMCH(NRMZ,N,'QTINVT',0) CALL QINCR ('QTINVT') C*THE NEXT STEP WAS ADDED TO ALLOW QTINVT TO HANDLE THIS TRIVIAL SITUATI IF(N.GT.1)GOTO1 Z(1,1)=1.D0 IZINCR = IZINCR - 1 RETURN 1 CONTINUE C*MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING THE RELATIVE C* *PRECISION OF FLOATING POINT ARITHMETIC MACHEP=.5D-14 IERR=0 TAG=0 ORDER=1.D0-E2(1) Q=0 C*ESTABLISH AND PROCESS NEXT SUBMATRIX 100 P=Q+1 DO 120Q=P,N IF(Q.EQ.N)GOTO140 IF(E2(Q+1).EQ.0.D0)GOTO140 120 CONTINUE C*FIND VECTORS BY INVERSE ITERATION 140 TAG=TAG+1 S=0 DO 920R=1,M IF(IND(R).NE.TAG)GOTO920 ITS=1 X1=W(R) IF(S.NE.0)GOTO510 C*CHECK FOR ISOLATED ROOT XU=1.D0 IF(P.NE.Q)GOTO490 RV6(P)=1.D0 GOTO 870 490 NORM=DABS(D(P)) IP=P+1 DO 500I=IP,Q 500 NORM=NORM+DABS(D(I))+DABS(E(I)) C*EPS2 IS THE CRITERION FOR GROUPING, EPS3 REPLACES 0 PIVOTS AND EQUAL C* *ROOTS ARE MODIFIED BY EPS3, EPS4 IS TAKEN VERY SMALL TO AVOID C* *OVERFLOW EPS2=1.D-3*NORM EPS3=MACHEP*NORM UK=Q-P+1 EPS4=UK*EPS3 UK=EPS4/DSQRT(UK) S=P 505 GROUP=0 GOTO 520 C*LOOK FOR CLOSE OR COINCIDENT ROOTS 510 IF(DABS(X1-X0).GE.EPS2)GOTO505 GROUP=GROUP+1 IF(ORDER*(X1-X0).LE.0.D0)X1=X0+ORDER*EPS3 C*ELIMINATION WITH INTERCHANGES AND INITIALIZATION OF VECTOR 520 V=0.D0 DO 580I=P,Q RV6(I)=UK IF(I.EQ.P)GOTO560 IF(DABS(E(I)).LT.DABS(U))GOTO540 C*WARNING - A DIVIDE CHECK MAY OCCUR HERE IF E2 ARRAY HAS NOT BEEN C* *SPECIFIED CORRECTLY. XU=U/E(I) RV4(I)=XU RV1(I-1)=E(I) RV2(I-1)=D(I)-X1 RV3(I-1)=0.D0 IF(I.NE.Q)RV3(I-1)=E(I+1) U=V-XU*RV2(I-1) V=-XU*RV3(I-1) GOTO 580 540 XU=E(I)/U RV4(I)=XU RV1(I-1)=U RV2(I-1)=V RV3(I-1)=0.D0 560 U=D(I)-X1-XU*V IF(I.NE.Q)V=E(I+1) 580 CONTINUE IF(U.EQ.0.D0)U=EPS3 RV1(Q)=U RV2(Q)=0.D0 RV3(Q)=0.D0 C*BACK SUBSTITUTION; FOR I=Q STEP -1 UNTIL P DO 600 DO 620II=P,Q I=P+Q-II RV6(I)=(RV6(I)-U*RV2(I)-V*RV3(I))/RV1(I) V=U U=RV6(I) 620 CONTINUE C*ORTHOGONALIZE WITH RESPECT TO PREVIOUS MEMBERS OF GROUP IF(GROUP.EQ.0)GOTO700 J=R DO 680JJ=1,GROUP 630 J=J-1 IF(IND(J).NE.TAG)GOTO630 XU=0.D0 DO 640I=P,Q 640 XU=XU+RV6(I)*Z(I,J) DO 660I=P,Q 660 RV6(I)=RV6(I)-XU*Z(I,J) 680 CONTINUE 700 NORM=0.D0 DO 720I=P,Q 720 NORM=NORM+DABS(RV6(I)) IF(NORM.GE.1.D0)GOTO840 C*FORWARD SUBSTITUTION IF(ITS.EQ.5)GOTO830 IF(NORM.NE.0.D0)GOTO740 RV6(S)=EPS4 S=S+1 IF(S.GT.Q)S=P GOTO 780 740 XU=EPS4/NORM DO 760I=P,Q 760 RV6(I)=RV6(I)*XU C*ELIMINATION OPERATIONS ON NEXT VECTOR; ITERATE 780 DO 820I=IP,Q U=RV6(I) C*IF RV1(I-1).EQ.E(I), A ROW INTERCHANGE WAS PERFORMED EARLIER IN THE C* *TRIANGULARIZATION PROCESS. IF(RV1(I-1).NE.E(I))GOTO800 U=RV6(I-1) RV6(I-1)=RV6(I) 800 RV6(I)=U-RV4(I)*RV6(I-1) 820 CONTINUE ITS=ITS+1 GOTO 600 C*SET ERROR - NONCONVERGED EIGENVECTOR 830 IERR=-R XU=0.D0 GOTO 870 C*NORMALIZE SO THAT SUM OF SQUARES IS 1 AND EXPAND TO FULL ORDER. 840 U=0.D0 DO 860I=P,Q 860 U=U+RV6(I)**2 XU=1.D0/DSQRT(U) 870 DO 880I=1,N 880 Z(I,R)=0.D0 DO 900I=P,Q 900 Z(I,R)=RV6(I)*XU X0=X1 920 CONTINUE IF(Q.LT.N)GOTO100 IZINCR = IZINCR - 1 RETURN END SUBROUTINE TITLE(CH,N) C TITLE: TITLE PRINTER C PROGRAMMER: CARL FINKBEINER C DATE: JUNE, 1974 C LOCATION: UNIVERSITY OF ILLINOIS C PURPOSE: TO PRINT A CHARACTER STRING SPECIFIED BY THE USER. C SYMBOLS: C CH = A CHARACTER STRING OF LENGTH N C N = THE LENGTH OF CH INCLUDING BLANKS AND PUNCTUATION. NOTE C THAT THERE IS NO LIMIT ON THE LENGTH OF CH. C LOGICAL*1 CH(1) IMPLICIT DOUBLE PRECISION(A-H,O-Z) LOGICAL CH(1) COMMON/ZINCR/IZINCR COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG C BOOKKEEPING CALL QINCR ('TITLE') C WRITE CHARACTER STRING N1=N/5 IF(N.NE.N1*5) N1=N1+1 WRITE(LUWN,1)(CH(I),I=1,N1) 1 FORMAT(('0',26A5)) IZINCR = IZINCR - 1 RETURN END SUBROUTINE XTRACE(NRMX,X,N,T) C TITLE: TRACE OF A MATRIX C PROGRAMMER: CARL FINKBEINER C DATE: JUNE, 1974 C LOCATION: UNIVERSITY OF ILLINOIS C PURPOSE: TO FIND THE TRACE (SUM OF THE MAIN DIAGONAL ELEMENTS) OF A C MATRIX C SYMBOLS: C X = INPUT MATRIX C NRMX = NUMBER OF ROWS OF X IN MAIN C N = # OF COLS OF X OR # OF ROWS OF X, WHICH EVER IS SMALLER. C T = TRACE OF X, A SCALAR. C QINCR = PROGRAM COUNTER SUBROUTINE C QDIMCH = DIMENSIONALITY CHECK SUBROUTINE C I = SCRATCH VARIABLE C PROCEDURE: THE MAIN DIAGONAL ELEMENTS OF X ARE SUMMED AND STORED IN T. C NOTE THAT T MAY BE A MATRIX IN THE CALLING PROGRAM, IN WHICH CASE, T C IS RETURNED IN THE FIRST ELEMENT OF THE MATRIX. IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION X(NRMX,1),T COMMON/ZINCR/IZINCR C BOOKKEEPING CALL QINCR ('XTRACE') CALL XDIMCH(NRMX,N,'TRACE',0) C SUM THE DIAGONAL OF X T=0.0D0 DO 1I=1,N 1 T=T+X(I,I) IZINCR = IZINCR - 1 RETURN END SUBROUTINE XTRNSP(NRMX,NRMY,X,Y,NR,NC) C TITLE: MATRIX TRNSPOSE C PROGRAMMER: CARL FINKBEINER C DATE: JUNE,1974 C LOCATION: UNIVERSITY OF ILLINOIS C PURPOSE: TO TRNSPOSE A MATRIX. C SYMBOLS: C X = MATRIX TO BE TRNSPOSED. C Y = THE TRNSPOSE OF X. C NRMX = NUMBER OF ROWS OF X IN MAIN C NRMY = NUMBER OF ROWS OF Y IN MAIN C NR = # OF ROWS OF X = # OF COLS OR Y. C NC = # OF COLS OF X = # OF ROWS OR Y. C QINCR = PROGRAM COUNTER SUBROUTINE. C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE. C I,J = SCRATCH VARIABLES. C PROCEDURE: ROW AND COLUMN SUBSCRIPTS IN X ARE INTERCHANGED WHILE C COPYING THE DATA FROM X TO Y. X AND Y MUST BE DIFFERENT C MATRICES CALLING PROGRAM. IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION X(NRMX,1),Y(NRMY,1) COMMON/ZINCR/IZINCR C BOOKKEEPING CALL QINCR ('XTRNSP') CALL XDIMCH(NRMX,NR,'TRNSP',0) CALL XDIMCH(NRMY,NC,'TRNSP',0) C INTERCHANGE ROW AND COLUMN SUBSCRIPTS FROM X TO Y. DO 1I=1,NR DO 1J=1,NC 1 Y(J,I)=X(I,J) IZINCR = IZINCR - 1 RETURN END SUBROUTINE QTRBAK(NRMA,NRMZ,N,A,E,M,Z) C*TITLE: BACK TRANSFORMATION OF EIGENVECTORS C*PROGRAMMER: EISPACK C*DATE: MAY,1972 C*LOCATION: ARGONNE NATIONAL LABORATORY C*PURPOSE: FORMS THE EIGENVECTORS OF A REAL SYMMETRIC MATRIX FROM THE C* *EIGENVECTORS OF THAT SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRED C*SYMBOLS: C* *NRMA = ROW DIMENSION OF A IN THE CALLING PROGRAM C* *NRMZ = ROW DIMENSION OF Z IN THE CALLING PROGRAM C* *N = ORDER 0F A, MUST BE .LE. NRMA. C* *A = SOME INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS USED IN C* * * *THE REDUCTION TO TRIDIAGONAL FORM (CONTAINED ONLY IN THE STRICT C* * * *LOWER TRIANGLE OF A). THE REMAINDER OF THE MATRIX IS ARBITRARY. C* *E = THE SAME AS E IN QIMTQL. C* *M = THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED. C* *Z = AT LEAST M EIGENVECTORS TO BE BACK TRANSFORMED. ON OUTPUT, Z C* * * *CONTAINS THE TRANSFORMED EIGENVECTORS. C*PROCEDURE: LET C BE THE ORIGINAL SYMMETRIC MATRIX WHICH HAS BEEN C* *REDUCED TO TRIDIAGONAL FORM - F - BY THE SIMILARITY TRANSFORMATION: C* *F=Q'CQ, WHERE Q IS A PRODUCT OF ORTHOGONAL SYMMETRIC MATRICES. THEN C* *GIVEN AN ARRAY Z OF COLUMN VECTORS WHICH ARE EIGENVECTORS OF F, C* *QTRBAK COMPUTES QZ WHIXH IS THE EIGENVECTORS OF C. THIS SUBROUTINE C* *SHOULD BE USED IN CONJUNCTION WITH THE SUBROUTINE QTRED1. C*REFERENCE: WILKINSON,J.H., & REINSCH,C. HANDBOOK FOR AUTOMATIC C* *COMPUTATION, VOL. II, LINEAR ALGEBRA, PART2, SPRINGER-VERLAG, C* *NEW YORK, HEIDELBERG, BERLIN, 1971. C* *SMITH,B.T., BOYLE,J.M., GARBOW,B.S., IKEBE,Y., KLEMA,V.C., MOLER,C. C* *MATRIX EIGENSYSTEM ROUTINES - EISPACK GUIDE IN LECTURE NOTES IN C* *COMPUTER SCIENCE, EDITED BY G. GOOS AND J. HARTMANIS, SPRINGER- C* *VERLAG, NEW YORK, BERLIN, HEIDELBERG, 1974. IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION A(NRMA,1),E(1),Z(NRMZ,1) COMMON/ZINCR/IZINCR CALL QINCR ('QTRBAK') CALL XDIMCH(NRMA,N,'QTRBAK',0) CALL XDIMCH(NRMZ,N,'QTRBAK',0) IF(N.EQ.1)IZINCR = IZINCR - 1 IF(N.EQ.1)RETURN DO 140 I=2,N L=I-1 C*H BELOW IS NEGATIVE OF H FORMED IN QTRED1. H=E(I)*A(I,L) IF(H.EQ.0.D0) GOTO 140 DO 130 J=1,M S=0.D0 DO 110 K=1,L 110 S=S+A(I,K)*Z(K,J) S=S/H DO 120 K=1,L 120 Z(K,J)=Z(K,J)+S*A(I,K) 130 CONTINUE 140 CONTINUE IZINCR = IZINCR - 1 RETURN END SUBROUTINE QTRED1(NRMA,N,A,D,E,E2) C*TITLE: REDUCE A SYMMETRIC MATRIX TO TRIDIAGONAL FORM C*PROGRAMMER: EISPACK C*DATE: MAY,1972 C*LOCATION: ARGONNE NATIONAL LABORATORY C*SYMBOLS: C* *NRMA = ROW DIMENSION OF A IN THE CALLING PROGRAM. C* *N = THE ORDER OF A, MUST BE .LE. NRMA. C* *A = THE REAL SYMMETRIC MATRIX TO BE TRIDIAGONALIZED. ONLY THE FULL C* * * *LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. ON OUTPUT,THE C* * * *STRICT LOWER TRIANGLE CONTAINS TRANSFORMATION INFORMATION. C* *D = THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL FORM OF A. C* *E = SAME AS E IN QIMTQL. C* *E2 = SAME AS E2 IN QTINVT. E2 AND E NEED NOT BE DISTINCT IN THE C* * * * CALLING PROGRAM. C*PROCEDURE: THE TRIDIAGONAL REDUCTION IS PERFORMED IN THE FOLLOWING C* *WAY. STARTING WITH J=N, THE ELEMENTS IN THE JTH ROW TO THE LEFT OF C* *THE DIAGONAL ARE FIRST SCALED, TO AVOID POSSIBLE UNDERFLOW IN THE C* *TRANSFORMATION THAT MIGHT RESULT IN SEVERE DEPARTURE FROM C* *ORTHOGONALITY. THE SUM OF SQUARES SIGMA OF THESE SCALED ELEMENTS IS C* *NEXT FORMED. THEN, A VECTOR U AND A SCALAR H (=U'U/2) DEFINE AN C* *OPERATOR P (=I-UU'/H) WHICH IS ORTHOGONAL AND SYMMETRIC AND FOR C* *WHICH THE SIMILARITY TRANSFORMATION PAP ELIMINATES THE ELEMENTS IN C* *THE JTH ROW OF A TO THE LEFT OF THE SUBDIAGONAL AND THE SYMMETRICAL C* *ELEMENTS IN THE JTH COLUMN. THE NON-ZERO COMPONENTS OF U ARE THE C* *ELEMENTS OF THE JTH ROW TO THE LEFT OF THE DIAGONAL WITH THE LAST C* *OF THEM AUGMENTED BY THE SQUARE ROOT OF SIGMA PREFIXED BY THE SIGN C* *OF THE SUBDIAGONAL ELEMENT. BY STORING THE TRANSFORMED SUBDIAGONAL C* *ELEMENT IN E(J) AND NOT OVERWRITING THE ROW ELEMENTS ELIMINATED IN C* *THE TRANSFORMATION, FUL INFORMATION ABOUT P IS SAVED FOR LATER USE C* *IN QTRBAK. THE TRANSFORMATION SETS E2(J) EQUAL TO SIGMA AND E(J) C* *EQUAL TO THE SQUARE ROOT OF SIGMA PREFIXED BY SIGN OPPOSITE TO THAT C* *OF THE REPLACED SUBDIAGONAL ELEMENT. THE ABOVE STEPS ARE REPEATED C* *ON FURTHER ROWS OF THE TRANSFORMED A IN REVERSE ORDER UNTIL A IS C* *REDUCED TO TRIDIAGONAL FORM; THAT IS, REPEATED J=N-1,N-2,...,3. C* *ONLY THE ELEMENTS IN THE LOWER TRIANGLE OF A ARE ACCESSED, AND C* *ALTHO THE DIAGONAL ELEMENTS ARE MODIFIED IN THE ALGORITHM, THEY ARE C* *RESTORED TO THEIR ORIGINAL CONTENTS BY THE END OF THE SUBROUTINE, C* *THUS PRESERVING THE FULL UPPER TRIANGLE OF A. C*REFERENCES: SAME AS QTRBAK IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION A(NRMA,1),D(1),E(1),E2(1) COMMON/ZINCR/IZINCR CALL QINCR ('QTRED1') CALL XDIMCH(NRMA,N,'QTRED1',0) DO 100I=1,N 100 D(I)=A(I,I) C*FOR I=N STEP -1 UNTIL 1 DO DO 300II=1,N I=N+1-II L=I-1 H=0.D0 SCALE=0.D0 IF(L.LT.1)GOTO130 C*SCALE ROW DO 120K=1,L 120 SCALE=SCALE+DABS(A(I,K)) IF(SCALE.NE.0.D0)GOTO140 130 E(I)=0.D0 E2(I)=0.D0 GOTO 290 140 DO 150K=1,L A(I,K)=A(I,K)/SCALE H=H+A(I,K)*A(I,K) 150 CONTINUE E2(I)=SCALE*SCALE*H F=A(I,L) G=-SIGN(DSQRT(H),F) E(I)=SCALE*G H=H-F*G A(I,L)=F-G IF(L.EQ.1)GOTO270 F=0.D0 DO 240J=1,L G=0.D0 C*FORM ELEMENT OF A*U DO 180K=1,J 180 G=G+A(J,K)*A(I,K) JP1=J+1 IF(L.LT.JP1)GOTO220 DO 200K=JP1,L 200 G=G+A(K,J)*A(I,K) C*FORM ELEMENT OF P 220 E(J)=G/H F=F+E(J)*A(I,J) 240 CONTINUE H=F/(H+H) C*FORM REDUCED A DO 260J=1,L F=A(I,J) G=E(J)-H*F E(J)=G DO 260K=1,J A(J,K)=A(J,K)-F*E(K)-G*A(I,K) 260 CONTINUE 270 DO 280K=1,L 280 A(I,K)=SCALE*A(I,K) 290 H=D(I) D(I)=A(I,I) A(I,I)=H 300 CONTINUE IZINCR = IZINCR - 1 RETURN END SUBROUTINE XVRTCL(NRMX,NRMY,NRMZ,X,Y,Z,NRX,NRY,NC) C TITLE: VERTICAL CONCATENATE C PROGRAMMER: CARL FINKBEINER C DATE: JUNE,1974 C LOCATION: UNIVERSITY OF ILLINOIS C PURPOSE: TO CONCATENATE THE COLS 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) C (Y) C NRMX = NUMBER OF ROWS OF X IN MAIN C NRMY = NUMBER OF ROWS OF Y IN MAIN C NRMZ = NUMBER OF ROWS OF Z IN MAIN C NRY = # OF ROWS IN Y C NRX = # OF ROWS IN X C NC = # OF COLS IN BOTH X AND Y C QINCR = PROGRAM COUNTER SUBROUTINE C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE C NR = # OF ROWS IN Z = NRX + NRY C I,J,K = SCRATCH VARIABLES C PROCEDURE: C A MATRIX Z = (X..Y) IS CREATED BY STORING X IN COLS 1 THRU NC, C ROWS 1 THRU NRX OF Z AND THEN BY STORING Y IN COLS 1 THRU NC AND ROW C (NRX+1) THRU NR OF Z. THUS Z WILL BE NR X NC. ALL 3 MATRICES, X,Y, C AND Z MAY BE THE SAME MATRIX IN THE CALLING PROGRAM; MATRICES C X AND Y MAY BE THE SAME BUT DIFFERENT FROM Z; MATRICES X AND Z MAY C BE THE SAME BUT DIFFERENT FROM Y; HOWEVER, MATRICES Y AND Z MAY C 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 ('XVRTCL') CALL XDIMCH(NRMX,NRX,'VRTCL',0) CALL XDIMCH(NRMY,NRY,'VRTCL',0) NR=NRX+NRY CALL XDIMCH(NRMZ,NR,'VRTCL',1) C LOOP FOR SETTING UPPER MOST PORTION OF Z EQUAL TO X. DO 1I=1,NRX DO 1J=1,NC 1 Z(I,J)=X(I,J) C LOOP FOR SETTING NEXT PORTION DOWN IN Z EQUAL TO Y. K=NRX DO 2I=1,NRY K=K+1 DO 2J=1,NC 2 Z(K,J)=Y(I,J) IZINCR = IZINCR - 1 RETURN END SUBROUTINE QWRGDM(NR,NC) C TITLE: ERROR MESSAGE ROUTINE FOR INVRS,SIMLN, AND DET C PROGRAMMER: CARL FINKBEINER C DATE: JUNE, 1974 C LOCATION: UNIVERSITY OF ILLINOIS C SYMBOLS: C NR = # OF ROWS OF THE MATRIX INPUT TO INVRS. C NC = # OF COLS OF THE MATRIX INPUT TO INVRS. C A = LABELLED COMMON AREA CONTAINING IPGM. C IPGM = PROGRAM COUNTER C PROCEDURE: EITHER NC OR BOTH NR AND NC MAY BE SET AT -9999 BY INVRS. C SUCCESSIVE TESTS ON NC AND NR DETERMINES THE APPROPRIATE ERROR C MESSAGE PRINTED OUT AND EXECUTION TERMINATES. IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/A/IPGM,ISUBRU COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG COMMON/ZINCR/IZINCR CALL QINCR ('QWRGDM') IF(NC.EQ.-9999)GOTO2 C IF BOTH DIMENSIONS ARE CORRECTLY PRESENT, PRINT MESSAGE CONCERNING NR C .GT. NC AND TERMINATE. CALL END (1) WRITE(LUWN,1)NR,NC 1 FORMAT('0* * * * * THE NUMBER OF ROWS (=',I4,') OF THE INPUT MATRI 2X TO SIMLN MUST BE LESS THAN THE NUMBER OF COLUMNS (=',I4,').') CALL END (2) CALL END (3) 2 IF(NR.EQ.-9999)GOTO4 C IF ONLY NR IS CORRECTLY PRESENT, PRINT MESSAGE CONCERNING NR BEING .EQ C AND TERMINATE. CALL END (1) WRITE(LUWN,3) 3 FORMAT('0* * * * * THE DIMENSIONALITY OF THE MATRIX MUST BE GREATE 3R THAN 1') CALL END (2) CALL END (3) C IF BOTH NR AND NC ARE NOT CORRECTLY PRESENT, PRINT MESSAGE CONCERNING C SINGULARITY AND TERMINATE. 4 CALL END (1) WRITE(LUWN,5) 5 FORMAT('0* * * * * THE MATRIX PASSED WAS SINGULAR NO INVERSE ', 2'WAS POSSIBLE.') CALL END (2) CALL END (3) END SUBROUTINE XXTX(NRMX,NRMY,X,Y,NR,NC) C TITLE: MATRIX CROSS-PRODUCTS C PROGRAMMER: CARL FINKBEINER C DATE: JUNE, 1974 C LOCATION: UNIVERSITY OF ILLINOIS C PURPOSE: TO FIND A FAST SOLUTION TO THE PRE-MULTIPLICATION OF A MATRIX C TRNSPOSE. C SYMBOLS: C X = INPUT MATRIX C Y = RESULT MATRIX = X'X C NRMX = NUMBER OF ROWS OF X IN MAIN C NRMY = NUMBER OF ROWS OF Y IN MAIN C NR = # OF ROWS OF X C NC = # OF COLS OF X = # OF ROWS AND COLS OF Y, MUST BE .LT. NRMX & NRMY C A = LABELLED COMMON AREA CONTAINING IPGM C IPGM = PROGRAM COUNTER C QINCR = PROGRAM COUNTER SUBROUTINE C XDIMCH = DIMENSIONALITY CHECK SUBROUTINE C I,J,K,N = SCRATCH VARIABLES C PROCEDURE: THE FIRST COL OF X IS DOT MULTIPLIED BY EACH SUCCESSIVE COL C OF X AND THE RESULTS ARE STORED IN Z. Z IS COPIED INTO THE FIRST CO C THE 2ND COL OF X IS DOT MULTIPLIED BY EACH OF THE SUCCEEDING COLS OF C THE RESULTS ARE STORED IN THE LAST (NC-1) ELEMENTS OF Z. THESE ELEM C ARE COPIED INTO THE LAST (NC-1) ELEMENTS C OF THE 2ND COL OF Y AND THE FIRST ELEMENT OF THE 2ND COL IS MADE A C C THE CORRESPONDING ELEMENT ON THE OTHER SIDE OF THE DIAGONAL. THIS P C IS REPEATED UNTIL THE PRODUCT IS COMPLETED. THE PROCEDURE USED HERE C AS FAST AS IT MIGHT BE BECAUSE OF THE COPY FROM Z TO Y, BUT IT IS CE C FASTER THAN THE SUCCESSIVE USE OF TRNSP AND MATML. THE PRESENT PR C WAS CHOSEN BECAUSE IT ALLOWS THE FURTHER POTENTIAL EFFICIENCY OF X A C BEING THE SAME MATRIX IN THE CALLING PROGRAM. IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION X(NRMX,NC),Y(NRMY,NC) COMMON/A/IPGM,ISUBRU COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG COMMON/ZINCR/IZINCR C BOOKKEEPING NCP1 = NC+1 CALL QINCR ('XXTX') CALL XDIMCH(NRMY,NCP1,'XTX',0) CALL XDIMCH(NRMX,NR,'XTX',0) CALL XDIMCH(NRMX,NC,'XTX',0) IF(NC.EQ.NRMX)GOTO4 C LOOP FOR FINDING CROSS-PRODUCTS DO 1I=1,NC C LOOP FOR STORING CROSS-PRODUCTS OF A COL WITH ITS SUCCESSORS IN Z DO 2J=I,NC X(J,NRMX)=0.D0 DO 2K=1,NR 2 X(J,NRMX)=X(J,NRMX)+X(K,I)*X(K,J) C LOOP FOR COPYING Z INTO LOWER TRIANGULAR PORTION OF A COL OF Y DO 3J=I,NC 3 Y(J,I)=X(J,NRMX) C LOOP FOR COPYING LOWER TRIANGLE INTO UPPER TRIANGULAR PORTION UP TO TH C CURRENT COLUMN. N=I-1 DO 1J=1,N 1 Y(J,I)=Y(I,J) DO 25 J = 1,NC Y(NCP1,J)=0.0D0 DO 25 I = 1,NR 25 Y(NCP1,J)=Y(NCP1,J)+X(I,J) IZINCR = IZINCR - 1 RETURN 4 CALL END (1) WRITE(LUWN,5)NC,NRMX 5 FORMAT('0* * * * * THE NUMBER OF COLUMNS (=',I4,') OF THE INPUT MA 1TRIX '/20X,' FOR WHICH YOU ARE FINDING CROSS-PRODUCTS MUST BE LESS 2 THAN THE DIMENSION OF THE MATRIX (=',I4,').') CALL END (2) CALL END (3) END FUNCTION QRAND1(T) C PROVIDES A PSEUDO RANDOM NORMAL DEVIATE (MEAN = 0, SD = 1) C BY INVERSE TRANSFORMATION OF A PSEUDO RANDOM P PROVIDED BY FUNCTION C RAN; FUNCTION QDVNRM IS USED FOR THE INVERSE TRANSFORMATION, QDVNRM C GIVES A QUICK APPROXIMATION TO THE INVERSE TRANSFORMATION C A 'SEED' INTEGER IS TO BE GIVEN FUNCTION RAN BEFORE THE FIRST USE C OF RAN BY SUBROUTINE SETRAN(IRN) WHERE IRN IS A 9 DIGIT INTEGER. C THE FINAL INTEGER MAY BE OBTAINED BY SUBROUTINE C SAVRAN(IRN). C T IS A DUMMY ARGUMENT. P = RAN(T) QRAND1 = QDVNRM(P) RETURN END SUBROUTINE RXTX(P,NIN,NAT) C FOR A MATRIX X READ FROM CARDS COMPUTES PRODUCT MATRIX P = X'X C COLUMN SUMS OF X ARE IN ROW NAT+1 OF P C PROGRAM WRITTEN BY LEDYARD R TUCKER, JANUARY 1978 C ROW DIMENSION OF ARRAY FOR P IS NRMP IN COMMON BLOCK B AS IN FORMAL C AND MUST BE AT LEAST NAT + 1 C COLUMN DIMENSION OF ARRAY FOR P MUST BE AT LEAST NAT + 1 C P IS MATRIX CONTAINING RESULTS; OUTPUT C NIN IS NUMBER OF OBJECTS (ROWS OF X); READ FROM HEADER CARD FOR C DATA DECK; OUTPUT C NAT IS NUMBER OF ATTRIBUTES (COLUMNS OF X); READ FROM HEADER CARD C FOR DATA DECK; OUTPUT C MATRIX X IS TO BE CONTAINED IN A DECK OF DATA CARDS WITH A RECORD OF C ONE OR MORE CARDS FOR EACH OBJECT (ROW OF X) C A HEADER CARD IS TO PRECEDE THE DATA DECK C COLUMNS 1 - 4: NIN C COLUMNS 5 - 8: NAT C COLUMNS 9 - 80: FORMAT OF DATA CARDS IN PARENTHESES C DATA TITLE MAY FOLLOW THE FORMAT AND WILL BE PRINTED IMPLICIT DOUBLE PRECISION (A-H, O-Z) DOUBLE PRECISION P(I0,1) COMMON/B/I0 COMMON/ZINCR/IZINCR NRMP = I0 CALL QINCR ('RXTX') CALL XRXTX(NRMP,P,NIN,NAT) IZINCR = IZINCR - 1 RETURN END SUBROUTINE XRXTX(NRMP,P,NIN,NAT) IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION P(NRMP,1),MFI(18) COMMON/ZINCR/IZINCR COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG CALL QINCR ('XRXTX') READ(LURN,900)NIN,NAT,(MFI(I),I=1,18) 900 FORMAT(2I4,18A4) C INITIALIZE NATP1 = NAT + 1 NMAT = NMAT + 1 CALL XDIMCH(NRMP,NATP1,'RXTX',0) DO 5 J = 1,NATP1 DO 5 K = 1,NAT 5 P(J,K) = 0.0D0 C READ EACH ROW OF X IN SUCCESSION UPDATING P AND SUMS DO 10 I = 1,NIN READ(LURN,MFI)(P(J,NATP1),J=1,NAT) DO 10 K = 1,NAT P(NATP1,K) = P(NATP1,K) + P(K,NATP1) DO 10 J = 1,K 10 P(J,K) = P(J,K) + P(J,NATP1)*P(K,NATP1) DO 15 K = 1,NAT DO 15 J = 1,K 15 P(K,J) = P(J,K) WRITE(LUWN,910) NMAT, LURN, NIN, NAT, MFI, XINPUT(1) 910 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) IZINCR = IZINCR - 1 RETURN END SUBROUTINE XCHSYM(NRMA,A,N) C CHECKS SYMMETRY OF MATRIX A C PROGRAM WRITTEN BY LEDYARD R TUCKER, APRIL 1978 C NRMA IS NUMBER OF ROWS IN MAIN OF ARRAY CONTAINING A; INPUT/OUTPUT C A IS MATRIX TO BE CHECKED; INPUT/OUTPUT C N IS ORDER OF MATRIX; INPUT/OUTPUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION A(NRMA,1) COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG COMMON/ZINCR/IZINCR CALL QINCR ('XCHSYM') CALL XDIMCH(NRMA,N,'CHSYM',0) 910 FORMAT('0','CHECK ON SYMMETRY OF A MATRIX') 911 FORMAT('0','THE ENTRIES IN THE FOLLOWING CELLS ARE NOT EQUAL') 912 FORMAT('0','CELL',I4,',',I4,5X,'WITH CELL',I4,',',I4) 913 FORMAT('0','NO DISCREPANCIES WERE LOCATED') WRITE(LUWN,910) DO 5 I = 2,N IM = I - 1 DO 5 J = 1,IM IF(A(I,J).NE.A(J,I)) GO TO 10 5 CONTINUE WRITE(LUWN,913) GO TO 20 10 WRITE(LUWN,911) DO 15 I = 2,N IM = I - 1 DO 15 J = 1,IM IF(A(I,J).NE.A(J,I)) WRITE(LUWN,912) I,J,J,I 15 CONTINUE 20 IZINCR = IZINCR - 1 RETURN END SUBROUTINE PRDIF(V,N,MD) C PRINTS A COLUMN OF VALUES AND FIRST DIFFERENCES C THESE VALUES MAY BE EITHER ENTRIES IN A VECTOR OR DIAGONAL ENTRIES OF A C MATRIX. A COMMON APPLICATION IS PRINTING A SERIES OF EIGENVALUES. C PROGRAM WRITTEN BY LEDYARD R TUCKER, APRIL 1978 C V IS ARRAY CONTAINING VALUES TO BE PRINTED; INPUT/OUTPUT C N IS NUMBER OF VALUES; INPUT/OUTPUT C MD IS A FLAG CONCERNING LOCATION OF VALUES IN V; INPUT/OUTPUT C = 0 FOR VECTOR (INCLUDES A COLUMN OF A TWO DIMENSIONAL ARRAY) C = NUMBER OF ROWS IN MAIN OF ARRAY V WHEN THE DIAGONAL ENTRIES C OF V ARE TO BE PRINTED. IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION V(1) COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG COMMON/ZINCR/IZINCR CALL QINCR ('PRDIF') 910 FORMAT('0',14X,'VALUE',6X,'DIFFERENCE') 911 FORMAT(' ',I6,F14.5) 912 FORMAT(' ',20X,F14.5) WRITE(LUWN,910) WRITE(LUWN,911) IC = 0 IP = 1 5 I = IP IC = IC + 1 WRITE(LUWN,911)IC,V(I) IF(IC.GE.N) GO TO 10 IP = IP + MD + 1 T = V(I) - V(IP) WRITE(LUWN,912) T GO TO 5 10 IZINCR = IZINCR - 1 RETURN END SUBROUTINE QLNEQT(NRMA,A,B,E,S,D,F,BP,N,MS,IERR,NCYC) C SOLVES FOR B IN LINEAR EQUATION A*B = C C BY AN ITERATIVE PROCEDURE C MAY BE USED TO IMPROVE AN APPROXIMATE SOLUTION C A IS AN N X N MATRIX OF COEFFICIENTS; INPUT/OUTPUT C B IS A VECTOR OF VALUES OF VARIABLES; C INITIAL VALUES INPUT/RESULTING VALUES OUTPUT C (WHEN AN INITIAL VECTOR IS NOT AVAILABLE, SET B = 0) C E IS A VECTOR TO CONTAIN C ON INPUT C CONTAINS DISCREPANCIES (E = C - AB) ON OUTPUT C S IS A VECTOR OF SUMS OF SQUARES OF COLUMNS OF A C NATURE OF INPUT INDICATED BY MS C MS = 0 WHEN S IS NOT INPUT (S TO BE COMPUTED) C = 1 WHEN S IS INPUT C S IS OUTPUT AND MS = 1 C D IS A SCRATCH VECTOR C F IS A SCRATCH VECTOR C BP IS A SCRATCH VECTOR C N IS ORDER OF A, B, E, S; INPUT/OUTPUT C NRMA IS NUMBER OF ROWS IN MAIN OF ARRAY CONTAINING A; INPUT/OUTPUT C IERR IS COMPLETION INDICATOR; OUPUT C = 0 WHEN SOLUTION IS OBTAINED WITHIN NCYC CYCLES C = 1 OTHERWISE C NCYC IS NUMBER OF ITERATIVE CYCLES; INPUT/OUTPUT C DEFAULT = 20 C IN CASE A IS SINGULAR, A LEAST SQUARES SOLUTION IS OBTAINED C IERR = 0 FOR AN EXACT SOLUTION C IERR = 2 FOR A LEAST SQUARES APPROXIMATION TO C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION A(NRMA,1),B(1),E(1),S(1),D(1),F(1),BP(1) COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG COMMON/ZINCR/IZINCR CALL QINCR ('QLNEQT') CALL XDIMCH(NRMA,N,'QLNEQT',0) NC = NCYC IF(NC.EQ.0) NC = 20 NSC = 4 EPS = 1.0D-10 EPS2 = 1.0D-35 IERR = 0 IC = 0 ISC = 0 C COMPUTE S WHEN NECESSARY IF(MS.NE.0) GO TO 10 MS = 1 DO 5 J = 1,N S(J) = 0.0D0 DO 5 I = 1,N 5 S(J) = S(J) + A(I,J)**2 10 CONTINUE C STANDARDIZE C CMAX = 0.0D0 DO 15 I = 1,N 15 IF(ABS(E(I)).GT.CMAX) CMAX = ABS(E(I)) DO 20 I = 1,N B(I) = B(I)/CMAX 20 E(I) = E(I)/CMAX C INITIAL E DO 25 I = 1,N DO 25 J = 1,N IF(B(J).EQ.0.0D0) GO TO 25 E(I) = E(I) - A(I,J)*B(J) 25 CONTINUE C START OF ITERATION CYCLE 100 IC = IC + 1 DO 102 I = 1,N 102 BP(I) = B(I) C START OF SUBCYCLE 103 ISC = ISC + 1 C TEST E DO 105 I = 1,N IF(ABS(E(I)).GE.EPS) GO TO 108 105 CONTINUE GO TO 125 108 CONTINUE IF(ISC.GT.NSC) GO TO 115 C GRADIENT STEP DO 112 I = 1,N D(I) = 0.0D0 DO 111 J = 1,N 111 D(I) = D(I) + A(J,I)*E(J) 112 D(I) = D(I)/S(I) GO TO 118 115 DO 116 I = 1,N 116 D(I) = B(I) - BP(I) ISC = 0 118 T1 = 0.0D0 T2 = 0.0D0 DO 120 I = 1,N F(I) = 0.0D0 DO 119 J = 1,N 119 F(I) = F(I) + A(I,J)*D(J) T1 = T1 + F(I)**2 120 T2 = T2 + F(I)*E(I) IF(T1.LT.EPS2) GO TO 123 T2 = T2/T1 DO 122 I = 1,N B(I) = B(I) + T2*D(I) 122 E(I) = E(I) - T2*F(I) IF(ISC.NE.0) GO TO 103 C TEST IC AGAINST NC AND RETURN TO START OF CYCLE IF(IC.LT.NC) GO TO 100 IERR = 1 WRITE(LUWN,900) NC 900 FORMAT('-','WARNING: CONVERGENCE NOT OBTAINED IN',I4,3X, 1'CYCLES OF IMPROVEMENT OF SOLUTION') GO TO 125 123 IERR = 2 WRITE(LUWN,902) 902 FORMAT('0','A LEAST SQUARES SOLUTION WAS OBTAINED') 125 DO 130 I = 1,N B(I) = B(I)*CMAX 130 E(I) = E(I)*CMAX IZINCR = IZINCR - 1 RETURN END SUBROUTINE RNKORD(V,S,IRP,IRO,N) C DETERMINES RANK ORDER (LOW TO HIGH) OF ENTRIES IN VECTOR V C PROGRAMMED BY LEDYARD R TUCKER; NOVEMBER, 1977 C V IS A REAL VECTOR; INPUT/OUTPUT C S IS A REAL VECTOR CONTAINING ON OUTPUT THE ENTRIES IN V C SORTED FROM LOW TO HIGH C IRP IS AN INTEGER VECTOR CONTAINING ON OUTPUT THE RANK POSITION C OF ENTRIES IN V C IRO IS AN INTEGER VECTOR CONTAINING ON OUTPUT ENTRY IDENTIFICATIONS C WHEN ENTRIES ARE SORTED INTO ASCENDING ORDER C N IS THE NUMBER OF ENTRIES IN V ; INPUT/OUTPUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION V(1),S(1),IRP(1),IRO(1) COMMON/ZINCR/IZINCR CALL QINCR ('RNKORD') C SORT ENTRIES IN V INTO ASCENDING ORDER IN S C WITH ORIGINAL POSITION IN IRO S(N) = V(N) IRO(N) = N NM1 = N - 1 DO 20 M = 1,NM1 MB = N - M T1 = V(MB) IT1 = MB DO 15 K = MB,NM1 KP1 = K + 1 IF(T1.GT.S(KP1)) GO TO 10 S(K) = T1 IRO(K) = IT1 GO TO 20 10 S(K) = S(KP1) 15 IRO(K) = IRO(KP1) S(N) = T1 IRO(N) = IT1 20 CONTINUE C INSERT RANK POSITION IN IRP DO 25 I = 1,N K = IRO(I) 25 IRP(K) = I IZINCR = IZINCR - 1 RETURN END SUBROUTINE XXXT(NRMX,NRMY,X,Y,NR,NC) C OBTAINS MATRIX PRODUCT Y = XX' C NRMX IS NUMBER OF ROWS IN MAIN OF ARRAY CONTAINING X; INPUT/OUTPUT. C NRMY IS NUMBER OF ROWS IN MAIN OF ARRAY TO CONTAIN Y; INPUT/OUTPUT. C X IS DATA MATRIX; INPUT/OUTPUT. C Y IS PRODUCT MATRIX; OUPUT. C NR IS NUMBER OF ROWS OF X AND ORDER OF Y; INPUT/OUTPUT. C NC IS NUMBER OF COLUMNS OF X; INPUT/OUTPUT. IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION X(NRMX,NC),Y(NRMY,NR) COMMON/ZINCR/IZINCR CALL QINCR ('XXXT') CALL XDIMCH(NRMX,NR,'XXT',0) CALL XDIMCH(NRMY,NR,'XXT',0) DO 10 I = 1,NR DO 10 K = 1,I Y(I,K) = 0.0D0 DO 5 J = 1,NC 5 Y(I,K) = Y(I,K) + X(I,J)*X(K,J) 10 Y(K,I) = Y(I,K) IZINCR = IZINCR - 1 RETURN END SUBROUTINE RNDVC(RVEC,ISRT,N,IRN1,IRN2) C GENERATES A VECTOR OF RANDOM NORMAL DEVIATES, MEAN = 0, S.D. = 1 C RANDOMLY SORTED C RVEC CONTAINS THE RANDOM NORMAL DEVIATES ON OUTPUT C ISRT IS A VECTOR OF RANDOM INTEGERS USED IN SORTING C N IS THE NUMBER OF ENTRIES IN RVEC AND ISRT C IRN1 AND IRN2 ARE TWO 9 DIGIT, ODD INTEGERS USED IN C GENERATING THE RANDOM DEVIATES AND INTEGERS C THEIR INITIAL VALUES ARE TO BE INPUT, THEIR VALUES ARE C CHANGED RANDOMLY BEFORE OUTPUT IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION RVEC(1) INTEGER ISRT(1) COMMON/ZINCR/IZINCR CALL QINCR ('RNDVC') C GENERATE RANDOM NORMAL DEVIATES IN RVEC DO 18 I=1,N IRN1 = IRN1*65539 IF(IRN1)12,13,13 12 IRN1 = IRN1 + 2147483647 + 1 13 RANP = IRN1 RANP = RANP*.4656613D-9 RANP = RANP - .5 TEMP = DABS(RANP) IF(TEMP - .450)14,14,15 14 RVEC(I) = TEMP*(2.45D0-3.746585D0*TEMP)/(1.0D0-1.757734D0*TEMP) GO TO 16 15 RVEC(I) = -2.074847D0+TEMP*(12.120613D0-15.794816D0*TEMP) RVEC(I) = RVEC(I)/(1.0D0-1.977724D0*TEMP) 16 IF(RANP)17,18,18 17 RVEC(I) = -RVEC(I) 18 CONTINUE C GENERATE RANDOM INTEGERS IN ISRT DO 23 I=1,N IRN2 = IRN2*65539 IF(IRN2)22,23,23 22 IRN2 = IRN2 + 2147483647 + 1 23 ISRT(I) = IRN2 C ADVANCE IRN2 BY 5 LOCATIONS TO PREVENT SYNCHRONOUS C REPETITION WITH IRN1 FOR SUBSEQUENT CYCLES OF GENERATOR DO 26 I=1,5 IRN2 = IRN2*65539 IF(IRN2)25,26,26 25 IRN2 = IRN2 + 2147483647 + 1 26 CONTINUE C SORT RVEC BY ORDER OF ISRT NRMZ = N-1 DO 34 I=1,NRMZ IB = N-I TEMP = RVEC(IB) ITEM = ISRT(IB) DO 32 J=IB,NRMZ JP = J+1 IF(ITEM.LE.ISRT(JP)) GO TO 33 RVEC(J) = RVEC(JP) 32 ISRT(J) = ISRT(JP) RVEC(N) = TEMP ISRT(N) = ITEM GO TO 34 33 RVEC(J) = TEMP ISRT(J) = ITEM 34 CONTINUE IZINCR = IZINCR - 1 RETURN END SUBROUTINE QMINVZ(NRMAA,NRMA,NRMIDU,AA,N,DET,IDET,A,IDUM,IER,V) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION AA(NRMAA,1),A(NRMA,1),IDUM(NRMIDU,2),V(1) COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG COMMON/ZINCR/IZINCR CALL QINCR ('QMINVZ') CALL XDIMCH(NRMAA,N,'QMINVZ',0) CALL XDIMCH(NRMA,N,'QMINVZ',0) CALL XDIMCH(NRMIDU,N,'QMINVZ',0) C----------------------------------------------------------------------- C COMPUTES THE INVERSE OF THE N BY N REAL, GENERAL INPUT MATRIX BY C GAUSS ELIMINATION. THE DETERMINANT IS CALCULATED AS A BY-PRODUCT C OF THE COMPUTATION. C PARAMETERS ARE: C AA - THE N BY N COEFFICIENT MATRIX, C WITH N.LE.EACH OF NRMAA,NRMA,AND NRMIDU C DET - THE DETERMINANT IS RETURNED IN THE FORM: C IDET DET*10.0**IDET, WHERE 1.0.LE.DET.LT.10.0 OR DET=0.0 . C WHEN DET=0.0 THEN THE INPUT MATRIX IS SINGULAR AND NO C INVERSE IS COMPUTED. C A - SCRATCH ARRAY USED BY QMINVZ. IT IS PERMISSIBLE FOR A AND AA C TO SHARE THE SAME STORAGE IN WHICH CASE AA IS OVERWRITTEN BY C THE INVERSE. IF A AND AA ARE DIFFERENT IN STORAGE, THEN AA C IS RETAINED IN THE COMPUTATION AND A CONTAINS THE INVERSE. C IDUM- SCRATCH STORAGE. 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 OR A PIVOT ELEMENT WAS = 0.D0. DET IS SET TO C 0.D0 AND NO INVERSE IS COMPUTED. C = K QMINVZ DETECTED A LOSS OF SIGNIFICANT DIGITS IN THE COM- C PUTED INVERSE. THE K'TH PIVOT ELEMENT (K=1,N-1) WAS C LESS THAN 1.D-15 TIMES THE LARGEST ELEMENT OF THE INPUT C MATRIX. THE INPUT MATRIX IS PROBABLY SINGULAR. C----------------------------------------------------------------------- IF(N.LE.1)GOTO140 IER1=IER IER=0 DO 10I=1,N DO 10J=1,N 10 A(I,J)=AA(I,J) C----------------------------------------------------------------------- C USE GAUSSIAN ELIMINATION TO DECOMPOSE THE INPUT MATRIX INTO THE C PRODUCT OF A LOWER AND UPPER TRIANGULAR MATRIX (LU DECOMPOSITION). C ROW INTERCHANGES ARE ALWAYS USED (PARTIAL PIVOTING). IN CASE OF C NUMERICAL DIFFICULTIES QMINVZ SWITCHES TO ROW AND COLUMN INTER- C CHANGES (COMPLETE PIVOTING). IDUM( ,1) AND IDUM( ,2) CONTAIN THE C PIVOT ROW AND COLUMN SUBSCRIPTS, RESPECTIVELY, IN THE ORDER CHOSEN C DURING THE ELIMINATION. THE DETERMINANT IS THE PRODUCT OF THE C PIVOT ELEMENTS. IT IS SET =0 WHEN ANY PIVOT ELEMENT IS =0. C----------------------------------------------------------------------- IPIV=1 JPIV=1 PIV=A(1,1) C-----IN DO LOOP 20 SELECT PIVOT=LARGEST OF ALL A(I,J). DO 20I=1,N DO 20J=1,N IF(DABS(A(I,J)).LE.DABS(PIV))GOTO20 IPIV=I JPIV=J PIV=A(I,J) 20 CONTINUE RNORM=DABS(PIV) TOL=1.D-14*RNORM C-----INITIALIZE DETERMINANT. DET=1.D0 IDET=0 H=0.D0 DO 130K=1,N C-----TEST ON SINGULARITY AND LOSS OF SIGNIFICANCE. IF(PIV.EQ.0.D0)GOTO140 IF(IER.NE.0)GOTO25 IF(DABS(PIV).GT.TOL)GOTO25 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 1E INVERSE PRODUCED BY QMINVZ.'/1X,'THE PIVOT ELEMENT AT STEP ',I3, 2' WAS LESS THAN 10**(-15) TIMES THE LARGEST ELEMENT OF THE INPUT M 3ATRIX.'/1X,'THE INPUT MATRIX IS PROBABLY SINGULAR.') CALL END (2) IER=K C-----SAVE ROW AND COLUMN SUBSCRIPTS. 25 IDUM(K,1)=IPIV IDUM(K,2)=JPIV C-----INTERCHANGE ROWS. IF(IPIV.EQ.K)GOTO40 DET=-DET DO 30J=1,N T=A(K,J) A(K,J)=A(IPIV,J) 30 A(IPIV,J)=T C-----INTERCHANGE COLUMNS. 40 IF(JPIV.EQ.K)GOTO60 DET=-DET DO 50I=1,N T=A(I,K) A(I,K)=A(I,JPIV) 50 A(I,JPIV)=T C-----UPDATE DETERMINANT. 60 DET=DET*PIV C-----IE MEASURES TO WHAT POWER OF 10 'DET' WAS. C-----FOR EXAMPLE, LOG 10=1, LOG 100=2, ETC. ZB=DLOG10(DABS(DET)) IE=ZB IF(ZB.LT.0.D0)IE=IE-1 DET=DET*10.D0**(-IE) IDET=IDET+IE K1=K+1 IF(K.EQ.N)GOTO150 DO 70J=K1,N C-----H IS MAX. ELEMENT SEEN IN PIVOT ROWS SO FAR. H=DMAX1(H,DABS(A(K,J))) C-----DIVIDE PIVOT ROW BY PIVOT ELEMENT. 70 A(K,J)=A(K,J)/PIV IPIV=K1 JPIV=K1 PIV=0.D0 C-----IF H >> ORIGINAL MAXIMUM 'A' ELEMENT THEN DO COMPLETE PIVOTING. IF(RNORM+(N-1)*H.GE.8.D0*N*RNORM)GOTO100 C-----ELIMINATION WITH PARTIAL PIVOTING AND SELECTION OF NEXT PIVOT. DO 90I=K1,N DO 80L=K1,N 80 A(I,L)=A(I,L)-A(I,K)*A(K,L) IF(DABS(A(I,K1)).LE.DABS(PIV))GOTO90 IPIV=I PIV=A(I,K1) 90 CONTINUE GOTO 130 C-----ELIMINATION WITH COMPLETE PIVOTING AND SELECTION OF NEXT PIVOT. 100 DO 120I=K1,N S=0.D0 DO 110L=K1,N A(I,L)=A(I,L)-A(I,K)*A(K,L) T=DABS(A(I,L)) IF(T.LE.S)GOTO110 S=T J=L 110 CONTINUE IF(S.EQ.0.D0)GOTO140 IF(DABS(A(I,J)).LE.DABS(PIV))GOTO120 IPIV=I JPIV=J PIV=A(I,J) 120 CONTINUE 130 CONTINUE C-----RETURN FOR SINGULAR CASE. 140 CALL END (1) IF(IER1.EQ.1)WRITE(LUWN,2002) 2002 FORMAT(1X,'NO INVERSE WAS COMPUTED BECAUSE N.LE.1 OR A PIVOT', 1 ' ELEMENT WAS EQUAL TO 0.D0.') CALL END (2) IER=-1 DET=0.D0 IDET=1 IZINCR = IZINCR - 1 RETURN C-----DECOMPOSITION SUCCESSFUL. CALL QMILU TO COMPUTE INVERSE ELEMENTS. 150 CALL QMILU(NRMA,NRMIDU,A,N,IDUM,V) IZINCR = IZINCR - 1 RETURN END SUBROUTINE QMILU(NRMA,NRMIDU,A,N,IDUM,V) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C----------------------------------------------------------------------- C QMILU FINDS THE INVERSE OF THE MATRIX WHOSE LU DECOMPOSITION IS C FURNISHED IN A. C----------------------------------------------------------------------- DIMENSION A(NRMA,1),IDUM(NRMIDU,2),V(1) COMMON/ZINCR/IZINCR CALL QINCR ('QMILU') CALL XDIMCH(NRMA,N,'QMILU',0) CALL XDIMCH(NRMIDU,N,'QMILU',0) DO 80KK=1,N K=N+1-KK K1=K+1 IF(K.EQ.N)GOTO30 DO 20JJ=K1,N J=N+K1-JJ C-----STORE INVERSE ELEMENTS WHICH WERE CALCULATED IN DO LOOP 60. A(J,K1)=V(J) S=0.D0 DO 10L=K1,N 10 S=S-A(K,L)*A(L,J) C-----STORE A**(-1)(K,J) IN V(J). 20 V(J)=S 30 R=A(K,K) IF(K.NE.N)GOTO40 V(K)=1.D0/R GOTO 80 40 DO 60JJ=K1,N J=N+K1-JJ C-----STORE INVERSE ELEMENTS WHICH WERE CALCULATED IN DO LOOP 20. A(K,J)=V(J) S=0.D0 DO 50L=K1,N 50 S=S-A(J,L)*A(L,K) C-----STORE A**(-1)(J,K) IN V(J). 60 V(J)=S/R S=0.D0 DO 70L=K1,N 70 S=S+A(K,L)*A(L,K) C-----STORE DIAGONAL ELEMENTS OF INVERSE IN V(K). V(K)=(1.D0-S)/R 80 CONTINUE NN=2*N C-----STORE FIRST COLUMN OF INVERSE ELEMENTS. DO 90J=1,N 90 A(J,1)=V(J) N1=N-1 C-----PERMUTE COLUMNS. FOR EVERY ROW INTERCHANGE IN THE DECOMPOSITION, C-----MUST INTERCHANGE LIKE COLUMNS OF THE INVERSE. DO 110KK=1,N1 K=N-KK IRK=IDUM(K,1) IF(IRK.EQ.K)GOTO110 DO 100I=1,N R=A(I,IRK) A(I,IRK)=A(I,K) 100 A(I,K)=R 110 CONTINUE C-----PERMUTE ROWS. FOR EVERY COLUMN INTERCHANGE IN THE DECOMPOSITION, C-----MUST INTERCHANGE LIKE ROWS OF THE INVERSE. DO 130KK=1,N K=N+1-KK ICK=IDUM(K,2) IF(ICK.EQ.K)GOTO130 DO 120J=1,N R=A(ICK,J) A(ICK,J)=A(K,J) 120 A(K,J)=R 130 CONTINUE IZINCR = IZINCR - 1 RETURN END SUBROUTINE QDIMCH(N,X,IFLAG) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/B/IO COMMON/ZINCR/IZINCR CALL QINCR ('QDIMCH') CALL XDIMCH(I0,N,X,IFLAG) IZINCR = IZINCR - 1 RETURN END SUBROUTINE QIMTQL(N,D,E,IERR) C*TITLE: EIGENVALUES OF TRIDIAGONAL MATRIX C*PROGRAMMER: EISPACK C*DATE: MAY, 1972 C*LOCATION: ARGONNE NATIONAL LABORATORY C*PURPOSE: DETERMINES EIGENVALUES OF A SYMMETRIC TRIDIAGONAL MATRIX C* *USING THE IMPLICIT QL METHOD C*SYMBOLS: C* *N = ORDER OF THE TRIDIAGONAL MATRIX C* *D = A ONE DIMENSIONAL ARRAY CONTAINING THE DIAGONAL ELEMENTS OF THE C* * * *TRIDIAGONAL MATRIX. ON OUTPUT, IT CONTAINS THE EIGENVALUES IN C* * * *ASCENDING ORDER. C* *E = A ONE DIMENSIONAL ARRAY CONTAINING THE SUBDIAGONAL ELEMENTS OF C* * * *THE TRIDIAGONAL MATRIX IN ITS LAST N-1 ELEMENTS. E(1) IS C* * * *ARBITRARY. QIMTQL DESTROYS THE CONTENTS OF E. C* *IERR = CONTAINS ERROR CODE ON OUTPUT. IF MORE THAN 30 ITERATIONS C* * * * * ARE REQUIRED, QIMTQL TERMINATES IMMEDIATELY WITH IERR SET TO C* * * * * THE INDEX OF THE EIGENVALUES FOR WHICH THE FAILURE OCCURS. C* * * * * IF EVERYTHING IS OK, IERR IS RETURNED AS 0. C*PROCEDURE: THE EIGENVALUES ARE DETERMINED BY THE IMPLICIT QL METHOD. C* *THE ESSENCE OF THIS METHOD IS A PROCESS WHEREBY A SEQUENCE OF C* *SYMMETRIC TRIDIAGONAL MATRICES, UNITARILY SIMILAR TO THE ORIGINAL C* *TRIDIAGONAL MATRIX, IS FORMED WHICH CONVERGES TO A DIAGONAL MATRIX. C* *THE RATE OF CONVERGENCE IS IMPROVED BY IMPLICITLY SHIFTING THE C* *ORIGIN AT EACH ITERATION. BEFORE EACH ITERATION, THE TRIDIAGONAL C* *MATRIX IS CHECKED FOR A POSSIBLE SPLITTING INTO SUBMATRICES. IF A C* *SPLIT OCCURS, ONLY THE UPPERMOST SUBMATRIX IS USED IN THE NEXT C* *ITERATION. THE EIGENVALUES ARE ORDERED IN ASCENDING ORDER AS THEY C* *ARE FOUND. THE ORIGIN SHIFT AT EACH ITERATION IS THE EIGENVALUE OF C* *THE CURRENT UPPERMOST 2 X 2 PRINCIPAL MINOR CLOSER TO THE FIRST C* *DIAGONAL ELEMENT OF THIS MINOR. WHEN THE UPPERMOST 1 X 1 SUBMATRIX C* *FINALLY SPLITS OFF, IT IS TAKEN TO BE AN EIGENVALUE AND EXECUTION C* *PROCEEDS WITH THE REMAINING SUBMATRIX. THIS PROCESS IS CONTINUED C* *UNTIL THE MATRIX HAS SPLIT COMPLETELY INTO SUBMATRICES OF ORDER 1. C* *THE TOLERANCES IN THE SPLITTING TEST ARE A TINY PROPORTION OF THE C* *DIAGONAL ELEMENTS. C*REFERENCES: WILKINSON,J.H., & REINSCH,C. HANDBOOK FOR AUTOMATIC C* *COMPUTATION, VOL. II, LINEAR ALGEBRA, PART 2, SPRINGER-VERLAG, C* *NEW YORK, HEIDELBERG, BERLIN, 1971. AND SAME AS QTRBAK. IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION MACHEP,D(1),E(1) COMMON/ZINCR/IZINCR CALL QINCR ('QIMTQL') C*MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING THE RELATIVE C* *PRECISION OF FLOATING POINT ARITHMETIC. MACHEP=.5D-14 IERR=0 IF(N.EQ.1)GOTO1001 DO 100I=2,N 100 E(I-1)=E(I) E(N)=0.D0 DO 290L=1,N J=0 C*LOOK FOR SMALL SUB-DIAGONAL ELEMENT 105 DO 110M=L,N IF(M.EQ.N)GOTO120 IF(DABS(E(M)).LE.MACHEP*(DABS(D(M))+DABS(D(M+1))))GOTO120 110 CONTINUE 120 P=D(L) IF(M.EQ.L)GOTO215 IF(J.EQ.30)GOTO1000 J=J+1 C*FORM SHIFT G=(D(L+1)-P)/(2.D0*E(L)) R=DSQRT(G*G+1.D0) G=D(M)-P+E(L)/(G+SIGN(R,G)) S=1.D0 C=1.D0 P=0.D0 MML=M-L C*FOR I=M-1 STEP -1 UNTIL L DO DO 200II=1,MML I=M-II F=S*E(I) B=C*E(I) IF(DABS(F).LT.DABS(G))GOTO150 C=G/F R=DSQRT(C*C+1.D0) E(I+1)=F*R S=1.D0/R C=C*S GOTO 160 150 S=F/G R=DSQRT(S*S+1.D0) E(I+1)=G*R C=1.D0/R S=S*C 160 G=D(I+1)-P R=(D(I)-G)*S+2.D0*C*B P=S*R D(I+1)=G+P G=C*R-B 200 CONTINUE D(L)=D(L)-P E(L)=G E(M)=0.D0 GOTO 105 C*ORDER EIGENVALUES 215 IF(L.EQ.1)GOTO250 C*FOR I=L STEP -1 UNTIL 2 DO DO 230II=2,L I=L+2-II IF(P.GE.D(I-1))GOTO270 D(I)=D(I-1) 230 CONTINUE 250 I=1 270 D(I)=P 290 CONTINUE GOTO 1001 C*SET ERROR - NO CONVERGENCE TO AN EIGENVALUE AFTER 30 ITERATIONS 1000 IERR=L 1001 IZINCR = IZINCR - 1 RETURN END SUBROUTINE VARMX(A,F,H,NV,NF) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A(IO,IO),F(IO,IO),H(IO) COMMON/B/IO COMMON/ZINCR/IZINCR CALL QINCR ('VARMX') NRMA = IO NRMF = IO NCMF = IO KRAW = 1 C SET TO DO NORMALIZED VARIMAX CALL XVARMX(NRMA,NRMF,NCMF,A,F,NV,NF,H,KRAW) IZINCR = IZINCR - 1 RETURN END SUBROUTINE XVARMX(NRMA,NRMF,NCMF,A,F,NAT,NF,H,MRAW) C PERFORMS VARIMAX ROTATION USING SUBROUTINE QVARMX C A IS ORIGINAL LOADING MATRIX; INPUT/OUTPUT C F IS ROTATED LOADING MATRIX; OUTPUT C IZ IS THE NUMBER OR ROWS AND COLUMNS OF A AND F IN MAIN C ALSO NUMBER OF ELEMENTS IN H C NAT IS NUMBER OF ATTRIBUTES; INPUT/OUTPUT C NF IS NUMBER OF FACTORS; INPUT/OUTPUT C NI IS THE NUMBER OF ITERATIONS; OUTPUT C H IS VECTOR CONTAINING COMMUNALITIES OF ATTRIBUTES; OUTPUT C MUST HAVE IZ ELEMENTS C MRAW IS AN INTEGER CONCERNING CASE OF VARIMAX; INPUT/OUTPUT C = 0 FOR UNWEIGHTED (RAW) SOLUTION C = 1 FOR NORMALIZED SOLUTION IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION A(NRMF,1),H(1),D(51),F(NRMF,1) COMMON/ZINCR/IZINCR CALL QINCR ('XVARMX') CALL XDIMCH(NRMA,NAT,'VARMX',0) CALL XDIMCH(NRMF,NAT,'VARMX',0) CALL XDIMCH(NCMF,NAT,'VARMX',0) C COPY A INTO F DO 5 I = 1,NAT DO 5 J = 1,NF 5 F(I,J) = A(I,J) KRAW = MRAW C CALL QVARMX CALL QVARMX(NRMA,NCMF,NAT,NF,F,NI,D,H,KRAW) IZINCR = IZINCR - 1 RETURN END SUBROUTINE QVARMX(NRMAV,NCMA,M,K,A,NC,TV,H,KRAW) C PERFORMS VARIMAX ROTATION C REVISION OF VARMX FROM IBM SSP C NCMA IS THE NUMBER OF COLUMNS OF A IN MAIN C M IS NUMBER OF ATTRIBUTES; INPUT/OUTPUT C K IS NUMBER OF FACTORS; INPUT/OUTPUT C A IS FACTOR MATRIX TO BE ROTATED ON INPUT, C ON OUTPUT CONTAINS THE ROTATED FACTOR MATRIX C NC NUMBER OF ITERATIONS; OUTPUT C TV IS A SCRATCH VECTOR WITH 51 ELEMENTS C ON OUTPUT CONTAINS ITERATION VARIANCES OF VARIMAX SOLUTION C H IS VECTOR CONTAINING COMMUNALITIES ON OUTPUT C MUST HAVE M ELEMENTS C KRAW IS AN INTEGER CONCERNING CASE OF VARIMAX SOLUTION C = 0 FOR UNWEIGHTED (RAW) SOLUTION C = 1 FOR NORMALIZED SOLUTION IMPLICIT DOUBLE PRECISION (A-H, O-Z) DOUBLE PRECISION TV(51),H(1),A(1) COMMON/ZINCR/IZINCR CALL QINCR ('QVARMX') ICOUNT = 0 NC = NCMA MAX=NC RAW=KRAW C INITIALIZATION EPS=0.00116D0 TVLT=0.0D0 LL=K-1 NV=1 NC=0 FN=M FFN=FN*FN CONS=0.7071066D0 C CALCULATE ORIGINAL COMMUNALITIES DO 110 I=1,M H(I)=0.0D0 DO 110 J=1,K L=MAX*(J-1)+I 110 H(I)=H(I)+A(L)*A(L) C CALCULATE NORMALIZED FACTOR MATRIX DO 120 I=1,M 115 H(I)= SQRT(H(I)) DO 120 J=1,K L=MAX*(J-1)+I C THIS CARD HAS BEEN MODIFIED TO ALLOW FOR ROTATION WITHOUT WEIGHTING BY COMMUN 120 A(L)=A(L)/(RAW*H(I)+1.D0-RAW) GO TO 132 C CALCULATE VARIANCE FOR FACTOR MATRIX 130 NV=NV+1 TVLT=TV(NV-1) 132 TV(NV)=0.0D0 DO 150 J=1,K AA=0.0D0 BB=0.0D0 LB=MAX*(J-1) DO 140 I=1,M L=LB+I CC=A(L)*A(L) AA=AA+CC 140 BB=BB+CC*CC 150 TV(NV)=TV(NV)+(FN*BB-AA*AA)/FFN IF(NV-51) 160, 430, 430 C PERFORM CONVERGENCE TEST 160 IF((TV(NV)-TVLT)-(1.E-7)) 170, 170, 190 170 NC=NC+1 IF(NC-3) 190, 190, 430 C ROTATION OF TWO FACTORS CONTINUES UP TO C THE STATEMENT 120. 190 DO 420 J=1,LL L1=MAX*(J-1) II=J+1 C CALCULATE NUM AND DEN DO 420 K1=II,K L2=MAX*(K1-1) AA=0.0D0 BB=0.0D0 CC=0.0D0 DD=0.0D0 DO 230 I=1,M L3=L1+I L4=L2+I U=(A(L3)+A(L4))*(A(L3)-A(L4)) T=A(L3)*A(L4) T=T+T CC=CC+(U+T)*(U-T) DD=DD+2.0D0*U*T AA=AA+U 230 BB=BB+T T=DD-2.0D0*AA*BB/FN B=CC-(AA*AA-BB*BB)/FN C COMPARISON OF NUM AND DEN IF(T-B) 280, 240, 320 240 IF((T+B)-EPS) 420, 250, 250 C NUM + DEN IS GREATER THAN OR EQUAL TO THE C TOLERANCE FACTOR 250 COS4T=CONS SIN4T=CONS GO TO 350 C NUM IS LESS THAN DEN 280 TAN4T= ABS(T)/ ABS(B) IF(TAN4T-EPS) 300, 290, 290 290 COS4T=1.0D0/ SQRT(1.0D0+TAN4T*TAN4T) SIN4T=TAN4T*COS4T GO TO 350 300 IF(B) 310, 420, 420 310 SINP=CONS COSP=CONS GO TO 400 C NUM IS GREATER THAN DEN 320 CTN4T= ABS(T/B) IF(CTN4T-EPS) 340, 330, 330 330 SIN4T=1.0D0/ SQRT(1.0D0+CTN4T*CTN4T) COS4T=CTN4T*SIN4T GO TO 350 340 COS4T=0.0D0 SIN4T=1.0D0 C DETERMINE COS THETA AND SIN THETA 350 COS2T= SQRT((1.0D0+COS4T)/2.0D0) SIN2T=SIN4T/(2.0D0*COS2T) 355 COST= SQRT((1.0D0+COS2T)/2.0D0) SINT=SIN2T/(2.0D0*COST) C DETERMINE COS PHI AND SIN PHI IF(B) 370, 370, 360 360 COSP=COST SINP=SINT GO TO 380 370 COSP=CONS*COST+CONS*SINT 375 SINP= ABS(CONS*COST-CONS*SINT) 380 IF(T) 390, 390, 400 390 SINP=-SINP C PERFORM ROTATION 400 DO 410 I=1,M L3=L1+I L4=L2+I AA=A(L3)*COSP+A(L4)*SINP A(L4)=-A(L3)*SINP+A(L4)*COSP 410 A(L3)=AA 420 CONTINUE GO TO 130 C DENORMALIZE VARIMAX LOADINGS 430 DO 440 J=1,K KEEP = MAX*(J-1) CSUM = 0.D0 DO 435 I=1,M L = KEEP + I C THIS CARD HAS BEEN MODIFIED TO ALLOW FOR ROTATION WITHOUT WEIGHTING BY COMMUN A(L)=A(L)*(RAW*H(I)+1.-RAW) 435 CSUM = CSUM + A(L) C CHECK FOR NEGATIVE COLUMN SUMS - IF FOUND, PERFORM 180 DEGREE ROTATE IF (CSUM)436,440,440 436 DO 438 I=1,M L = KEEP + I 438 A(L) = -A(L) 440 CONTINUE C RECONSITUTE COMMUNALITIES NC=NV-1 DO 450 I=1,M 450 H(I)=H(I)*H(I) IZINCR = IZINCR - 1 RETURN END FUNCTION CHIVLT (P,NDF) C VALUE OF CHI-SQUARE WITH NDF DEGREES OF FREEDOM C HAVING PROPORTION P ABOVE C REQUIRES FUNCTION CHIPRA C USES HYPERBOLIC INTERPOLATION C PROGRAM WRITTEN BY: LEDYARD R TUCKER C UNIVERSITY OF ILLINOIS C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DF = NDF SSZ = 2.0D0*DSQRT(2.D0*DF) EPSX = 1.0D-10*SSZ EPSY = 1.0D-10 C INTERVAL SEARCH SECTION YDF = CHIPRA (DF,DF) IF (DABS(YDF - P).LE.EPSY) GOTO 30 IF (P.LT.YDF) GOTO 20 C LOWER SEARCH SECTION XT = DF YT = YDF 10 XB = XT - SSZ IF (XB.LE.EPSX) GOTO 15 YB = CHIPRA (XB,DF) IF (DABS(YB - P).LE.EPSY) GOTO 35 IF (YB.GT.P) GOTO 50 XT = XB YT = YB GOTO 10 15 XB = 0.0D0 YB = 1.0D0 GOTO 50 C UPPER SEARCH SECTION 20 XB = DF YB = YDF 25 XT = XB + SSZ YT = CHIPRA (XT,DF) IF (DABS(YT - P).LE.EPSY) GOTO 40 IF (YT.LE.P) GOTO 50 XB = XT YB = YT GOTO 25 C ESTABLISH CHIVLT FOR SPECIAL CASES 30 CHIVLT = DF GOTO 75 35 CHIVLT = XB GOTO 75 40 CHIVLT = XT GOTO 75 C INTERPOLATION SECTION 50 XM = (XB + XT)/2.0D0 YM = CHIPRA(XM,DF) 55 XTB = XB - XM XTT = XT - XM YTB = YB - YM YTT = YT - YM YTN = P - YM DX = XTT - XTB DY = YTT - YTB T = YTN*(XTT*YTB - XTB*YTT) T = YTT*YTB*DX - T XTN = XTT*XTB*YTN*DY/T XN = XTN + XM IF (XN.LT.(XB + EPSX)) XN = XB + EPSX IF (XN.GT.(XT - EPSX)) XN = XT - EPSX IF (DABS(XN - XM).LE.EPSX) GOTO 70 YN = CHIPRA (XN,DF) IF (DABS(YN - P).LE.EPSY) GOTO 70 IF (XN.GT.XM) GOTO 60 XT = XM YT = YM GOTO 65 60 XB = XM YB = YM 65 XM = XN YM = YN GOTO 55 70 CHIVLT = XN 75 RETURN END FUNCTION QDVNRM (P) C THIS FUNCTION WILL INPUT A RANDOM UNIFORM DEVIATE AND C TRANSFORM IT INTO A NORMAL DEVIATE. C IT WORKS VIA THE INVERSE PROBABILITY FUNCTION. C IN ESSENCE, IF P IS THE CUMULATIVE AREA UNDER A NORMAL CURVE, THEN C QDVNRM WILL BE THE NORMAL DEVIATE WHICH CUTS THE NORMAL C DISTRIBUTION INTO THIS AREA. C DISCREPANCIES ARE IN THE INTERVAL -0.001 TO +0.001 C IN THE RANGE OF P FROM 0.0001 TO 0.9999 C C THE PROGRAM WAS WRITTEN BY LEDYARD TUCKER C UNIVERSITY OF ILLINOIS C IMPLICIT DOUBLE PRECISION (A-H, O-Z) PM = P - 0.5D0 U = DABS(PM) IF (U.GT.0.18D0) GOTO 5 X = 8.24233014D0*U/(3.35216041D0 - U) GOTO 25 5 IF (U.GT.0.32D0) GOTO 10 X = (0.05909112D0 + 2.21627402D0*U)/(1.15930392D0 - U) GOTO 25 10 IF (U.GT.0.405D0) GOTO 15 X = (0.17385088D0 + 0.6838065D0*U)/(0.74895890D0 - U) GOTO 25 15 IF (U.GT.0.48D0) GOTO 20 T1 = -0.20175681D0 + U*(2.36282085D0 - 3.64101574D0*U) X = T1/(0.52553108D0 - U) GOTO 25 20 T1 = 0.5D0 - U IF (T1.LT.1.0D-15) T1 = 1.0D-15 T1 = -DLOG(T1) T2 = DSQRT(T1) X = -1.34882754D0 - 0.04797507D0*T1 + 1.81520099D0*T2 25 QDVNRM = X IF (PM.LT.0.0D0) QDVNRM = -X RETURN END SUBROUTINE XDIMCH(ICHECK,N,X,IFLAG) C TITLE: DIMENSIONALITY CHECK C PROGRAMMER: CARL FINKBEINER C DATE: JUNE, 1974 C LOCATION: UNIVERSITY OF ILLINOIS C PURPOSE: C THIS PROGRAM TESTS A MATRIX DIMENSION SPECIFIED BY THE USER AND C THE MAXIMUM ALLOWABLE DIMENSION AND OUTPUTS AN APPROPRIATE ERROR C MESSAGE C SYMBOLS: C N = USER SPECIFIED DIMENSION C X = A DOUBLE PRECISION VARIABLE CONTAINING THE SUBROUTINE NAME C (UP TO 6 CHARACTERS) C WHICH IS CALLING QDIMCH C IFLAG = 0 IF DESIRE ERROR MESSAGE IN FORMAT 1 C 1 IF DESIRE ERROR MESSAGE IN FORMAT 2 C IPGM = PROGRAM # OF CALLING SUBROUTINE C ICHECK = MAXIMUM ALLOWABLE DIMENSIONALITY C A = LABELED COMMON AREAS CONTAINING IPGM C PROCEDURE: N IS TESTED AGAINST ICHECK AND IF N .LE. ICHECK, A NORMAL RETURN OC C IF ICHECK .LE. 0, A MESSAGE IS PRINTED AND PROCESSING IS TERMINATED. C IF N .GT. ICHECK, ONE OF THE ERROR MESSAGES IN FORMAT STATEMENTS 1 OR 2 C PRINTED OUT AND EXECUTION IS TERMINATED. NOTE THAT THIS PROGRAM CAN C EASILY BE MODIFIED TO PRINT OUT DIFFERENT ERROR MESSAGES BY PUTTING C ADDITIONAL TESTS ON IFLAG AND INSERTING APPROPRIATE WRITE STATEMENTS IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON/A/IPGM,ISUBRU COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG COMMON/ZINCR/IZINCR CALL QINCR ('XDIMCH') C TEST N AGAINST 0 IF(N.LE.0)GOTO 3 IF(ICHECK.LE.0)GOTO3 C TEST N AGAINST ICHECK IF(N.LE.ICHECK)IZINCR = IZINCR - 1 IF(N.LE.ICHECK)RETURN C IF TEST FAILS, CHOOSE APPROPRIATE ERROR MESSAGE IF(IFLAG.EQ.0)GOTO5 C THIS ERROR MESSAGE IS, SO FAR, USED ONLY BY THE VERTICAL AND HORIZONTA C CONCATENATION SUBROUTINES CALL END (1) WRITE(LUWN,2) X 2 FORMAT('0* * * * * T I L T !!! THE VALUES YOU SPECIFIED FOR THE 1DIMENSIONS OF THE TWO MATRICES TO BE CONCATENATED BY ',A6,' WOULD' 2//25X,'HAVE PRODUCED A SUPER MATRIX OF ORDER GREATER THAN THE MAXI 3MUM SPECIFIED FOR YOUR SUPERMATRIX.') CALL END (2) CALL END (3) C THIS IS THE USUAL ERROR MESSAGE AND IS USED BY MOST OF THE SUBROUTINES 5 CALL END (1) WRITE(LUWN,1) ICHECK, X 1 FORMAT('0 YOU HAVE SPECIFIED A VALUE GREATER THAN',I3,' FOR AT LEA 2ST ONE OF THE DIMENSIONS OF A MATRIX, IN PROGRAM ',A6,'.') CALL END (2) CALL END (3) C THIS ERROR MESSAGE IS FOR NEGATIVE DIMENSIONS 3 CALL END (1) WRITE(LUWN,4)N,ICHECK,X 4 FORMAT('0YOU HAVE SPECIFIED A VALUE EQUAL OR LESS THAN 0 FOR AT LE 1AST ONE OF THE DIMENSIONS OF A MATRIX.', 2//24X,' YOU WANTED TO USE A DIMENSIONALITY OF',I5, 3' WHEN YOU SAID THAT THE MATRIX SIZE TO BE ',I5,//11X, 4' WHICH WAS A CALL TO THE ',A6,' PROGRAM') CALL END (2) CALL END (3) END SUBROUTINE QINCR (PRGNAM) C TITLE: PROGRAM COUNTER C PROGRAMMER: CARL FINKBEINER C DATE: JUNE, 1974 C LOCATION: UNIVERSITY OF ILLINOIS C PURPOSE: THIS SUBROUTINE IS CALLED BY ALL OF THE FORMAL C SUBROUTINES IN ORDER TO KEEP TRACK OF THE NUMBER OF SUBROUTINES C CALL BY THE USER. THIS NUMBER MAY BE PRINTED OUT WHEN AN ERROR C CONDITION HAS OCCURED. THE FIRST TIME QINCR IS CALLED, A FORMAL C HEADING IS PRINTED OUT, AND DATA INPUT AND OUTPUT FILES ARE C OPENED TO DEFAULTS OF 'FRMLIN.DAT' AND 'FRMLOT.DAT', RESPECTIVELY. C SYMBOLS: C A = LABELED COMMON AREA CONTAINING IPGM C IPGM = THE CURRENT NUMBER OF FORMAL CALL STATEMENTS MADE FROM THE C MAIN PROGRAM. C ISUB = THE NUMBER OF SUBROUTINES USED INDIRECTLY SINCE C THE LAST CALL IN MAIN, WILL BE ZERO AT CALLS FROM C MAIN. C IZINCR = WILL BE 1 AT ALL FORMAL CALLS MADE FROM THE C MAIN PROGRAM. 2 WILL BE AT ALL CALLS MADE ONE C SUBROUTINE FROM THE MAIN, ETC. C PRGNAM = THE NAME OF THE CURRENT SUBROUTINE C FRSTPG = ORIGINAL CALL TO FORMAL BY USER (I.E., WHEN IZINCR WAS 1) C PROCEDURE: QINCR ADDS 1 TO IPGM, TESTS IPGM TO SEE IF THIS IS THE FIRST C QINCR AND IF SO, IT PRINTS OUT A HEADING. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON/A/IPGM,ISUB COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XINPUT(5),ZOUTPT(5),FRSTPG COMMON/ZINCR/IZINCR IZINCR = IZINCR + 1 IF (IZINCR.EQ.1) FRSTPG = PRGNAM IF (IZINCR.GE.1.AND.IPGM.GE.0) GOTO 1 CALL END (1) WRITE (LUWN,4) 4 FORMAT (' STOP. SOMETHING IS WRONG WITH SUBROUTINE "QINCR"', 2 ' SEE DR. FLEISHMAN') CALL END (2) CALL END (3) 1 IF (IZINCR.EQ.1.AND.IPGM.LT.1) GOTO 3 IF (IZINCR.EQ.1) GOTO 2 ISUB = ISUB + 1 RETURN C INCREASE IPGM BY 1 2 IPGM=IPGM+1 ISUB = 0 C IF THIS IS THE FIRST TIME QINCR HAS BEEN CALLED BY A SUBROUTINE, C THE HEADING WRITTEN OUT. 3 IF (IPGM.LE.1) CALL IOPUT ('FRMLIN.DAT 2 ','FRMLOT.DAT ' 3 ,3) 5 RETURN END SUBROUTINE IOPUT (XINPUT, OUTPUT, IIFLAG) C TITLE: Input Output controller C Programmer: Allen I. Fleishman C Date: August, 1979 C Location: LEDERLE LABS C Purpose: This subroutine will change the name of the input C and/or output file to 'XINPUT' and 'OUTPUT', respectively. C IFLAG will indicate the type of change: C 1 input file - 'XINPUT' only will be opened; a close C to the old input file will be done if necessary. C 2 output file - 'OUTPUT' only will be opened; a message C to this newly opened file will be generated; if any C FORMAL calls were made previous to the first use of C this program, then the previous output file (default C is 'FRMLOT.DAT') will be closed and a message will C be made where future output can be found; the old C file will also be closed. C 3 both input and output file will be opened as by options C 1 and 2 above. C Symbols: C 'XINPUT' = name of file containing input data (in the form C of filename.ext - ext will default to DAT) C 'OUTPUT' = name of file which will contain the output data (in the C form of filename.ext - ext will default to DAT) C IFLAG = flag denoting whether input and/or output will be C done. C A = labeled common area containing IPGM C C = labeled common area containing DSNR (name of input device and NMAT C IPGM = program counter to be incremented by this subroutine C Procedure: Adds 1 to IPGM as in QINCR, tests to see if this is C the the first usage, if so defaults will be C 'FRMLIN.DAT' and 'FRMLOT.DAT'. Files will be opened with the C appropriate header. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DOUBLE PRECISION XINPUT(5),OUTPUT(5),OLDIN(5),OLDOUT(5) COMMON/A/IPGM,ISUB COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XIN(5),OUT(5),FRSTPG COMMON/ZINCR/IZINCR COMMON/KMISNR/KMISNR DATA (XIN(I),I=1,5)/' ',' ',' ', 2' ',' '/ DATA (OUT(I),I=1,5)/' ',' ',' ', 2' ',' '/ DATA (OLDIN(I),I=1,5)/' ',' ',' ', 2' ',' '/ DATA (OLDOUT(I),I=1,5)/' ',' ',' ', 2' ',' '/ DATA IPGM/0/,ISUB/0/,LURN/2/,LUWN/3/,LUEX/25/,NMAT/0/,IZINCR/0/ DATA KMISNR/0/,LUPN/1/ C C THIS BLOCK OF DATA WILL INITIALIZE: C THE PROGRAM COUNTER - IPGM = 0 C THE SUBROUTINE COUNTER - ISUB = 0 C THE (LOGICAL) READ UNIT - LURN = 2 OR THE CARD READER C THE (LOGICAL) WRITE UNIT - LUWN = 3 OR THE LINE PRINTER C THE (LOGICAL) EXTRA UNIT - LUEX = 25 (A TEMPORARY UNIT) C THE (LOGICAL) PUNCH UNIT - LUPN = 1 OR THE DISK AREA C THE MATRIX BEING READ - NMAT = 0 C THE SUBROUTINE 'DEPTH' - IZINCR = 0 C THE MISNR SWITCH - KMISNR = 0 (SEE MISNR FOR DETAILS) C IFLAG=IIFLAG IF(IZINCR.LT.0.AND.IPGM.NE.0.OR.IZINCR.GT.1) CALL END (1) IF(IZINCR.LT.0.AND.IPGM.NE.0.OR.IZINCR.GT.1) WRITE (LUWN,1003) 1003 FORMAT ('1STOP. SOMETHING IS WRONG WITH SUBROUTINE "IO"', 2 'SEE DR. FLEISHMAN') IF(IZINCR.LT.0.AND.IPGM.NE.0.OR.IZINCR.GT.1) CALL END (2) IF(IZINCR.LT.0.AND.IPGM.NE.0.OR.IZINCR.GT.1) CALL END (3) IF (IPGM.GE.1) GOTO 30 IF (IFLAG.EQ.1) OUTPUT(1) = 'FRMLOT.DAT' IF (IFLAG.EQ.1) OLDOUT(1) = 'FRMLOT.DAT' IF (IFLAG.EQ.2) XINPUT(1) = 'FRMLIN.DAT' IF (IFLAG.EQ.2) OLDIN(1) = 'FRMLIN.DAT' IF (IFLAG.EQ.1.OR.IFLAG.EQ.2) IFLAG = 3 30 IF (IFLAG.GE.1.AND.IFLAG.LE.3) GOTO 99 WRITE (LUWN, 98) IFLAG 98 FORMAT (' ',13X, '* * * * * T I L T * *' 2' * * *'/'0',13X, 'YOU DID NOT GIVE THE INPUT/OUTPUT' 3 ' PROGRAM "IO" A VALID FLAG.'/'0', 13X, 'THE FLAG SHOULD BE' 4 ' IN THE RANGE OF 1 TO 3. YOU GAVE IT A VALUE OF ',I9) CALL END (2) CALL END (3) 99 IF (IFLAG.EQ.2) GOTO 2 IF (IPGM.LT.1) GOTO 4 IF (XINPUT(1).EQ.OLDIN(1)) GOTO 54 CLOSE (UNIT=LURN,DIALOG=OLDIN) 4 OPEN (UNIT=LURN,DIALOG=XINPUT,ACCESS='SEQIN') CALL SET (OLDIN,XINPUT) 54 IF (IFLAG.EQ.1) GOTO 999 2 IF (IPGM.LT.1) GOTO 5 IF (OUTPUT(1).EQ.OLDOUT(1)) GOTO 999 29 CONTINUE WRITE (LUWN,3) (OUTPUT(I),I=1,5) 3 FORMAT ('1', 13X, '* * * * THIS FILE IS BEING CLOSED', 2 ' * * * *'/'0',15X, 'FOR FUTURE OUTPUT SEE FILE ',5A10) CLOSE (UNIT=LUWN,DIALOG=OLDOUT) 5 OPEN (UNIT=LUWN,DIALOG=OUTPUT,ACCESS='SEQOUT') IF (IPGM.LT.1) GOTO 18 WRITE (LUWN, 19) (OLDOUT(I), I=1,5) 19 FORMAT ('1', 13X, '* * * * THIS FILE IS BEING OPENED', 2 ' * * * *'/'0',15X, 'FOR POSSIBLE PAST OUTPUT SEE FILE ',5A10) 18 CALL SET (OLDOUT,OUTPUT) CALL TIME (Y,Z) CALL DATE (X) WRITE (LUWN,1) X, Y 1 FORMAT('0',58X,'* F O R M A L *'/50X,'(FORTRAN MATRIX ALGEBRA LIBR 1ARY)'/' ',51X, 'DATE: ',A9,2X,'TIME: ',A5/'0') C OPEN (UNIT=LUPN,DIALOG='PUNCH',ACCESS='APPEND') 999 IPGM=IPGM+1 ISUB = 0 IF (IFLAG.EQ.1.OR.IFLAG.EQ.3) CALL SET (XIN,XINPUT) IF (IFLAG.EQ.2.OR.IFLAG.EQ.3) CALL SET (OUT,OUTPUT) RETURN END SUBROUTINE SET (X,Y) DOUBLE PRECISION X(5), Y(5) DO 1 I = 1,5 1 X(I) = Y(I) RETURN END SUBROUTINE END (I) IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON/A/IPGM,ISUBRU COMMON/C/LURN,LUWN,LUEX,LUPN,NMAT,XIN(5),OUT(5),FRSTPG COMMON/ZINCR/IZINCR GOTO (10, 20, 30), I 10 WRITE (LUWN,1) 1 FORMAT (' ',30X, '* * * * * T I L T * *', 2' * * *'//) RETURN 20 WRITE(LUWN,4) FRSTPG,IPGM,ISUBRU,IZINCR 4 FORMAT(//20X,' FATAL ERROR IN THE ',A6,' PROGRAM. THIS WAS THE ', 2 I4,'-TH CALL MADE BY YOU TO THE FORMAL LIBRARY.'/20X, 3' ACTUALLY ',I4,' MORE CALLS WERE MADE TO THE FORMAL PACKAGE (INTE 4RNALLY), CURRENT DEPTH: ',I4/) RETURN 30 WRITE (LUWN, 3) 3 FORMAT (//11X, 2' * * * * * P R O C E S S I N G M U S T T E R M I N A T E 3 P R E M A T U R E L Y * * * * *') STOP END