Trailing-Edge
-
PDP-10 Archives
-
decuslib10-12
-
43,50550/forml2.for
There are no other files named forml2.for in the archive.
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