home *** CD-ROM | disk | FTP | other *** search
/ High Voltage Shareware / high1.zip / high1 / DIR9 / ACMALG01.ZIP / ACM633.FOR < prev    next >
Text File  |  1991-02-18  |  109KB  |  3,541 lines
  1.  
  2. C     ALGORITHM 633 COLLECTED ALGORITHMS FROM ACM.
  3. C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.11, NO. 2,
  4. C     JUN., 1985, PP. 170-182.
  5.       INTEGER I, IPRINT, J, KRES, LDL, LDX, LMAX, LOCIW, LOCW,
  6.      *     N, P, P2
  7.       INTEGER ITITLE(72), IWORK(50), LVAR(6,10)
  8.       DOUBLE PRECISION PCERR, WORK(500), X(20,6), XC(20,6)
  9.       LDX = 20
  10.       LOCW = 500
  11.       LDL = 6
  12.       LOCIW = 50
  13.       READ (5,99999) (ITITLE(I),I=1,72)
  14.       WRITE (6,99996) (ITITLE(I),I=1,72)
  15.       DO 60 K=1,2
  16.            READ (5,99998) N, P, JOBAB, PCERR, IPRINT, KRES
  17.            IF (K.GE.2) GO TO 20
  18.            DO 10 I=1,N
  19.                 READ (5,99997) (X(I,J),J=1,P)
  20.    10      CONTINUE
  21.    20      CONTINUE
  22.            DO 40 I=1,N
  23.                 DO 30 J=1,P
  24.                      XC(I,J) = X(I,J)
  25.    30           CONTINUE
  26.    40      CONTINUE
  27.            LMAX = 1
  28.            IF (K.EQ.1) GO TO 50
  29.            P2 = 2
  30.            LMAX = 6
  31.    50      CALL LDA(XC, LDX, N, P, JOBAB, PCERR, P2, LVAR, LDL,
  32.      *          LMAX, KRES, IPRINT, WORK, LOCW, IWORK, LOCIW)
  33.    60 CONTINUE
  34.       STOP
  35. 99999 FORMAT (72A1)
  36. 99998 FORMAT (3I3, F3.0, 10I3)
  37. 99997 FORMAT (F5.1, F7.0, F5.0, F5.0, F7.0, F5.0, F6.0)
  38. 99996 FORMAT (1H1, 72A1)
  39.       END
  40. C
  41. C --------------------------------------------------------------
  42. C
  43.       SUBROUTINE LDA(X, LDX, N, P, JOBAB, PCERR, P2, LVAR, LDL,
  44.      *     LMAX, KRES, IPRINT, WORK, LOCW, IWORK, LOCIW)
  45. C ***
  46. C
  47. C     LDA PERFORMS A LINEAR DEPENDENCY ANALYSIS, WHICH ATTEMPTS TO
  48. C        IDENTIFY THE NUMBER OF LINEAR DEPENDENCIES AND THE BEST SUBSETS
  49. C        OF ESTIMATED AND PREDICTOR VARIABLES.
  50. C
  51. C     DESCRIPTION OF ARGUMENTS:
  52. C
  53. C        X       REAL(LDX,P).        LDX .GE. MAX(N,P).
  54. C                MATRIX CONTAINING THE N OBSERVATIONS IN ROWS WITH THE
  55. C                P VARIABLES IN COLUMNS. X IS DESTROYED BY LDA.
  56. C
  57. C        LDX     INTEGER.
  58. C                LDX IS THE LEADING DIMENSION OF THE ARRAY X. LDX MUST
  59. C                BE GREATER THAN OR EQUAL TO THE MAXIMUM OF N AND P.
  60. C
  61. C        N       INTEGER.            N .GT. 1
  62. C                N IS THE NUMBER OF ROWS (OBSERVATIONS) IN THE MATRIX X.
  63. C
  64. C        P       INTEGER.            P .GT. 1
  65. C                P IS THE NUMBER OF COLUMNS (VARIABLES) IN THE MATRIX X.
  66. C
  67. C        JOBAB   INTEGER.
  68. C                JOBAB CONTROLS THE TYPE OF COMPUTATIONS PERFORMED.
  69. C                IT HAS THE DECIMAL EXPANSION  AB  WITH THE FOLLOWING
  70. C                MEANING:
  71. C
  72. C                   A  =  1    DETERMINE NUMBER OF DEPENDENCIES FROM
  73. C                              PCERR.
  74. C                      =  2    NUMBER OF DEPENDENCIES IS SPECIFIED BY
  75. C                              INPUT PARAMETER P2.
  76. C                      =  3    NUMBER OF DEPENDENCIES IS SPECIFIED BY
  77. C                              P2 AND SPECIFIED SETS OF ESTIMATED
  78. C                              VARAIABLES GIVEN IN LVAR.
  79. C
  80. C                   B  =  0    COMPUTE ONLY BENCHMARK SOLUTION (USUALLY
  81. C                              CLOSE TO OPTIMAL SOLUTION)
  82. C                      =  1    TRY ALL POSSIBLE SUBSETS OF VARIABLES
  83. C                              TO FIND SETS WITH SMALLEST "REDUCED
  84. C                              SPACED" RESIDUALS.
  85. C                FOR EXAMPLE: IF JOBAB = 10, LDA WILL USE PCERR TO
  86. C                   DETERMINE THE NUMBER OF DEPENDENCIES AND ONLY
  87. C                   COMPUTE THE BENCHMARK SOLUTION.
  88. C
  89. C        PCERR   REAL.     PCERR .GE. 0  AND  PCERR .LE. 100.
  90. C                WHEN A = 1 IN JOBAB, PCERR IS THE MAXIMUM PERCENTAGE
  91. C                   OF UNEXPLAINED ERROR VARIATION WHICH IS ACCEPTABLE,
  92. C                   ASSUMING A DIAGONAL STRUCTURE IN THE ERROR MATRIX.
  93. C                IF A = 2 OR 3 IN JOBAB, PCERR IS NOT REFERENCED.
  94. C
  95. C        P2      INTEGER.      P2 .GT. 0  AND  P2 .LT. P.
  96. C                IF A = 2 OR 3 IN JOBAB, P2 IS THE NUMBER OF DEPENDENT
  97. C                   VARIABLES FOR THE ANALYSIS.
  98. C                IF A = 1 IN JOBAB, P2 MUST APPEAR AS AN ARGUMENT BUT
  99. C                   MAY BE UNSPECIFIED. THE COMPUTED VALUE OF P2 IS
  100. C                   RETURNED.
  101. C
  102. C        LVAR    INTEGER(LDL,LMAX).      LDL .GE. P2.
  103. C                IF A = 1 OR 2 IN JOBAB, LVAR MUST APPEAR AS AN ARGUMENT
  104. C                   BUT MAY BE UNSPECIFIED.
  105. C                IF A = 3 IN JOBAB, LVAR CONTAINS LMAX SETS OF ESTIMATED
  106. C                   VARIABLES TO ANALYZE. EACH SET CONTAINS P2
  107. C                   VARIABLES.
  108. C
  109. C        LDL     INTEGER.
  110. C                LDL IS THE LEADING DIMENSION OF THE ARRAY LVAR. LDL
  111. C                MUST BE GREATER THAN OR EQUAL TO THE INPUT VALUE OF P2
  112. C                (A=2,3) OR ITS EXPECTED COMPUTED VALUE (A=1).
  113. C
  114. C        LMAX    INTEGER.
  115. C                FOR A = 1 OR 2 IN JOBAB:
  116. C                   IF B = 1 IN JOBAB, LMAX IS THE MAXIMUM NUMBER OF
  117. C                      COMPETING DEPENDENCIES TO OUTPUT FOR USER
  118. C                      INSPECTION.
  119. C                   IF B = 0 IN JOBAB, LMAX MUST BE SET TO 1.
  120. C                FOR A = 3 IN JOBAB:
  121. C                   LMAX IS THE NUMBER OF SETS OF ESTIMATED VARIABLES
  122. C                   IN LVAR TO ANALYZE.
  123. C
  124. C        KRES    INTEGER.
  125. C                IF B = 1 IN JOBAB, A SET OF ESTIMATED VARIABLES WILL
  126. C                   NOT BE ANALYZED FURTHER IF THE "REDUCED SPACE"
  127. C                   RESIDUAL IS 10**KRES TIMES LARGER THAN THE SMALLEST
  128. C                   CURRENT "REDUCED SPACE" RESIDUAL COMPUTED (BENCHMARK
  129. C                   SET EVALUATED FIRST). IF KRES IS LESS THAN ZERO, THE
  130. C                   DEFAULT VALUE OF 2 IS USED.
  131. C                IF B = 0 IN JOBAB, KRES IS NOT REFERENCED.
  132. C
  133. C        IPRINT  INTEGER.     IPRINT .GE. 1  AND  IPRINT .LE. 4.
  134. C                IPRINT CONTROLS THE AMOUNT OF OUTPUT. FOR MOST ANALYSIS
  135. C                IPRINT = 1 WILL PROBABLY BE SATISFACTORY. AS IPRINT
  136. C                INCREASES, THE AMOUNT OF OUTPUT INCREASES.
  137. C
  138. C        WORK    REAL(LOCW).
  139. C                SCRATCH STORAGE USED BY LDA. USER NEED NOT SPECIFY
  140. C                ANY ELEMENTS IN WORK.
  141. C
  142. C        LOCW    INTEGER.
  143. C                THE NUMBER OF ELEMENTS IN WORK. LOCW MUST BE LARGER
  144. C                THAN    N*P + 3*P**2 + 4*P + N + 5*LMAX.
  145. C
  146. C        IWORK   INTEGER(LOCIW).
  147. C                SCRATCH STORAGE USED BY LDA. USER NEED NOT SPECIFY
  148. C                ANY ELEMENTS IN IWORK.
  149. C
  150. C        LOCIW   INTEGER.
  151. C                THE NUMBER OF ELEMENTS IN IWORK. LOCIW MUST BE LARGER
  152. C                THAN    4*P.
  153. C
  154. C     THIS VERSION OF LDA IS DATED 10/1/84.
  155. C
  156. C     REFERENCES:
  157. C
  158. C        V. E. KANE, R. C. WARD, AND G. J. DAVIS, "ASSESSMENT OF LINEAR
  159. C           DEPENDENCIES IN MULTIVARIATE DATA," SIAM J. SCI. STAT.
  160. C           COMPUT., TO APPEAR
  161. C
  162. C        R. C. WARD, G. J. DAVIS, AND V. E. KANE, "AN ALGORITHM FOR
  163. C           LINEAR DEPENDENCY ANALYSIS OF MULTIVARIATE DATA,"
  164. C           ACM TOMS, SUBMITTED.
  165. C
  166. C     THE LDA SUPPLIED CODE CONSISTS OF THE FOLLOWING SUBROUTINES:
  167. C
  168. C        BETA, CHKEV, DEPEND, LDA, PICVAR, PRINTM, RESFLY, SCPROD.
  169. C
  170. C     THE LDA SUPPLIED CODE CALLS THE FOLLOWING ADDITIONAL SOFTWARE:
  171. C
  172. C        LINPACK:   DPODI, DPOFA, DPOSL, DSVDC
  173. C        EISPACK:   RS
  174. C        BLAS:      DAXPY, DCOPY, DDOT, DNRM2, DSCAL, DROT, DROTG, DSWAP
  175. C        FORTRAN:   DABS,         DLOG, DMAX1, DSIGN, DSQRT, IFIX, MIN0,
  176. C                   SNGL
  177. C
  178. C ***
  179.       INTEGER I, I1, I10, I2, I3, I4, I5, I6, I7, I8, I9, ID,
  180.      *     IERR, IPRINT, J, J1, J2, J3, J4, JOB, JOBAB, KALL,
  181.      *     KRES, LDL, LDX, LMAX, LOCIW, LOCW, MAX0, N, P, P2
  182.       INTEGER IWORK(LOCIW), LVAR(LDL,LMAX), IDUM(1,1)
  183.       DOUBLE PRECISION PCERR, WORK(LOCW), X(LDX,P), DUM(1,1)
  184.       IERR = 0
  185.       ID = 1
  186.       CALL PRINTM(1, DUM, ID, ID, ID, IDUM, ID, ID, ID)
  187.       WRITE (6,99998)
  188.       WRITE (6,99997) N, P, LDX, JOBAB, IPRINT
  189.       IF (LDX.GE.MAX0(N,P)) GO TO 10
  190.       WRITE (6,99996)
  191.       IERR = 1
  192.    10 IF (MIN0(N,P).GE.2) GO TO 20
  193.       WRITE (6,99995)
  194.       IERR = 1
  195.    20 IF (IPRINT.GE.1 .AND. IPRINT.LE.4) GO TO 30
  196.       WRITE (6,99994)
  197.       IERR = 1
  198.    30 JOB = JOBAB/10
  199.       KALL = JOBAB - JOB*10
  200.       IF (JOB.EQ.3) KALL = 1
  201.       IF (JOB.GE.1 .AND. JOB.LE.3) GO TO 40
  202.       WRITE (6,99993)
  203.       GO TO 150
  204.    40 IF (KALL.EQ.0 .OR. KALL.EQ.1) GO TO 50
  205.       WRITE (6,99992)
  206.       IERR = 1
  207.    50 IF (JOB.NE.1) GO TO 60
  208.       WRITE (6,99991) PCERR
  209.       IF (PCERR.GE.0.D0.AND.PCERR.LT.100.D0) GO TO 60
  210.       WRITE (6,99990)
  211.       IERR = 1
  212.    60 IF (JOB.NE.2) GO TO 70
  213.       WRITE (6,99989) P2
  214.       IF (P2.GT.0 .AND. P2.LT.P) GO TO 70
  215.       WRITE (6,99988)
  216.       IERR = 1
  217.    70 IF (JOB.NE.3) GO TO 120
  218.       WRITE (6,99987) P2, LMAX
  219.       DO 80 I=1,P2
  220.            WRITE (6,99986) (LVAR(I,J),J=1,LMAX)
  221.    80 CONTINUE
  222.       IF (P2.GT.0 .AND. P2.LT.P) GO TO 90
  223.       WRITE (6,99988)
  224.       IERR = 1
  225.    90 DO 110 I=1,P2
  226.            DO 100 J=1,LMAX
  227.                 K = LVAR(I,J)
  228.                 IF (K.GE.1 .AND. K.LE.P) GO TO 100
  229.                 WRITE (6,99999) I, J
  230.                 IERR = 1
  231.   100      CONTINUE
  232.   110 CONTINUE
  233.       GO TO 130
  234.   120 IF (KALL.EQ.1) WRITE (6,99985) KRES, LMAX
  235.       IF (KALL.NE.0) GO TO 130
  236.       IF (LMAX.EQ.1) GO TO 130
  237.       WRITE (6,99984)
  238.       IERR = 1
  239.   130 IF (KALL.EQ.1 .AND. KRES.LT.0) KRES = 2
  240.       I1 = 1 + N*P
  241.       I2 = I1 + P
  242.       I3 = I2 + P*P
  243.       I4 = I3 + P*P
  244.       I5 = I4 + P
  245.       I6 = I5 + P*P
  246.       I7 = I6 + N
  247.       I8 = I7 + P
  248.       I9 = I8 + LMAX*5
  249.       I10 = I9 + P - 1
  250.       J1 = P + 1
  251.       J2 = J1 + P
  252.       J3 = J2 + P
  253.       J4 = J3 + P - 1
  254.       WRITE (6,99983) LOCW, I10
  255.       WRITE (6,99982) LOCIW, J4
  256.       IF (I10.LE.LOCW .AND. J4.LE.LOCIW) GO TO 140
  257.       WRITE (6,99981)
  258.       IERR = 1
  259.   140 IF (IERR.NE.0) GO TO 150
  260.       WRITE (6,99980)
  261.       CALL DEPEND(X, LDX, N, P, JOB, KALL, PCERR, P2, LVAR,
  262.      *     LDL, LMAX, KRES, WORK(1), WORK(I1), WORK(I2),
  263.      *     WORK(I3), WORK(I4), WORK(I5), WORK(I6), WORK(I7),
  264.      *     WORK(I8), WORK(I9), IWORK(1), IWORK(J1), IWORK(J2),
  265.      *     IWORK(J3), IPRINT)
  266.   150 RETURN
  267. 99999 FORMAT (/6H LVAR(, I3, 1H,, I3, 22H) LESS THAN 1 OR GREAT,
  268.      *     9HER THAN P//36H **** PROGRAM ABORTING AFTER CHECKIN,
  269.      *     22HG OTHER ARGUMENTS ****)
  270. 99998 FORMAT (1X, 16HINPUT PARAMETERS/)
  271. 99997 FORMAT (4H N =, I4/4H P =, I3/6H LDX =, I4/8H JOBAB =,
  272.      *     I3/9H IPRINT =, I3)
  273. 99996 FORMAT (/45H LDX NOT GREATER THAN OR EQUAL TO MAXIMUM OF ,
  274.      *     7HN AND P//38H **** PROGRAM ABORTING AFTER CHECKING ,
  275.      *     20HOTHER ARGUMENTS ****)
  276. 99995 FORMAT (/26H N OR P NOT GREATER THAN 1//14H **** PROGRAM ,
  277.      *     44HABORTING AFTER CHECKING OTHER ARGUMENTS ****)
  278. 99994 FORMAT (/27H IPRINT NOT BETWEEN 1 AND 4//13H **** PROGRAM,
  279.      *     45H ABORTING AFTER CHECKING OTHER ARGUMENTS ****)
  280. 99993 FORMAT (/34H A IN JOBAB IS NOT BETWEEN 1 AND 3//7H **** P,
  281.      *     20HROGRAM ABORTING ****)
  282. 99992 FORMAT (/25H B IN JOBAB IS NOT 0 OR 1//15H **** PROGRAM A,
  283.      *     43HBORTING AFTER CHECKING OTHER ARGUMENTS ****)
  284. 99991 FORMAT (8H PCERR =, F7.2)
  285. 99990 FORMAT (/30H PCERR NOT BETWEEN 0. AND 100.//10H **** PROG,
  286.      *     48HRAM ABORTING AFTER CHECKING OTHER ARGUMENTS ****)
  287. 99989 FORMAT (5H P2 =, I3)
  288. 99988 FORMAT (/23H P2 NOT BETWEEN 0 AND P//17H **** PROGRAM ABO,
  289.      *     41HRTING AFTER CHECKING OTHER ARGUMENTS ****)
  290. 99987 FORMAT (5H P2 =, I3/7H LMAX =, I3/16H LVAR(P2,LMAX) =)
  291. 99986 FORMAT (1X, 40I3)
  292. 99985 FORMAT (7H KRES =, I3/7H LMAX =, I3)
  293. 99984 FORMAT (/43H B IN JOBAB IS ZERO, LMAX MUST BE SET TO 1.//
  294.      *     49H **** PROGRAM ABORTING AFTER CHECKING OTHER ARGUM,
  295.      *     9HENTS ****)
  296. 99983 FORMAT (7H LOCW =, I7, 10X, 26H(ACTUAL WORK SPACE NEEDED ,
  297.      *     1H=, I7, 1H))
  298. 99982 FORMAT (8H LOCIW =, I6, 10X, 25H(ACTUAL IWORK SPACE NEEDE,
  299.      *     3HD =, I6, 1H))
  300. 99981 FORMAT (/43H USER DID NOT ALLOCATE ENOUGH WORKING SPACE//
  301.      *     27H **** PROGRAM ABORTING ****)
  302. 99980 FORMAT (////)
  303.       END
  304. C
  305. C --------------------------------------------------------------
  306. C
  307.       SUBROUTINE DEPEND(X, LDX, N, P, JOB, KALL, PCERR, P2,
  308.      *     LVAR, LDL, LMAX, KRES, XSAVE, W, V, CORR, XMU, R,
  309.      *     WK, WK2, OUTPUT, COVD, IVAR, INT, DVAR, PVAR, IPRINT)
  310. C ***
  311. C     SUBROUTINE DEPEND DRIVES THE COMPUTATIONS IN THE LINEAR
  312. C        DEPENDENCY ANALYSIS CALLING ON VARIOUS SUBROUTINES TO
  313. C        COMPUTE NEEDED QUANTITIES.
  314. C ***
  315.       INTEGER I, ID, IE, IERR, IFIX, II, INFO, IOPT, IP,
  316.      *     IPRINT, ISC, IT1, IT2, IT3, J, JE, JOB, JP, K, KALL,
  317.      *     KE, KP1, KRES, LAST, LDL, LDX, LIST, LKNT, LMAX,
  318.      *     LOC, LRES, MATZ, N, NO, P, P1, P12, P2, P2PJ
  319.       INTEGER DVAR(P), INT(P), IVAR(P), LVAR(LDL,LMAX),
  320.      *     PVAR(P), IDUM(1,1)
  321.       DOUBLE PRECISION BOUND, CRES, ENU2, PCERR, PROD, PSIK,
  322.      *     BIG, PSIMIN, RLKDR, RN, RNM1, RVAL, XLAMAX, DUMS,
  323.      *     TEMP
  324.       DOUBLE PRECISION CORR(P,P), COVD(P), DET(2),
  325.      *     OUTPUT(LMAX,5), R(P,P), V(P,P), W(P), WK(N), WK2(P),
  326.      *     X(LDX,P), XMU(P), XSAVE(N,P), DUM(1,1)
  327.       DOUBLE PRECISION         DLOG, DSIGN, DSQRT
  328.       REAL SNGL
  329. C ***
  330. C     PROLOG.
  331. C ***
  332.       ID = 1
  333.       IF (IPRINT.EQ.4) CALL PRINTM(4, X, LDX, N, P, IDUM, ID,
  334.      *     ID, ID)
  335.       RN = N
  336.       RNM1 = N-1
  337.       MATZ = 0
  338.       BIG = 1.D35
  339. C ***
  340. C     COMPUTE MEANS.
  341. C ***
  342.       DO 20 J=1,P
  343.            XMU(J) = 0.D0
  344.            DO 10 I=1,N
  345.                 XMU(J) = XMU(J) + X(I,J)
  346.    10      CONTINUE
  347.            XMU(J) = XMU(J)/RN
  348.    20 CONTINUE
  349.       IF (IPRINT.GE.2) CALL PRINTM(2, XMU, LDX, P, 1, IDUM, ID,
  350.      *     ID, ID)
  351. C ***
  352. C     COMPUTE COVARIANCE MATRIX.
  353. C ***
  354.       DO 50 I=1,P
  355.            DO 40 J=1,P
  356.                 CORR(I,J) = 0.D0
  357.                 DO 30 K=1,N
  358.                      CORR(I,J) = CORR(I,J) + X(K,I)*X(K,J)
  359.    30           CONTINUE
  360.                 CORR(I,J) = (CORR(I,J)-RN*XMU(I)*XMU(J))/RNM1
  361.    40      CONTINUE
  362.    50 CONTINUE
  363.       IF (IPRINT.GE.3) CALL PRINTM(3, CORR, P, P, P, IDUM, ID,
  364.      *     ID, ID)
  365. C ***
  366. C     CHECK FOR CONSTANT VARIABLES
  367. C ***
  368.       K = 0
  369.       DO 60 I=1,P
  370.            IF (CORR(I,I).NE.0.D0) GO TO 60
  371.            IDUM(1,1) = I
  372.            CALL PRINTM(36, DUM, ID, ID, ID, IDUM, 1, 1, 1)
  373.            K = 1
  374.    60 CONTINUE
  375.       IF (K.EQ.1) GO TO 600
  376. C ***
  377. C     COMPUTE CENTERED AND SCALED X, AND SAVE.
  378. C ***
  379.       DO 80 I=1,N
  380.            DO 70 J=1,P
  381.                 X(I,J) = (X(I,J)-XMU(J))/DSQRT(CORR(J,J))
  382.                 XSAVE(I,J) = X(I,J)
  383.    70      CONTINUE
  384.    80 CONTINUE
  385.       IDUM(1,1) = 1
  386.       IF (IPRINT.EQ.4) CALL PRINTM(9, X, LDX, N, P, IDUM, 1, 1,
  387.      *     1)
  388. C ***
  389. C     SAVE COVARIANCE DIAGONALS AND COMPUTE CORRELATION MATRIX.
  390. C ***
  391.       DO 90 I=1,P
  392.            COVD(I) = CORR(I,I)
  393.    90 CONTINUE
  394.       DO 110 I=1,P
  395.            DO 100 J=I,P
  396.                 CORR(I,J) = CORR(I,J)/DSQRT(COVD(I)*COVD(J))
  397.                 CORR(J,I) = CORR(I,J)
  398.   100      CONTINUE
  399.   110 CONTINUE
  400.       IF (IPRINT.GE.3) CALL PRINTM(6, CORR, P, P, P, IDUM, ID,
  401.      *     ID, ID)
  402. C ***
  403. C     COMPUTE SVD OF CENTERED, SCALED DATA MATRIX.
  404. C ***
  405.       IOPT = 1
  406.       CALL DSVDC(X, LDX, N, P, W, R, DUM, 1, X, LDX, WK, IOPT,
  407.      *     IERR)
  408.       IF (IERR.EQ.0) GO TO 120
  409.       CALL PRINTM(13, DUM, ID, ID, ID, IDUM, ID, ID, ID)
  410.       GO TO 600
  411.   120 DO 130 I=1,P
  412.            R(I,1) = W(I)
  413.            R(I,2) = (W(I)**2)/RNM1
  414.            W(I) = R(I,2)
  415.            R(I,3) = 100.D0
  416.            IF (W(I).GE..5D0) GO TO 130
  417.            R(I,3) = 100.D0*DSQRT(W(I)/(1.D0-W(I)))
  418.   130 CONTINUE
  419.       CALL PRINTM(7, R, P, P, 3, IDUM, ID, ID, ID)
  420.       TEMP = BIG
  421.       IF (R(P,1).NE.0.D0) TEMP = R(1,1)/R(P,1)
  422.       R(1,1) = TEMP
  423.       R(1,2) = BIG
  424.       IF (W(P).NE.0.D0) R(1,2) = W(1)/W(P)
  425.       CALL PRINTM(5, R, P, 1, 2, IDUM, ID, ID, ID)
  426.       IF (IPRINT.EQ.4) CALL PRINTM(8, X, LDX, P, P, IDUM, ID,
  427.      *     ID, ID)
  428.       IF (JOB.NE.1) GO TO 160
  429. C ***
  430. C     COMPUTE P2 FROM PER CENT ERROR.
  431. C ***
  432.       ENU2 = (PCERR/100.D0)**2
  433.       BOUND = ENU2/(1.D0+ENU2)
  434.       DO 140 I=1,P
  435.            II = P + 1 - I
  436.            IF (W(II).GT.BOUND) GO TO 150
  437.   140 CONTINUE
  438.   150 P2 = I - 1
  439.       R(1,1) = PCERR
  440.       R(2,1) = BOUND
  441.   160 IF (P2.GT.0) GO TO 170
  442.       CALL PRINTM(10, DUM, ID, ID, ID, IDUM, ID, ID, ID)
  443.       GO TO 600
  444.   170 IF (P2.LT.P) GO TO 180
  445.       CALL PRINTM(11, DUM, ID, ID, ID, IDUM, ID, ID, ID)
  446.       GO TO 600
  447.   180 P1 = P - P2
  448.       CALL PRINTM(12, R, P, 2, 1, IDUM, LMAX, JOB, P2)
  449.       IF (JOB.EQ.3) GO TO 310
  450. C ***
  451. C     COMPUTE NEARLY OPTIMAL SET OF ESTIMATED VARIABLES.
  452. C ***
  453.       DO 210 I=1,P2
  454.            II = P + 1 - I
  455.            DO 190 J=1,P
  456.                 V(I,J) = X(J,II)
  457.   190      CONTINUE
  458.            DO 200 J=1,P2
  459.                 R(I,J) = 0.D0
  460.   200      CONTINUE
  461.            R(I,I) = DSQRT(RNM1*W(II))
  462.   210 CONTINUE
  463.       DO 220 I=1,P
  464.            IVAR(I) = I
  465.   220 CONTINUE
  466.       IOPT = 0
  467.       IF (IPRINT.EQ.4) CALL PRINTM(29, V, P, P2, P, IDUM, ID,
  468.      *     ID, ID)
  469.       CALL PICVAR(P2, P, V, P, IOPT, IVAR, R, P, WK, WK2, RVAL)
  470. C ***
  471. C     COMPUTE ORDERED PREDICTOR VARIABLES.
  472. C ***
  473.       DO 240 J=1,P1
  474.            P2PJ = P2 + J
  475.            K = IVAR(P2PJ)
  476.            DO 230 I=1,P1
  477.                 V(I,J) = X(K,I)
  478.   230      CONTINUE
  479.   240 CONTINUE
  480.       IOPT = 1
  481.       CALL PICVAR(P1, P1, V, P, IOPT, IVAR(P2+1), DUM, 1, WK,
  482.      *     WK2, DUMS)
  483.       DUM(1,1) = RVAL
  484.       CALL PRINTM(14, DUM, 1, 1, 1, IVAR, P, P2, 1)
  485.       IF (KALL.EQ.1) GO TO 270
  486.       CALL PRINTM(34, DUM, ID, ID, ID, IDUM, ID, ID, ID)
  487.       CRES = RVAL
  488.       DO 250 I=1,P
  489.            PVAR(I) = IVAR(I)
  490.   250 CONTINUE
  491.       DO 260 I=1,P2
  492.            IT1 = P2 + 1 - I
  493.            DVAR(I) = IVAR(IT1)
  494.   260 CONTINUE
  495.       LKNT = 0
  496.       GO TO 380
  497. C ***
  498. C     ORDER VARIABLES.
  499. C ***
  500.   270 IF (P1.EQ.1) GO TO 290
  501.       P12 = P1/2
  502.       DO 280 I=1,P12
  503.            IT1 = P2 + I
  504.            IT2 = P + 1 - I
  505.            K = IVAR(IT1)
  506.            IVAR(IT1) = IVAR(IT2)
  507.            IVAR(IT2) = K
  508.   280 CONTINUE
  509. C ***
  510. C     START SEQUENTIAL VARIABLE REMOVAL PROCESS.
  511. C ***
  512.   290 DO 300 I=1,P2
  513.            INT(I) = P + 1 - I
  514.   300 CONTINUE
  515.       IDUM(1,1) = 1
  516.       CALL PRINTM(32, DUM, ID, ID, ID, IDUM, 1, 1, 1)
  517.   310 LKNT = 0
  518.   320 CONTINUE
  519.       DO 330 I=1,P2
  520.            IT1 = INT(I)
  521.            DVAR(I) = IVAR(IT1)
  522.            IF (JOB.EQ.3) DVAR(I) = LVAR(I,LKNT+1)
  523.   330 CONTINUE
  524.       DO 360 I=1,P2
  525.            II = P + 1 - I
  526.            DO 340 J=1,P
  527.                 V(I,J) = X(J,II)
  528.   340      CONTINUE
  529.            DO 350 J=1,P2
  530.                 R(I,J) = 0.D0
  531.   350      CONTINUE
  532.            R(I,I) = DSQRT(RNM1*W(II))
  533.   360 CONTINUE
  534.       DO 370 I=1,P
  535.            PVAR(I) = I
  536.   370 CONTINUE
  537.       IF (IPRINT.GE.3) CALL PRINTM(15, DUM, ID, ID, ID, DVAR,
  538.      *     P, P2, 1)
  539.       CALL RESFLY(P2, P, V, P, PVAR, DVAR, JOB, KRES, R, P, WK,
  540.      *     WK2, RVAL, LRES, IPRINT)
  541.       CRES = WK(1)
  542.       IF (LRES.EQ.0 .AND. CRES.LT.RVAL) RVAL = CRES
  543.       IF (IPRINT.EQ.4) CALL PRINTM(26, DUM, ID, ID, ID, PVAR,
  544.      *     P, P, 1)
  545.       DUM(1,1) = CRES
  546.       IF (IPRINT.GE.3 .AND. LRES.EQ.0) CALL PRINTM(16, DUM, 1,
  547.      *     1, 1, IDUM, ID, ID, ID)
  548. C ***
  549. C     CHECK EIGENVALUE INEQUALITIES FOR SMALL RESIDUALS AND PROCEED.
  550. C ***
  551.       IF (JOB.EQ.3) GO TO 380
  552.       IF (LRES.NE.0) GO TO 540
  553.   380 DO 400 I=1,P
  554.            IE = P + 1 - I
  555.            IP = PVAR(IE)
  556.            DO 390 J=1,P
  557.                 JE = P + 1 - J
  558.                 JP = PVAR(JE)
  559.                 V(I,J) = CORR(IP,JP)
  560.   390      CONTINUE
  561.   400 CONTINUE
  562.       IOPT = 0
  563.       CALL CHKEV(V, P, P, P1, IOPT, R, WK, WK2, PSIMIN, XLAMAX,
  564.      *     PSIK, IPRINT)
  565. C ***
  566. C     COMPUTE LIKELIHOOD RATIO.
  567. C ***
  568.       RLKDR = 0.D0
  569.       IF (N.LE.P) GO TO 430
  570.       DO 420 I=1,P2
  571.            IE = P2 + 1 - I
  572.            IP = PVAR(IE)
  573.            DO 410 J=I,P2
  574.                 JE = P2 + 1 - J
  575.                 JP = PVAR(JE)
  576.                 R(I,J) = V(I,J)*DSQRT(COVD(IP)*COVD(JP))
  577.                 R(J,I) = R(I,J)
  578.   410      CONTINUE
  579.            WK(I) = 1.D0/R(I,I)
  580.   420 CONTINUE
  581.       IF (IPRINT.EQ.4) CALL PRINTM(35, R, P, P2, P2, IDUM, ID,
  582.      *     ID, ID)
  583.       CALL SCPROD(WK, P2, PROD, ISC)
  584.       CALL DPOFA(R, P, P2, INFO)
  585.       IDUM(1,1) = INFO
  586.       IF (INFO.NE.0) CALL PRINTM(33, DUM, ID, ID, ID, IDUM, 1,
  587.      *     1, 1)
  588.       CALL DPODI(R, P, P2, DET, 10)
  589.       PROD = PROD*DET(1)
  590.       ISC = ISC + IFIX(SNGL(DET(2)+DSIGN(.1D0,DET(2))))
  591.       RLKDR = -RN*DLOG(PROD*10.D0**ISC)
  592.   430 LKNT = LKNT + 1
  593. C ***
  594. C     UPDATE ORDERED OUTPUT ARRAYS.
  595. C ***
  596.       LAST = LKNT - 1
  597.       IF (LKNT.GT.1) GO TO 440
  598.       LOC = 1
  599.       GO TO 510
  600.   440 IF (LKNT.LE.LMAX) GO TO 450
  601.       IF (CRES.GT.OUTPUT(LMAX,1)) GO TO 540
  602.       LAST = LMAX
  603.   450 DO 500 K=1,LAST
  604.            IF (CRES.GE.OUTPUT(K,1)) GO TO 500
  605.            KE = LKNT
  606.            IF (LKNT.GE.LMAX) KE = LMAX
  607.            IF (K.EQ.LMAX) GO TO 490
  608.            KP1 = K + 1
  609.            DO 480 I=KP1,KE
  610.                 IT1 = KE + KP1 - I
  611.                 IT2 = IT1 - 1
  612.                 DO 460 J=1,5
  613.                      OUTPUT(IT1,J) = OUTPUT(IT2,J)
  614.   460           CONTINUE
  615.                 DO 470 J=1,P2
  616.                      LVAR(J,IT1) = LVAR(J,IT2)
  617.   470           CONTINUE
  618.   480      CONTINUE
  619.   490      LOC = K
  620.            GO TO 510
  621.   500 CONTINUE
  622.       LOC = LKNT
  623.   510 DO 520 I=1,P2
  624.            II = P2 + 1 - I
  625.            LVAR(I,LOC) = DVAR(II)
  626.   520 CONTINUE
  627.       OUTPUT(LOC,1) = CRES
  628.       OUTPUT(LOC,2) = PSIK
  629.       TEMP = BIG
  630.       IF (XLAMAX.NE.0.D0) TEMP = PSIMIN/XLAMAX
  631.       OUTPUT(LOC,3) = TEMP
  632.       TEMP = 1.D0
  633.       IF (W(P1).NE.0.D0) TEMP = 0.D0
  634.       IF (W(P1+1).NE.0.D0) TEMP = OUTPUT(LOC,3)/(W(P1)/W(P1+1))
  635.       IF (XLAMAX.EQ.0.D0) TEMP = 1.D0
  636.       OUTPUT(LOC,4) = TEMP
  637.       OUTPUT(LOC,5) = RLKDR
  638.       IF (JOB.EQ.3) GO TO 530
  639.       IF (KALL.EQ.0) GO TO 570
  640.       GO TO 540
  641.   530 IF (LKNT.LT.LMAX) GO TO 320
  642.       GO TO 570
  643. C ***
  644. C     DETERMINE NEXT SET OF ESTIMATED VARIABLES TO TRY.
  645. C ***
  646.   540 DO 560 I=1,P2
  647.            IT1 = P2 + 1 - I
  648.            IF (LRES.NE.0 .AND. I.EQ.1) IT1 = LRES
  649.            IF (INT(IT1).EQ.I) GO TO 560
  650.            INT(IT1) = INT(IT1) - 1
  651.            IF (IT1.EQ.P2) GO TO 320
  652.            IT2 = IT1 + 1
  653.            DO 550 J=IT2,P2
  654.                 IT3 = J + 1 - IT2
  655.                 INT(J) = INT(IT1) - IT3
  656.   550      CONTINUE
  657.            GO TO 320
  658.   560 CONTINUE
  659.       IDUM(1,1) = 2
  660.       CALL PRINTM(32, DUM, ID, ID, ID, IDUM, 1, 1, 1)
  661.   570 LIST = LMAX
  662.       IF (LIST.GT.LKNT) LIST = LKNT
  663.       CALL PRINTM(17, OUTPUT, LMAX, LIST, 5, LVAR, LDL, P2,
  664.      *     LIST)
  665. C ***
  666. C     COMPUTE BETA FOR BEST LMAX RESIDUAL CASES.
  667. C ***
  668.       DO 590 I=1,P
  669.            DO 580 J=I,P
  670.                 CORR(I,J) = CORR(I,J)*DSQRT(COVD(I)*COVD(J))
  671.                 CORR(J,I) = CORR(I,J)
  672.   580      CONTINUE
  673.   590 CONTINUE
  674.       IF (IPRINT.EQ.4) CALL PRINTM(3, CORR, P, P, P, IDUM, ID,
  675.      *     ID, ID)
  676.       NO = LIST
  677.       CALL BETA(LVAR, LDL, NO, XSAVE, N, P, CORR, P, P2, XMU,
  678.      *     X, LDX, V, P, W, WK, WK2, PVAR, IPRINT)
  679.   600 CONTINUE
  680.       RETURN
  681.       END
  682. C
  683. C --------------------------------------------------------------
  684. C
  685.       SUBROUTINE PICVAR(P2, P, V, MV, IOPT, PVAR, RES, MR, WO,
  686.      *     WN, RNORM)
  687. C ***
  688. C     SUBROUTINE PICVAR COMPUTES THE BENCHMARK SOLUTION BY A QR
  689. C        FACTORIZATION. IF IOPT = 0, THE ESTIMATED VARIABLES ARE
  690. C        COMPUTED; IF IOPT = 1, THE PREDICTOR VARIABLES ARE COMPUTED.
  691. C ***
  692.       INTEGER I, IOPT, J, JP, K, KM1, L, LEN, LP1, MAXJ, MR,
  693.      *     MV, P, P2, P2P1
  694.       INTEGER PVAR(P)
  695.       DOUBLE PRECISION CNORM, FAC, GNORM, ONEPS, PROD, RNORM
  696.       DOUBLE PRECISION RES(MR,P2), V(MV,P), WN(P), WO(P)
  697.       DOUBLE PRECISION DABS, DDOT, DMAX1, DNRM2, DSIGN, DSQRT
  698. C ***
  699. C     COMPUTE COLUMN NORMS.
  700. C ***
  701.       DO 10 J=1,P
  702.            WN(J) = DNRM2(P2,V(1,J),1)
  703.            WO(J) = WN(J)
  704.    10 CONTINUE
  705. C ***
  706. C     REDUCE TO TRAPEZOIDAL FORM.
  707. C ***
  708.       DO 70 L=1,P2
  709. C ***
  710. C        PIVOT COLUMNS.
  711. C ***
  712.            GNORM = 0.D0
  713.            MAXJ = L
  714.            DO 20 J=L,P
  715.                 IF (WN(J).LE.GNORM) GO TO 20
  716.                 GNORM = WN(J)
  717.                 MAXJ = J
  718.    20      CONTINUE
  719.            IF (MAXJ.EQ.L) GO TO 30
  720.            CALL DSWAP(P2, V(1,L), 1, V(1,MAXJ), 1)
  721.            WN(MAXJ) = WN(L)
  722.            WO(MAXJ) = WO(L)
  723.            JP = PVAR(MAXJ)
  724.            PVAR(MAXJ) = PVAR(L)
  725.            PVAR(L) = JP
  726.    30      CONTINUE
  727.            WN(L) = 0.D0
  728.            IF (L.EQ.P2) GO TO 70
  729. C ***
  730. C        REDUCE COLUMN L.
  731. C ***
  732.            LEN = P2 - L + 1
  733.            CNORM = DNRM2(LEN,V(L,L),1)
  734.            IF (CNORM.EQ.0.D0) GO TO 70
  735.            IF (V(L,L).NE.0.D0) CNORM = DSIGN(CNORM,V(L,L))
  736.            CALL DSCAL(LEN, 1.D0/CNORM, V(L,L), 1)
  737.            V(L,L) = 1.D0 + V(L,L)
  738. C ***
  739. C        APPLY TRANSFORMATION TO OTHER COLUMNS.
  740. C ***
  741.            LP1 = L + 1
  742.            DO 50 J=LP1,P
  743.                 PROD = -DDOT(LEN,V(L,L),1,V(L,J),1)/V(L,L)
  744.                 CALL DAXPY(LEN, PROD, V(L,L), 1, V(L,J), 1)
  745.                 IF (WN(J).EQ.0.D0) GO TO 50
  746.                 FAC = 1.D0 - (DABS(V(L,J))/WN(J))**2
  747.                 FAC = DMAX1(FAC,0.D0)
  748.                 ONEPS = 1.D0 + .05D0*FAC*(WN(J)/WO(J))**2
  749.                 IF (ONEPS.EQ.1.D0) GO TO 40
  750.                 WN(J) = WN(J)*DSQRT(FAC)
  751.                 GO TO 50
  752.    40           CONTINUE
  753.                 WN(J) = DNRM2(P2-L,V(L+1,J),1)
  754.                 WO(J) = WN(J)
  755.    50      CONTINUE
  756.            IF (IOPT.EQ.1) GO TO 70
  757. C ***
  758. C        APPLY TRANSFORMATION TO RESIDUAL MATRIX.
  759. C ***
  760.            DO 60 J=1,P2
  761.                 PROD = -DDOT(LEN,V(L,L),1,RES(L,J),1)/V(L,L)
  762.                 CALL DAXPY(LEN, PROD, V(L,L), 1, RES(L,J), 1)
  763.    60      CONTINUE
  764.            WN(L) = V(L,L)
  765.            V(L,L) = -CNORM
  766.    70 CONTINUE
  767.       IF (IOPT.EQ.1) GO TO 140
  768. C ***
  769. C     COMPUTE DEPENDENCIES.
  770. C ***
  771.       P2P1 = P2 + 1
  772.       DO 90 J=P2P1,P
  773.            V(P2,J) = V(P2,J)/V(P2,P2)
  774.            DO 80 L=2,P2
  775.                 K = P2 - L + 2
  776.                 KM1 = K - 1
  777.                 CALL DAXPY(KM1, -V(K,J), V(1,K), 1, V(1,J), 1)
  778.                 V(KM1,J) = V(KM1,J)/V(KM1,KM1)
  779.    80      CONTINUE
  780.    90 CONTINUE
  781. C ***
  782. C     APPLY INVERSE TO RESIDUAL MATRIX.
  783. C ***
  784.       DO 110 J=1,P2
  785.            RES(P2,J) = RES(P2,J)/V(P2,P2)
  786.            DO 100 L=2,P2
  787.                 K = P2 - L + 2
  788.                 KM1 = K - 1
  789.                 CALL DAXPY(KM1, -RES(K,J), V(1,K), 1, RES(1,J),
  790.      *               1)
  791.                 RES(KM1,J) = RES(KM1,J)/V(KM1,KM1)
  792.   100      CONTINUE
  793.   110 CONTINUE
  794. C ***
  795. C     COMPUTE RESIDUAL NORM.
  796. C ***
  797.       RNORM = 0.D0
  798.       DO 130 I=1,P2
  799.            DO 120 J=1,P2
  800.                 RNORM = RNORM + RES(I,J)**2
  801.   120      CONTINUE
  802.   130 CONTINUE
  803.       RNORM = DSQRT(RNORM)
  804.   140 CONTINUE
  805.       RETURN
  806.       END
  807. C
  808. C --------------------------------------------------------------
  809. C
  810.       SUBROUTINE RESFLY(P2, P, V, MV, PVAR, DVAR, JOB, KRES,
  811.      *     RES, MR, WN, VT, RMIN, LRES, IPRINT)
  812. C ***
  813. C     SUBROUTINE RESFLY COMPUTES THE REDUCED SPACE RESIDUAL FOR A
  814. C        PARTICULAR SET OF ESTIMATED VARIABLES BY A SEQUENTIAL VARIABLE
  815. C        DELECTION PROCESS. LARGE RESIDUALS ARE DETECTED EARLY IF
  816. C        POSSIBLE.
  817. C ***
  818.       INTEGER I, ID, IPRINT, J, JOB, JP, K, KRES, L, LEN, LL,
  819.      *     LM1, LM2, LMI, LMI1, LRES, MR, MV, NJ, P, P2
  820.       INTEGER DVAR(P2), PVAR(P), IDUM(1,1)
  821.       DOUBLE PRECISION CNORM, PROD, RMIN, RNORM, TEMP, TOOBIG
  822.       DOUBLE PRECISION RES(MR,P2), V(MV,P), VT(P2), WN(P),
  823.      *     DUM(1,1)
  824.       DOUBLE PRECISION DDOT, DNRM2, DSIGN, DSQRT
  825.       ID = 1
  826.       LRES = 0
  827.       IF (JOB.NE.3) TOOBIG = RMIN*(10.D0**KRES)
  828. C ***
  829. C     REDUCE TO TRAPEZOIDAL FORM WHILE UPDATING RESIDUALS.
  830. C ***
  831.       DO 160 L=1,P2
  832.            NJ = DVAR(L)
  833.            DO 10 LL=L,P
  834.                 IF (NJ.EQ.PVAR(LL)) GO TO 20
  835.    10      CONTINUE
  836.    20      NJ = LL
  837.            IF (NJ.EQ.L) GO TO 30
  838. C ***
  839. C        PIVOT COLUMNS.
  840. C ***
  841.            CALL DSWAP(P2, V(1,L), 1, V(1,NJ), 1)
  842.            JP = PVAR(NJ)
  843.            PVAR(NJ) = PVAR(L)
  844.            PVAR(L) = JP
  845.    30      CONTINUE
  846.            IF (L.EQ.1) GO TO 50
  847. C ***
  848. C      APPLY PREVIOUS TRANSFORMATIONS TO COLUMN L OF V AND RES.
  849. C ***
  850.            LM1 = L - 1
  851.            DO 40 K=1,LM1
  852.                 LEN = P2 - K + 1
  853.                 TEMP = V(K,K)
  854.                 V(K,K) = WN(K)
  855.                 PROD = -DDOT(LEN,V(K,K),1,V(K,L),1)/V(K,K)
  856.                 CALL DAXPY(LEN, PROD, V(K,K), 1, V(K,L), 1)
  857.                 PROD = -DDOT(LEN,V(K,K),1,RES(K,L),1)/V(K,K)
  858.                 CALL DAXPY(LEN, PROD, V(K,K), 1, RES(K,L), 1)
  859.                 V(K,K) = TEMP
  860.    40      CONTINUE
  861.    50      CONTINUE
  862.            IF (L.LT.P2) GO TO 60
  863.            IF (V(L,L).EQ.0.D0) GO TO 170
  864.            GO TO 80
  865. C ***
  866. C        REDUCE COLUMN L.
  867. C ***
  868.    60      LEN = P2 - L + 1
  869.            CNORM = DNRM2(LEN,V(L,L),1)
  870.            IF (CNORM.EQ.0.D0) GO TO 170
  871.            IF (V(L,L).NE.0.D0) CNORM = DSIGN(CNORM,V(L,L))
  872.            CALL DSCAL(LEN, 1.D0/CNORM, V(L,L), 1)
  873.            V(L,L) = 1.D0 + V(L,L)
  874. C ***
  875. C        APPLY TRANSFORMATION TO FIRST L COLUMNS OF RES.
  876. C ***
  877.            DO 70 J=1,L
  878.                 PROD = -DDOT(LEN,V(L,L),1,RES(L,J),1)/V(L,L)
  879.                 CALL DAXPY(LEN, PROD, V(L,L), 1, RES(L,J), 1)
  880.    70      CONTINUE
  881.            WN(L) = V(L,L)
  882.            V(L,L) = -CNORM
  883.    80      CONTINUE
  884. C ***
  885. C        COMPUTE FIRST L COMPONENTS IN COLUMN L OF RESIDUAL MATRIX.
  886. C ***
  887.            RES(L,L) = RES(L,L)/V(L,L)
  888.            IF (L.EQ.1) GO TO 130
  889.            CALL DAXPY(LM1, -RES(L,L), V(1,L), 1, RES(1,L), 1)
  890.            RES(LM1,L) = RES(LM1,L)/V(LM1,LM1)
  891.            VT(LM1) = V(LM1,L)/V(LM1,LM1)
  892.            IF (L.EQ.2) GO TO 110
  893.            LM2 = L - 2
  894.            DO 90 I=1,LM2
  895.                 VT(I) = V(I,L)
  896.    90      CONTINUE
  897.            DO 100 I=2,LM1
  898.                 LMI1 = LM1 - I + 2
  899.                 LMI = LMI1 - 1
  900.                 CALL DAXPY(LMI, -RES(LMI1,L), V(1,LMI1), 1,
  901.      *               RES(1,L), 1)
  902.                 RES(LMI,L) = RES(LMI,L)/V(LMI,LMI)
  903.                 CALL DAXPY(LMI, -VT(LMI1), V(1,LMI1), 1, VT(1),
  904.      *               1)
  905.                 VT(LMI) = VT(LMI)/V(LMI,LMI)
  906.   100      CONTINUE
  907.   110      CONTINUE
  908. C ***
  909. C        COMPUTE TOP LEFT L BY L-1 RESIDUAL SUBMATRIX.
  910. C ***
  911.            CALL DSCAL(LM1, 1.D0/V(L,L), RES(L,1), MR)
  912.            DO 120 I=1,LM1
  913.                 CALL DAXPY(LM1, -VT(I), RES(L,1), MR, RES(I,1),
  914.      *               MR)
  915.   120      CONTINUE
  916.   130      CONTINUE
  917. C ***
  918. C        COMPUTE SIZE OF CURRENT RESIDUAL MATRIX.
  919. C ***
  920.            RNORM = 0.D0
  921.            DO 150 I=1,L
  922.                 DO 140 J=1,L
  923.                      RNORM = RNORM + RES(I,J)**2
  924.   140           CONTINUE
  925.   150      CONTINUE
  926.            RNORM = DSQRT(RNORM)
  927.            IDUM(1,1) = L
  928.            DUM(1,1) = RNORM
  929.            IF (IPRINT.GE.3) CALL PRINTM(27, DUM, 1, 1, 1, IDUM,
  930.      *          1, 1, 1)
  931.            IF (JOB.EQ.3) GO TO 160
  932.            IF (RNORM.GT.TOOBIG) GO TO 180
  933.   160 CONTINUE
  934.       WN(1) = RNORM
  935.       GO TO 190
  936.   170 IF (JOB.NE.3) GO TO 180
  937. C ***
  938. C     THE FOLLOWING NUMBER REPRESENTS AN INFINITE RESIDUAL NORM. IT
  939. C         SHOULD NOT PRODUCE OVERFLOW ON MOST KNOWN COMPUTERS.
  940. C
  941.       WN(1) = 1.E+30
  942. C
  943. C ***
  944.       CALL PRINTM(37, DUM, ID, ID, ID, DVAR, P2, L, 1)
  945.   180 LRES = L
  946.       IDUM(1,1) = L
  947.       IF (IPRINT.GE.3) CALL PRINTM(28, DUM, ID, ID, ID, IDUM,
  948.      *     1, 1, 1)
  949.   190 CONTINUE
  950.       RETURN
  951.       END
  952. C
  953. C --------------------------------------------------------------
  954. C
  955.       SUBROUTINE PRINTM(KODE, A, LD, M, N, IA, LDI, IM, IN)
  956. C ***
  957. C     SUBROUTINE PRINTM PRINTS MOST OF THE OUTPUT FROM THE LINEAR
  958. C        DEPENDENCY ANALYSIS. (THE REMAINING OUTPUT IS PRINTED IN
  959. C        SUBROUTINES LDA AND BETA.)
  960. C ***
  961.       INTEGER I, IE, IEP1, IM, IMP, IN, ISKIP, ISYM, J, K,
  962.      *     KODE, LD, LDI, M, N, NCOLE, NCOLS, P1
  963.       INTEGER IA(LDI,IN)
  964.       DOUBLE PRECISION A(LD,N), CUM, SUM, XP
  965.       DATA IBLK /1H /, IAST /1H*/
  966.       GO TO (10, 20, 40, 50, 60, 70, 80, 100, 110, 120, 130,
  967.      *     140, 170, 180, 190, 200, 210, 210, 240, 250, 260,
  968.      *     270, 280, 290, 300, 310, 320, 330, 340, 350, 360,
  969.      *     370, 380, 390, 400, 410, 420), KODE
  970. C
  971.    10 WRITE (6,99998)
  972.       GO TO 470
  973. C
  974.    20 WRITE (6,99997)
  975.       DO 30 I=1,M
  976.            WRITE (6,99996) A(I,1)
  977.    30 CONTINUE
  978.       WRITE (6,99995)
  979.       GO TO 470
  980. C
  981.    40 WRITE (6,99994)
  982.       GO TO 430
  983. C
  984.    50 WRITE (6,99993)
  985.       GO TO 430
  986. C
  987.    60 WRITE (6,99992) A(1,1), A(1,2)
  988.       WRITE (6,99995)
  989.       GO TO 470
  990. C
  991.    70 WRITE (6,99991)
  992.       GO TO 430
  993. C
  994.    80 WRITE (6,99990)
  995.       XP = M
  996.       SUM = 0.D0
  997.       DO 90 I=1,M
  998.            SUM = SUM + A(I,2)
  999.            CUM = (SUM/XP)*100.D0
  1000.            WRITE (6,99989) A(I,1), A(I,2), A(I,3), CUM
  1001.    90 CONTINUE
  1002.       WRITE (6,99995)
  1003.       GO TO 470
  1004. C
  1005.   100 WRITE (6,99988)
  1006.       GO TO 430
  1007. C
  1008.   110 IF (IA(1,1).EQ.1) WRITE (6,99987)
  1009.       IF (IA(1,1).EQ.2) WRITE (6,99986)
  1010.       IF (IA(1,1).EQ.3) WRITE (6,99985) IM, IN
  1011.       GO TO 430
  1012. C
  1013.   120 WRITE (6,99984)
  1014.       GO TO 470
  1015. C
  1016.   130 WRITE (6,99983)
  1017.       GO TO 470
  1018. C
  1019.   140 WRITE (6,99982) IN
  1020.       IF (IM.NE.1) GO TO 150
  1021.       WRITE (6,99981) A(1,1), A(2,1)
  1022.       GO TO 470
  1023.   150 IF (IM.NE.2) GO TO 160
  1024.       WRITE (6,99980)
  1025.       GO TO 470
  1026.   160 WRITE (6,99979) LDI
  1027.       GO TO 470
  1028. C
  1029.   170 WRITE (6,99978)
  1030.       GO TO 470
  1031. C
  1032.   180 WRITE (6,99977) IM
  1033.       IE = IM
  1034.       IF (IM.GT.33) IE = 33
  1035.       WRITE (6,99976) (IA(I,1),I=1,IE)
  1036.       IF (IM.GT.33) WRITE (6,99975) (IA(I,1),I=34,IM)
  1037.       IMP = IM + 1
  1038.       IE = LDI
  1039.       P1 = LDI - IM
  1040.       IF (P1.GT.33) IE = IM + 33
  1041.       WRITE (6,99974) (IA(I,1),I=IMP,IE)
  1042.       IEP1 = IE + 1
  1043.       IF (P1.GT.33) WRITE (6,99975) (IA(I,1),I=IEP1,LDI)
  1044.       WRITE (6,99973) A(1,1)
  1045.       WRITE (6,99995)
  1046.       GO TO 470
  1047. C
  1048.   190 IE = IM
  1049.       IF (IM.GT.33) IE = 33
  1050.       WRITE (6,99972) (IA(I,1),I=1,IE)
  1051.       IF (IM.GT.33) WRITE (6,99975) (IA(I,1),I=34,IM)
  1052.       GO TO 470
  1053. C
  1054.   200 WRITE (6,99971) A(1,1)
  1055.       GO TO 470
  1056. C
  1057.   210 WRITE (6,99970)
  1058.       ISKIP = 1
  1059.       DO 230 I=1,M
  1060.            ISYM = IBLK
  1061.            IF (A(I,3).GT.1.D0) GO TO 220
  1062.            ISYM = IAST
  1063.            ISKIP = 0
  1064.   220      IE = IM
  1065.            IF (IM.GT.15) IE = 15
  1066.            WRITE (6,99969) (A(I,J),J=1,3), ISYM,
  1067.      *          (A(I,J),J=4,5), (IA(J,I),J=1,IE)
  1068.            IF (IM.GT.15) WRITE (6,99967) (IA(J,I),J=16,IM)
  1069.   230 CONTINUE
  1070.       IF (ISKIP.EQ.0) WRITE (6,99968)
  1071.       WRITE (6,99966)
  1072.       GO TO 470
  1073. C
  1074.   240 WRITE (6,99965) IA(1,1)
  1075.       GO TO 470
  1076. C
  1077.   250 WRITE (6,99964) IA(1,1)
  1078.       GO TO 470
  1079. C
  1080.   260 WRITE (6,99963)
  1081.       GO TO 430
  1082. C
  1083.   270 WRITE (6,99962)
  1084.       GO TO 430
  1085. C
  1086.   280 WRITE (6,99961)
  1087.       GO TO 430
  1088. C
  1089.   290 WRITE (6,99960)
  1090.       WRITE (6,99996) (A(I,1),I=1,M)
  1091.       GO TO 470
  1092. C
  1093.   300 WRITE (6,99959)
  1094.       WRITE (6,99996) (A(I,1),I=1,M)
  1095.       GO TO 470
  1096. C
  1097.   310 IE = IM
  1098.       IF (IM.GT.30) IE = 30
  1099.       WRITE (6,99958) (IA(I,1),I=1,IE)
  1100.       IF (IM.GT.30) WRITE (6,99957) (IA(I,1),I=31,IM)
  1101.       GO TO 470
  1102. C
  1103.   320 WRITE (6,99956) IA(1,1), A(1,1)
  1104.       GO TO 470
  1105. C
  1106.   330 WRITE (6,99955) IA(1,1)
  1107.       GO TO 470
  1108. C
  1109.   340 WRITE (6,99954)
  1110.       GO TO 430
  1111. C
  1112.   350 IF (IA(1,1).EQ.1) WRITE (6,99953)
  1113.       IF (IA(1,1).EQ.2) WRITE (6,99952)
  1114.       GO TO 430
  1115. C
  1116.   360 WRITE (6,99951) A(1,1)
  1117.       GO TO 470
  1118. C
  1119.   370 IF (IA(1,1).EQ.1) WRITE (6,99950)
  1120.       IF (IA(1,1).EQ.2) WRITE (6,99949)
  1121.       GO TO 470
  1122. C
  1123.   380 WRITE (6,99948) IA(1,1)
  1124.       GO TO 470
  1125. C
  1126.   390 WRITE (6,99947)
  1127.       GO TO 470
  1128. C
  1129.   400 WRITE (6,99946)
  1130.       GO TO 430
  1131. C
  1132.   410 WRITE (6,99999) IA(1,1)
  1133.       GO TO 470
  1134. C
  1135.   420 WRITE (6,99943)
  1136.       WRITE (6,99975) (IA(I,1),I=1,IM)
  1137.       GO TO 470
  1138. C ***
  1139. C     MATRIX PRINT UTILITY.
  1140. C ***
  1141.   430 DO 450 K=1,N,4
  1142.            NCOLS = K
  1143.            NCOLE = K + 3
  1144.            IF (NCOLE.GT.N) NCOLE = N
  1145.            WRITE (6,99945) NCOLS, NCOLE
  1146.            DO 440 I=1,M
  1147.                 WRITE (6,99944) (A(I,J),J=NCOLS,NCOLE)
  1148.   440      CONTINUE
  1149.            IF (NCOLE.EQ.N) GO TO 460
  1150.   450 CONTINUE
  1151.   460 WRITE (6,99995)
  1152.   470 RETURN
  1153. 99999 FORMAT (//9H VARIABLE, I3, 27H  IS CONSTANT FOR ALL OBSER,
  1154.      *     7HVATIONS, 38H *** CODE TERMINATES AFTER CHECKING OT,
  1155.      *     13HHER VARIABLES)
  1156. 99998 FORMAT (///1H1, 20X, 34H*** LINEAR DEPENDENCY ANALYSIS ***
  1157.      *     ///)
  1158. 99997 FORMAT (28H THE VECTOR OF MEANS IS...../)
  1159. 99996 FORMAT (1PD16.6)
  1160. 99995 FORMAT (//)
  1161. 99994 FORMAT (30H THE COVARIANCE MATRIX IS.....)
  1162. 99993 FORMAT (33H THE ORIGINAL DATA MATRIX IS.....)
  1163. 99992 FORMAT (19H CONDITION NUMBERS., 4X, 18HCENTERED, SCALED X,
  1164.      *     12H MATRIX....., 1PD12.5/30X, 19HCORRELATION MATRIX.,
  1165.      *     4H...., 1PD12.5/)
  1166. 99991 FORMAT (31H THE CORRELATION MATRIX IS.....)
  1167. 99990 FORMAT (1X, 18HSINGULAR VALUES OF, 8X, 14HEIGENVALUES OF,
  1168.      *     7X, 14HPERCENT ERRORS, 4X, 21HEIGENVALUE CUMULATIVE/
  1169.      *     1X, 18HCENTERED, SCALED X, 6X, 18HCORRELATION MATRIX,
  1170.      *     4X, 16H(DIAGONAL MODEL), 9X, 10HPERCENTAGE/)
  1171. 99989 FORMAT (2X, 1PD14.6, 10X, 1PD14.6, 12X, 0PF6.1, 15X, F6.1)
  1172. 99988 FORMAT (46H SINGULAR VECTORS OF CORRELATION MATRIX ARE...,
  1173.      *     2H..)
  1174. 99987 FORMAT (/44H THE CENTERED AND SCALED DATA MATRIX IS.....)
  1175. 99986 FORMAT (/45H THE PERMUTED CENTERED AND SCALED DATA MATRIX,
  1176.      *     8H IS.....)
  1177. 99985 FORMAT (/26H THE DATA ARRAY CONTAINING/4X, 10H THE FIRST,
  1178.      *     I4, 32H COLUMNS OF PREDICTOR VARIABLES,/4X, 6H THE L,
  1179.      *     4HAST , I4, 22H COLUMNS OF RESIDUALS.)
  1180. 99984 FORMAT (46H ERROR.  COMPUTED P2 LESS THAN OR EQUAL TO ZER,
  1181.      *     2HO.)
  1182. 99983 FORMAT (46H ERROR.  COMPUTED P2 GREATER THAN OR EQUAL TO ,
  1183.      *     2HP.)
  1184. 99982 FORMAT (/30H NUMBER OF DEPENDENCIES (P2) =, I3)
  1185. 99981 FORMAT (6X, 29HDETERMINED FROM USER SUPPLIED, 9H ALLOWABL,
  1186.      *     21HE PERCENTAGE ERROR OF, F6.1/6X, 8HYIELDING,
  1187.      *     34H A BOUND ON THE SINGULAR VALUES OF, 1PD14.6/)
  1188. 99980 FORMAT (6X, 17HSPECIFIED BY USER/)
  1189. 99979 FORMAT (6X, 28HSPECIFIED BY USER ALONG WITH, I3, 6H SETS ,
  1190.      *     30HOF ESTIMATED VARIABLES TO TEST/)
  1191. 99978 FORMAT (18H ERROR FROM DSVDC./)
  1192. 99977 FORMAT (/28H BENCHMARK SOLUTION FOR P2 =, I3//6X,
  1193.      *     13HTHE FOLLOWING, 31H TWO SETS OF VARIABLES ARE ORDE,
  1194.      *     4HRED./6X, 15HTHE MOST LIKELY, 18H ONES TO BE IN THA,
  1195.      *     23HT SET ARE LISTED FIRST./)
  1196. 99976 FORMAT (6X, 23HP2 ESTIMATED VARIABLES-, 2X, 33I3)
  1197. 99975 FORMAT (31X, 33I3)
  1198. 99974 FORMAT (6X, 23HP1 PREDICTOR VARIABLES-, 2X, 33I3)
  1199. 99973 FORMAT (/6X, 41HAPPROXIMATE (REDUCED SPACE) RESIDUAL IS..,
  1200.      *     1H., 2H.., 1PD16.6)
  1201. 99972 FORMAT (//29H***** FOR ESTIMATED VARIABLES, 2X, 33I3/)
  1202. 99971 FORMAT (31H REDUCED SPACE RESIDUAL IS....., 1PD16.8/)
  1203. 99970 FORMAT (/////1X, 13HREDUCED SPACE, 4X, 13HPSI CONDITION,
  1204.      *     22X, 10HPSI/LAMBDA/3X, 8HRESIDUAL, 11X, 6HNUMBER,
  1205.      *     6X, 14HPSI/LAMBDA GAP, 3X, 14HNORMALIZED GAP, 2X,
  1206.      *     16HLIKELIHOOD RATIO, 3X, 19HESTIMATED VARIABLES/)
  1207. 99969 FORMAT (1X, 2(1PD12.5, 5X), 1PD12.5, 1X, A1, 3X,
  1208.      *     2(1PD12.5, 5X), 15I3)
  1209. 99968 FORMAT (/45H *  DENOTES CASES WHERE LARGEST EIGENVALUE OF,
  1210.      *     1H , 42HLAMBDA EXCEEDS SMALLEST EIGENVALUE OF PSI.)
  1211. 99967 FORMAT (86X, 15I3)
  1212. 99966 FORMAT (///)
  1213. 99965 FORMAT (46H MATRIX NOT POSITIVE DEFINITE IN CHKEV.  IERR ,
  1214.      *     1H=, I4)
  1215. 99964 FORMAT (44H ERROR RETURN FROM EISPACK IN CHKEV.  IERR =,
  1216.      *     I4)
  1217. 99963 FORMAT (23H THE PSI MATRIX IS.....)
  1218. 99962 FORMAT (26H THE LAMBDA MATRIX IS.....)
  1219. 99961 FORMAT (46H THE PARTIALLY FACTORIZED SIGMA MATRIX IS.....)
  1220. 99960 FORMAT (35H THE EIGENVALUES OF LAMBDA ARE...../)
  1221. 99959 FORMAT (32H THE EIGENVALUES OF PSI ARE...../)
  1222. 99958 FORMAT (37H VARIABLE ORDERING OUT OF RESFLY....., 30I3)
  1223. 99957 FORMAT (37X, 30I3)
  1224. 99956 FORMAT (7H STEP =, I3, 10X, 10HRESIDUAL =, 1PD15.7)
  1225. 99955 FORMAT (32H LARGE RESIDUAL DETECTED AT STEP, I3)
  1226. 99954 FORMAT (46H THE SINGULAR VECTORS ENTERING PICVAR ARE.....)
  1227. 99953 FORMAT (40H THE PERMUTED CORRELATION MATRIX IS.....)
  1228. 99952 FORMAT (39H THE PERMUTED COVARIANCE MATRIX IS.....)
  1229. 99951 FORMAT (//17H RESIDUAL IS....., 1PD16.8/)
  1230. 99950 FORMAT (///43H BEGIN EXAMINING ALL POSSIBLE SUBSETS OF P2,
  1231.      *     10H VARIABLES/)
  1232. 99949 FORMAT (/46H ALL POSSIBLE SUBSETS OF P2 VARIABLES EXAMINED
  1233.      *     ///)
  1234. 99948 FORMAT (/8H ERROR =, I3, 3X, 25HFROM DPOFA IN LIKELIHOOD ,
  1235.      *     6HRATIO , 11HCOMPUTATION)
  1236. 99947 FORMAT (/45H ONLY BENCHMARK SOLUTION REQUESTED. ALL POSSI,
  1237.      *     3HBLE, 21H SUBSETS NOT EXAMINED///)
  1238. 99946 FORMAT (/35H THE UNSCALED LAMBDA MATRIX IS.....)
  1239. 99945 FORMAT (/8H COLUMNS, I5, 8H THROUGH, I5, 5H...../)
  1240. 99944 FORMAT (4(1PD16.6))
  1241. 99943 FORMAT (46H SINGULAR V22 RESULTS FROM ESTIMATED VARIABLES,
  1242.      *     5H.....)
  1243.       END
  1244. C
  1245. C --------------------------------------------------------------
  1246. C
  1247.       SUBROUTINE CHKEV(SIG, LDS, P, P1, IOPT, W, TEM1, TEM2,
  1248.      *     PSIMIN, XLAMAX, PSIK, IPRINT)
  1249. C ***
  1250. C     SUBROUTINE CHKEV COMPUTES THE PSI AND LAMBDA MATRICES FOR A
  1251. C        PARTICULAR SET OF ESTIMATED VARIABLES. THE MINIMUM EIGENVALUE
  1252. C        OF PSI AND THE MAXIMUM EIGENVALUE OF LAMBDA ARE COMPUTED.
  1253. C ***
  1254.       INTEGER I, ID, IERR, IOPT, IP1, IPRINT, J, JM1, JP1, K,
  1255.      *     LDS, LEN, MATZ, P, P1, P2
  1256.       INTEGER IDUM(1,1)
  1257.       DOUBLE PRECISION PSIK, PSIMIN, S, T, XLAMAX
  1258.       DOUBLE PRECISION SIG(LDS,P), TEM1(P), TEM2(P), W(P),
  1259.      *     DUM(1,1)
  1260.       DOUBLE PRECISION DDOT, DSQRT
  1261.       ID = 1
  1262.       P2 = P - P1
  1263.       IDUM(1,1) = 1
  1264.       IF (IPRINT.EQ.4 .AND. IOPT.EQ.0) CALL PRINTM(30, SIG,
  1265.      *     LDS, P, P, IDUM, 1, 1, 1)
  1266.       IDUM(1,1) = 2
  1267.       IF (IPRINT.EQ.4 .AND. IOPT.EQ.1) CALL PRINTM(30, SIG,
  1268.      *     LDS, P, P, IDUM, 1, 1, 1)
  1269.       MATZ = 0
  1270. C ***
  1271. C     COMPUTE PSI (CORRELATION MATRIX OF ORDER P1) AND
  1272. C             LAMBDA (SCHUR COMPLEMENT OF PSI).
  1273. C ***
  1274.       DO 30 J=1,P
  1275.            IERR = J
  1276.            S = 0.D0
  1277.            JM1 = J - 1
  1278.            IF (JM1.LT.1) GO TO 20
  1279.            DO 10 K=1,JM1
  1280.                 LEN = K - 1
  1281.                 IF (K.GT.P1) LEN = P1
  1282.                 T = SIG(K,J) - DDOT(LEN,SIG(1,K),1,SIG(1,J),1)
  1283.                 IF (K.LE.P1) T = T/SIG(K,K)
  1284.                 SIG(K,J) = T
  1285.                 IF (K.LE.P1) S = S + T*T
  1286.    10      CONTINUE
  1287.    20      S = SIG(J,J) - S
  1288.            SIG(J,J) = S
  1289.            IF (J.GT.P1) GO TO 30
  1290.            IF (S.LE.0.D0) GO TO 40
  1291.            SIG(J,J) = DSQRT(S)
  1292.    30 CONTINUE
  1293.       IERR = 0
  1294.       IF (IPRINT.EQ.4) CALL PRINTM(23, SIG, LDS, P, P, IDUM,
  1295.      *     ID, ID, ID)
  1296.    40 IF (IERR.EQ.0) GO TO 50
  1297.       IDUM(1,1) = IERR
  1298.       CALL PRINTM(19, DUM, ID, ID, ID, IDUM, 1, 1, 1)
  1299.       GO TO 110
  1300.    50 IF (IOPT.EQ.1) GO TO 110
  1301. C ***
  1302. C     COMPUTE EIGENVALUESS OF PSI.
  1303. C ***
  1304.       DO 60 I=1,P1
  1305.            SIG(I,I) = 1.D0
  1306.    60 CONTINUE
  1307.       IF (IPRINT.EQ.4) CALL PRINTM(21, SIG, LDS, P1, P1, IDUM,
  1308.      *     ID, ID, ID)
  1309.       CALL RS(LDS, P1, SIG, W, MATZ, DUM, TEM1, TEM2, IERR)
  1310.       IF (IERR.EQ.0) GO TO 70
  1311.       IDUM(1,1) = IERR
  1312.       CALL PRINTM(20, DUM, ID, ID, ID, IDUM, 1, 1, 1)
  1313.       GO TO 110
  1314.    70 PSIMIN = W(1)
  1315.       PSIK = W(P1)/W(1)
  1316.       IF (IPRINT.EQ.4) CALL PRINTM(25, W, P1, P1, 1, IDUM, ID,
  1317.      *     ID, ID)
  1318. C ***
  1319. C     COMPUTE EIGENVALUES OF LAMBDA.
  1320. C ***
  1321.       DO 90 I=1,P2
  1322.            IP1 = I + P1
  1323.            DO 80 J=I,P2
  1324.                 JP1 = J + P1
  1325.                 SIG(J,I) = SIG(IP1,JP1)
  1326.                 SIG(I,J) = SIG(J,I)
  1327.    80      CONTINUE
  1328.    90 CONTINUE
  1329.       IF (IPRINT.EQ.4) CALL PRINTM(22, SIG, LDS, P2, P2, IDUM,
  1330.      *     ID, ID, ID)
  1331.       CALL RS(LDS, P2, SIG, W, MATZ, DUM, TEM1, TEM2, IERR)
  1332.       IF (IERR.EQ.0) GO TO 100
  1333.       IDUM(1,1) = IERR
  1334.       CALL PRINTM(20, DUM, ID, ID, ID, IDUM, 1, 1, 1)
  1335.       GO TO 110
  1336.   100 XLAMAX = W(P2)
  1337.       IF (IPRINT.EQ.4) CALL PRINTM(24, W, P2, P2, 1, IDUM, ID,
  1338.      *     ID, ID)
  1339.   110 RETURN
  1340.       END
  1341. C
  1342. C --------------------------------------------------------------
  1343. C
  1344.       SUBROUTINE BETA(LVAR, LDL, NO, XSAVE, N, P, CORR, LDC,
  1345.      *     P2, XMU, X, LDX, V, LDV, W, WK, WK2, PVAR, IPRINT)
  1346. C ***
  1347. C     SUBROUTINE BETA COMPUTES AND OUTPUTS FOR THE SELECTED SETS OF
  1348. C        ESTIMATED VARIABLES THE BETA MATRIX OF COEFFICIENTS, THE BETA
  1349. C        VECTOR OF INTERCEPTS, 95 PER CENT CONFIDENCE INTERVALS, THE
  1350. C        ZERO/NONZERO STRUCTURE OF THE BETAS, AND THE FULL SPACE
  1351. C        RESIDUAL.
  1352. C ***
  1353.       INTEGER I, ID, IE, IERR, IL, IOPT, IP, IPRINT, IS, J,
  1354.      *     JIS, JJ, JL, JN, JP, JP1, K, KK, KP1, KT, L, LDC,
  1355.      *     LDL, LDV, LDX, LEN, N, NO, P, P1, P2
  1356.       INTEGER LSTR(36), LVAR(LDL,NO), PVAR(P), IDUM(1,1)
  1357.       DOUBLE PRECISION CVAR, QUAD, RESID, RN, SRTN, SUM, DUM1,
  1358.      *     DUM2, DUM3
  1359.       DOUBLE PRECISION CORR(LDC,P), V(LDV,P), W(P), WK(N),
  1360.      *     WK2(P), X(LDX,P), XMU(P), XSAVE(N,P), DUM(1,1)
  1361.       DOUBLE PRECISION DABS, DSQRT
  1362.       DATA IZERO /1H0/, IAST /1H*/
  1363.       ID = 1
  1364.       P1 = P - P2
  1365. C ***
  1366.       DO 350 K=1,NO
  1367.            WRITE (6,99999)
  1368. C ***
  1369. C     DETERMINE VARIABLE ORDER.
  1370. C ***
  1371.            JP = 0
  1372.            DO 20 J=1,P
  1373.                 DO 10 JJ=1,P2
  1374.                      IF (J.EQ.LVAR(JJ,K)) GO TO 20
  1375.    10           CONTINUE
  1376.                 JP = JP + 1
  1377.                 PVAR(JP) = J
  1378.    20      CONTINUE
  1379.            DO 30 J=1,P2
  1380.                 JL = P1 + J
  1381.                 JP = LVAR(J,K)
  1382.                 PVAR(JL) = JP
  1383.    30      CONTINUE
  1384. C ***
  1385. C     COPY PERMUTED COVARIANCE MATRIX.
  1386. C ***
  1387.            DO 50 I=1,P
  1388.                 IP = PVAR(I)
  1389.                 DO 40 J=1,P
  1390.                      JP = PVAR(J)
  1391.                      V(I,J) = CORR(IP,JP)
  1392.    40           CONTINUE
  1393.    50      CONTINUE
  1394. C ***
  1395. C     COMPUTE LAMBDA MATRIX DIAGONALS.
  1396. C ***
  1397.            IOPT = 1
  1398.            CALL CHKEV(V, LDV, P, P1, IOPT, W, WK2, WK, DUM1,
  1399.      *          DUM2, DUM3, IPRINT)
  1400.            DO 60 I=1,P2
  1401.                 IP = P1 + I
  1402.                 WK(I) = V(IP,IP)
  1403.    60      CONTINUE
  1404. C ***
  1405. C     COMPUTE BETA BY PSI(INV)*(PSI*BETA).
  1406. C ***
  1407.            DO 80 I=1,P
  1408.                 IP = PVAR(I)
  1409.                 DO 70 J=1,P
  1410.                      JP = PVAR(J)
  1411.                      V(I,J) = CORR(IP,JP)
  1412.    70           CONTINUE
  1413.    80      CONTINUE
  1414.            CALL DPOFA(V, LDV, P1, IERR)
  1415.            IF (IERR.NE.0) WRITE (6,99998) IERR
  1416.            DO 90 J=1,P2
  1417.                 JP = P1 + J
  1418.                 CALL DPOSL(V, LDV, P1, V(1,JP))
  1419.    90      CONTINUE
  1420. C ***
  1421. C     COMPUTE QUADRATIC FORM FOR CONFIDENCE INTERVALS.
  1422. C ***
  1423.            DO 100 I=1,P1
  1424.                 IP = PVAR(I)
  1425.                 W(I) = XMU(IP)
  1426.   100      CONTINUE
  1427.            CALL DPOSL(V, LDV, P1, W)
  1428.            QUAD = 0.D0
  1429.            DO 110 I=1,P1
  1430.                 IP = PVAR(I)
  1431.                 QUAD = QUAD + XMU(IP)*W(I)
  1432.   110      CONTINUE
  1433.            QUAD = QUAD + 1.D0
  1434. C ***
  1435. C     COMPUTE DIAGONALS OF THE PSI INVERSE MATRIX.
  1436. C ***
  1437.            DO 130 I=1,P1
  1438.                 DO 120 J=1,P1
  1439.                      W(J) = 0.D0
  1440.   120           CONTINUE
  1441.                 W(I) = 1.D0
  1442.                 CALL DPOSL(V, LDV, P1, W)
  1443.                 WK2(I) = W(I)
  1444.   130      CONTINUE
  1445. C ***
  1446. C     COMPUTE CONFIDENCE INTERVALS FOR BETA.
  1447. C ***
  1448.            RN = N
  1449.            SRTN = DSQRT(RN)
  1450.            DO 150 I=1,P2
  1451.                 IP = I + P1
  1452.                 DO 140 J=1,P1
  1453.                      V(IP,J) = 1.96D0*DSQRT(WK2(J)*WK(I))/SRTN
  1454.   140           CONTINUE
  1455.   150      CONTINUE
  1456. C ***
  1457. C    COMPUTE CONFIDENCE INTERVALS FOR BETA0.
  1458. C ***
  1459.            DO 160 I=1,P2
  1460.                 IP = I + P1
  1461.                 W(I) = 1.96D0*DSQRT(QUAD*WK(I))/SRTN
  1462.   160      CONTINUE
  1463. C ***
  1464. C    COMPUTE BETA0.
  1465. C ***
  1466.            DO 180 I=1,P2
  1467.                 SUM = 0.D0
  1468.                 IL = I + P1
  1469.                 DO 170 J=1,P1
  1470.                      JP = PVAR(J)
  1471.                      SUM = SUM + V(J,IL)*XMU(JP)
  1472.   170           CONTINUE
  1473.                 IP = PVAR(IL)
  1474.                 WK2(I) = XMU(IP) - SUM
  1475.   180      CONTINUE
  1476. C ***
  1477. C    OUTPUT RESULTS.
  1478. C ***
  1479.            WRITE (6,99997)
  1480.            KT = (P1+7)/8
  1481.            DO 220 KK=1,KT
  1482.                 IS = (KK-1)*8 + 1
  1483.                 IE = KK*8
  1484.                 IF (IE.GT.P1) IE = P1
  1485.                 IF (KK.GT.1) GO TO 200
  1486.                 WRITE (6,99996) (PVAR(I),I=IS,IE)
  1487.                 DO 190 I=1,P2
  1488.                      J = I + P1
  1489.                      WRITE (6,99995) PVAR(J), WK2(I),
  1490.      *                    (V(L,J),L=IS,IE)
  1491.                      WRITE (6,99994) W(I), (V(J,L),L=IS,IE)
  1492.   190           CONTINUE
  1493.                 GO TO 220
  1494.   200           WRITE (6,99993) (PVAR(I),I=IS,IE)
  1495.                 DO 210 I=1,P2
  1496.                      J = I + P1
  1497.                      WRITE (6,99992) PVAR(J), (V(L,J),L=IS,IE)
  1498.                      WRITE (6,99991) (V(J,L),L=IS,IE)
  1499.   210           CONTINUE
  1500.   220      CONTINUE
  1501.            WRITE (6,99990)
  1502.            KT = (P1+34)/35
  1503.            DO 280 KK=1,KT
  1504.                 IS = (KK-1)*35 + 1
  1505.                 IE = KK*35
  1506.                 IF (IE.GT.P1) IE = P1
  1507.                 IF (KK.GT.1) GO TO 250
  1508.                 WRITE (6,99989) (PVAR(I),I=IS,IE)
  1509.                 LEN = IE
  1510.                 KP1 = LEN + 1
  1511.                 DO 240 I=1,P2
  1512.                      IP = I + P1
  1513.                      LSTR(1) = IZERO
  1514.                      IF (DABS(WK2(I)).GT.W(I)) LSTR(1) = IAST
  1515.                      DO 230 J=1,LEN
  1516.                           JP1 = J + 1
  1517.                           LSTR(JP1) = IZERO
  1518.                           IF (DABS(V(J,IP)).GT.V(IP,J))
  1519.      *                         LSTR(JP1) = IAST
  1520.   230                CONTINUE
  1521.                      WRITE (6,99988) PVAR(IP), (LSTR(J),J=1,KP1)
  1522.   240           CONTINUE
  1523.                 GO TO 280
  1524.   250           WRITE (6,99987) (PVAR(I),I=IS,IE)
  1525.                 LEN = IE - IS + 1
  1526.                 DO 270 I=1,P2
  1527.                      IP = I + P1
  1528.                      DO 260 J=1,LEN
  1529.                           JIS = J + IS - 1
  1530.                           LSTR(J) = IZERO
  1531.                           IF (DABS(V(JIS,IP)).GT.V(IP,JIS))
  1532.      *                         LSTR(J) = IAST
  1533.   260                CONTINUE
  1534.                      WRITE (6,99986) PVAR(IP), (LSTR(J),J=1,LEN)
  1535.   270           CONTINUE
  1536.   280      CONTINUE
  1537. C ***
  1538. C     COMPUTE CENTERED AND SCALED BETA MATRIX.
  1539. C ***
  1540.            DO 300 I=1,P1
  1541.                 IP = PVAR(I)
  1542.                 CVAR = DSQRT(CORR(IP,IP))
  1543.                 DO 290 J=1,P2
  1544.                      JN = P1 + J
  1545.                      JP = PVAR(JN)
  1546.                      V(I,JN) = V(I,JN)*(CVAR/DSQRT(CORR(JP,JP)))
  1547.   290           CONTINUE
  1548.   300      CONTINUE
  1549. C ***
  1550. C     PERMUTE COLUMNS OF CENTERED AND SCALED X.
  1551. C ***
  1552.            DO 310 J=1,P
  1553.                 JP = PVAR(J)
  1554.                 CALL DCOPY(N, XSAVE(1,JP), 1, X(1,J), 1)
  1555.   310      CONTINUE
  1556.            IDUM(1,1) = 2
  1557.            IF (IPRINT.EQ.4) CALL PRINTM(9, X, LDX, N, P, IDUM,
  1558.      *          1, 1, 1)
  1559. C ***
  1560. C     COMPUTE RESIDUAL (FULL SPACE RESIDUAL).
  1561. C ***
  1562.            RESID = 0.D0
  1563.            DO 340 I=1,N
  1564.                 DO 330 J=1,P2
  1565.                      JP = J + P1
  1566.                      SUM = 0.D0
  1567.                      DO 320 L=1,P1
  1568.                           SUM = SUM + X(I,L)*V(L,JP)
  1569.   320                CONTINUE
  1570.                      X(I,JP) = X(I,JP) - SUM
  1571.                      RESID = RESID + X(I,JP)**2
  1572.   330           CONTINUE
  1573.   340      CONTINUE
  1574.            IDUM(1,1) = 3
  1575.            IF (IPRINT.EQ.4) CALL PRINTM(9, X, LDX, N, P, IDUM,
  1576.      *          1, P1, P2)
  1577.            RESID = DSQRT(RESID)
  1578.            DUM(1,1) = RESID
  1579.            CALL PRINTM(31, DUM, 1, 1, 1, IDUM, ID, ID, ID)
  1580.   350 CONTINUE
  1581.       RETURN
  1582. 99999 FORMAT (///1X, 120(1H-))
  1583. 99998 FORMAT (27H IERR FROM DPOFA IN BETA IS, I3)
  1584. 99997 FORMAT (///43H BETA0 - BETA MATRIX / 95 PER CENT CONFIDEN,
  1585.      *     3HCE , 9HINTERVALS//)
  1586. 99996 FORMAT (25X, 9HPREDICTOR/1X, 9HESTIMATED, 15X, 8HVARIABLE,
  1587.      *     5HS..../1X, 9HVARIABLES, 3X, 9HINTERCEPT, 1X, 8(6X,
  1588.      *     I3, 4X))
  1589. 99995 FORMAT (/4X, I3, 3X, 1P9D13.4)
  1590. 99994 FORMAT (10X, 1P9D13.4)
  1591. 99993 FORMAT (//25X, 9HPREDICTOR/1X, 9HESTIMATED, 15X, 6HVARIAB,
  1592.      *     3HLES, 2X, 14H(CONTINUED).../1X, 9HVARIABLES, 13X,
  1593.      *     8(6X, I3, 4X))
  1594. 99992 FORMAT (/4X, I3, 16X, 1P8D13.4)
  1595. 99991 FORMAT (23X, 1P8D13.4)
  1596. 99990 FORMAT (///1X, 39HDEPENDENCY STRUCTURE AT 95 PER CENT CON,
  1597.      *     8HFIDENCE , 5HLEVEL/11X, 24H(* = NONZERO,  0 = ZERO)/
  1598.      *     )
  1599. 99989 FORMAT (13X, 9HPREDICTOR/1X, 9HESTIMATED, 3X, 9HVARIABLES,
  1600.      *     4H..../1X, 9HVARIABLES, 3H  I, 35(I3))
  1601. 99988 FORMAT (/4X, I3, 3X, 36(2X, A1))
  1602. 99987 FORMAT (13X, 9HPREDICTOR/1X, 9HESTIMATED, 3X, 9HVARIABLES,
  1603.      *     1H , 14H(CONTINUED).../1X, 9HVARIABLES, 3X, 35I3)
  1604. 99986 FORMAT (/10X, 35(2X, A1))
  1605.       END
  1606. C
  1607. C --------------------------------------------------------------
  1608. C
  1609.       SUBROUTINE SCPROD(X, N, PROD, ISC)
  1610. C ***
  1611. C     SUBROUTINE SCPROD COMPUTES THE PRODUCT OF A SEQUENCE OF NUMBERS
  1612. C        USING A SCALING PROCEDURE WHICH COMPUTES THE ANSWER AS A
  1613. C        MANTISSA AND EXPONENT.
  1614. C ***
  1615.       INTEGER I, ISC, N
  1616.       DOUBLE PRECISION PROD, SC, X(N)
  1617.       DOUBLE PRECISION DABS
  1618.       PROD = 1.D0
  1619.       ISC = 0
  1620.       SC = 10.D0
  1621.       DO 30 I=1,N
  1622.            PROD = X(I)*PROD
  1623.            IF (PROD.EQ.0.D0) GO TO 40
  1624.    10      IF (DABS(PROD).GE.1.D0) GO TO 20
  1625.            PROD = SC*PROD
  1626.            ISC = ISC - 1
  1627.            GO TO 10
  1628.    20      IF (DABS(PROD).LT.SC) GO TO 30
  1629.            PROD = PROD/SC
  1630.            ISC = ISC + 1
  1631.            GO TO 20
  1632.    30 CONTINUE
  1633.    40 CONTINUE
  1634.       RETURN
  1635.       END
  1636.       SUBROUTINE DPODI(A,LDA,N,DET,JOB)                                 
  1637. 00000010
  1638.       INTEGER LDA,N,JOB
  1639.       DOUBLE PRECISION A(LDA,1)
  1640.       DOUBLE PRECISION DET(2)
  1641. C
  1642. C     DPODI COMPUTES THE DETERMINANT AND INVERSE OF A CERTAIN
  1643. C     DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE MATRIX (SEE BELOW)
  1644. C     USING THE FACTORS COMPUTED BY DPOCO, DPOFA OR DQRDC.
  1645. C
  1646. C     ON ENTRY
  1647. C
  1648. C        A       DOUBLE PRECISION(LDA, N)
  1649. C                THE OUTPUT  A  FROM DPOCO OR DPOFA
  1650. C                OR THE OUTPUT  X  FROM DQRDC.
  1651. C
  1652. C        LDA     INTEGER
  1653. C                THE LEADING DIMENSION OF THE ARRAY  A .
  1654. C
  1655. C        N       INTEGER
  1656. C                THE ORDER OF THE MATRIX  A .
  1657. C
  1658. C        JOB     INTEGER
  1659. C                = 11   BOTH DETERMINANT AND INVERSE.
  1660. C                = 01   INVERSE ONLY.
  1661. C                = 10   DETERMINANT ONLY.
  1662. C
  1663. C     ON RETURN
  1664. C
  1665. C        A       IF DPOCO OR DPOFA WAS USED TO FACTOR  A  THEN
  1666. C                DPODI PRODUCES THE UPPER HALF OF INVERSE(A) .
  1667. C                IF DQRDC WAS USED TO DECOMPOSE  X  THEN
  1668. C                DPODI PRODUCES THE UPPER HALF OF INVERSE(TRANS(X)*X)
  1669. C                WHERE TRANS(X) IS THE TRANSPOSE.
  1670. C                ELEMENTS OF  A  BELOW THE DIAGONAL ARE UNCHANGED.
  1671. C                IF THE UNITS DIGIT OF JOB IS ZERO,  A  IS UNCHANGED.
  1672. C
  1673. C        DET     DOUBLE PRECISION(2)
  1674. C                DETERMINANT OF  A  OR OF  TRANS(X)*X  IF REQUESTED.
  1675. C                OTHERWISE NOT REFERENCED.
  1676. C                DETERMINANT = DET(1) * 10.0**DET(2)
  1677. C                WITH  1.0 .LE. DET(1) .LT. 10.0
  1678. C                OR  DET(1) .EQ. 0.0 .
  1679. C
  1680. C     ERROR CONDITION
  1681. C
  1682. C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
  1683. C        A ZERO ON THE DIAGONAL AND THE INVERSE IS REQUESTED.
  1684. C        IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED CORRECTLY
  1685. C        AND IF DPOCO OR DPOFA HAS SET INFO .EQ. 0 .
  1686. C
  1687. C     LINPACK.  THIS VERSION DATED 08/14/78 .
  1688. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  1689. C
  1690. C     SUBROUTINES AND FUNCTIONS
  1691. C
  1692. C     BLAS DAXPY,DSCAL
  1693. C     FORTRAN MOD
  1694. C
  1695. C     INTERNAL VARIABLES
  1696. C
  1697.       DOUBLE PRECISION T
  1698.       DOUBLE PRECISION S
  1699.       INTEGER I,J,JM1,K,KP1
  1700. C
  1701. C     COMPUTE DETERMINANT
  1702. C
  1703.       IF (JOB/10 .EQ. 0) GO TO 70
  1704.          DET(1) = 1.0D0
  1705.          DET(2) = 0.0D0
  1706.          S = 10.0D0
  1707.          DO 50 I = 1, N
  1708.             DET(1) = A(I,I)**2*DET(1)
  1709. C        ...EXIT
  1710.             IF (DET(1) .EQ. 0.0D0) GO TO 60
  1711.    10       IF (DET(1) .GE. 1.0D0) GO TO 20
  1712.                DET(1) = S*DET(1)
  1713.                DET(2) = DET(2) - 1.0D0
  1714.             GO TO 10
  1715.    20       CONTINUE
  1716.    30       IF (DET(1) .LT. S) GO TO 40
  1717.                DET(1) = DET(1)/S
  1718.                DET(2) = DET(2) + 1.0D0
  1719.             GO TO 30
  1720.    40       CONTINUE
  1721.    50    CONTINUE
  1722.    60    CONTINUE
  1723.    70 CONTINUE
  1724. C
  1725. C     COMPUTE INVERSE(R)
  1726. C
  1727.       IF (MOD(JOB,10) .EQ. 0) GO TO 140
  1728.          DO 100 K = 1, N
  1729.             A(K,K) = 1.0D0/A(K,K)
  1730.             T = -A(K,K)
  1731.             CALL DSCAL(K-1,T,A(1,K),1)
  1732.             KP1 = K + 1
  1733.             IF (N .LT. KP1) GO TO 90
  1734.             DO 80 J = KP1, N
  1735.                T = A(K,J)
  1736.                A(K,J) = 0.0D0
  1737.                CALL DAXPY(K,T,A(1,K),1,A(1,J),1)
  1738.    80       CONTINUE
  1739.    90       CONTINUE
  1740.   100    CONTINUE
  1741. C
  1742. C        FORM  INVERSE(R) * TRANS(INVERSE(R))
  1743. C
  1744.          DO 130 J = 1, N
  1745.             JM1 = J - 1
  1746.             IF (JM1 .LT. 1) GO TO 120
  1747.             DO 110 K = 1, JM1
  1748.                T = A(K,J)
  1749.                CALL DAXPY(K,T,A(1,J),1,A(1,K),1)
  1750.   110       CONTINUE
  1751.   120       CONTINUE
  1752.             T = A(J,J)
  1753.             CALL DSCAL(J,T,A(1,J),1)
  1754.   130    CONTINUE
  1755.   140 CONTINUE
  1756.       RETURN
  1757.       END
  1758.       SUBROUTINE DPOFA(A,LDA,N,INFO)                                    
  1759. 00000010
  1760.       INTEGER LDA,N,INFO
  1761.       DOUBLE PRECISION A(LDA,1)
  1762. C
  1763. C     DPOFA FACTORS A DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE
  1764. C     MATRIX.
  1765. C
  1766. C     DPOFA IS USUALLY CALLED BY DPOCO, BUT IT CAN BE CALLED
  1767. C     DIRECTLY WITH A SAVING IN TIME IF  RCOND  IS NOT NEEDED.
  1768. C     (TIME FOR DPOCO) = (1 + 18/N)*(TIME FOR DPOFA) .
  1769. C
  1770. C     ON ENTRY
  1771. C
  1772. C        A       DOUBLE PRECISION(LDA, N)
  1773. C                THE SYMMETRIC MATRIX TO BE FACTORED.  ONLY THE
  1774. C                DIAGONAL AND UPPER TRIANGLE ARE USED.
  1775. C
  1776. C        LDA     INTEGER
  1777. C                THE LEADING DIMENSION OF THE ARRAY  A .
  1778. C
  1779. C        N       INTEGER
  1780. C                THE ORDER OF THE MATRIX  A .
  1781. C
  1782. C     ON RETURN
  1783. C
  1784. C        A       AN UPPER TRIANGULAR MATRIX  R  SO THAT  A = TRANS(R)*R
  1785. C                WHERE  TRANS(R)  IS THE TRANSPOSE.
  1786. C                THE STRICT LOWER TRIANGLE IS UNALTERED.
  1787. C                IF  INFO .NE. 0 , THE FACTORIZATION IS NOT COMPLETE.
  1788. C
  1789. C        INFO    INTEGER
  1790. C                = 0  FOR NORMAL RETURN.
  1791. C                = K  SIGNALS AN ERROR CONDITION.  THE LEADING MINOR
  1792. C                     OF ORDER  K  IS NOT POSITIVE DEFINITE.
  1793. C
  1794. C     LINPACK.  THIS VERSION DATED 08/14/78 .
  1795. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  1796. C
  1797. C     SUBROUTINES AND FUNCTIONS
  1798. C
  1799. C     BLAS DDOT
  1800. C     FORTRAN DSQRT
  1801. C
  1802. C     INTERNAL VARIABLES
  1803. C
  1804.       DOUBLE PRECISION DDOT,T
  1805.       DOUBLE PRECISION S
  1806.       INTEGER J,JM1,K
  1807. C     BEGIN BLOCK WITH ...EXITS TO 40
  1808. C
  1809. C
  1810.          DO 30 J = 1, N
  1811.             INFO = J
  1812.             S = 0.0D0
  1813.             JM1 = J - 1
  1814.             IF (JM1 .LT. 1) GO TO 20
  1815.             DO 10 K = 1, JM1
  1816.                T = A(K,J) - DDOT(K-1,A(1,K),1,A(1,J),1)
  1817.                T = T/A(K,K)
  1818.                A(K,J) = T
  1819.                S = S + T*T
  1820.    10       CONTINUE
  1821.    20       CONTINUE
  1822.             S = A(J,J) - S
  1823. C     ......EXIT
  1824.             IF (S .LE. 0.0D0) GO TO 40
  1825.             A(J,J) = DSQRT(S)
  1826.    30    CONTINUE
  1827.          INFO = 0
  1828.    40 CONTINUE
  1829.       RETURN
  1830.       END
  1831.       SUBROUTINE DPOSL(A,LDA,N,B)                                       
  1832. 00000010
  1833.       INTEGER LDA,N
  1834.       DOUBLE PRECISION A(LDA,1),B(1)
  1835. C
  1836. C     DPOSL SOLVES THE DOUBLE PRECISION SYMMETRIC POSITIVE DEFINITE
  1837. C     SYSTEM A * X = B
  1838. C     USING THE FACTORS COMPUTED BY DPOCO OR DPOFA.
  1839. C
  1840. C     ON ENTRY
  1841. C
  1842. C        A       DOUBLE PRECISION(LDA, N)
  1843. C                THE OUTPUT FROM DPOCO OR DPOFA.
  1844. C
  1845. C        LDA     INTEGER
  1846. C                THE LEADING DIMENSION OF THE ARRAY  A .
  1847. C
  1848. C        N       INTEGER
  1849. C                THE ORDER OF THE MATRIX  A .
  1850. C
  1851. C        B       DOUBLE PRECISION(N)
  1852. C                THE RIGHT HAND SIDE VECTOR.
  1853. C
  1854. C     ON RETURN
  1855. C
  1856. C        B       THE SOLUTION VECTOR  X .
  1857. C
  1858. C     ERROR CONDITION
  1859. C
  1860. C        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS
  1861. C        A ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES
  1862. C        SINGULARITY BUT IT IS USUALLY CAUSED BY IMPROPER SUBROUTINE
  1863. C        ARGUMENTS.  IT WILL NOT OCCUR IF THE SUBROUTINES ARE CALLED
  1864. C        CORRECTLY AND  INFO .EQ. 0 .
  1865. C
  1866. C     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
  1867. C     WITH  P  COLUMNS
  1868. C           CALL DPOCO(A,LDA,N,RCOND,Z,INFO)
  1869. C           IF (RCOND IS TOO SMALL .OR. INFO .NE. 0) GO TO ...
  1870. C           DO 10 J = 1, P
  1871. C              CALL DPOSL(A,LDA,N,C(1,J))
  1872. C        10 CONTINUE
  1873. C
  1874. C     LINPACK.  THIS VERSION DATED 08/14/78 .
  1875. C     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
  1876. C
  1877. C     SUBROUTINES AND FUNCTIONS
  1878. C
  1879. C     BLAS DAXPY,DDOT
  1880. C
  1881. C     INTERNAL VARIABLES
  1882. C
  1883.       DOUBLE PRECISION DDOT,T
  1884.       INTEGER K,KB
  1885. C
  1886. C     SOLVE TRANS(R)*Y = B
  1887. C
  1888.       DO 10 K = 1, N
  1889.          T = DDOT(K-1,A(1,K),1,B(1),1)
  1890.          B(K) = (B(K) - T)/A(K,K)
  1891.    10 CONTINUE
  1892. C
  1893. C     SOLVE R*X = Y
  1894. C
  1895.       DO 20 KB = 1, N
  1896.          K = N + 1 - KB
  1897.          B(K) = B(K)/A(K,K)
  1898.          T = -B(K)
  1899.          CALL DAXPY(K-1,T,A(1,K),1,B(1),1)
  1900.    20 CONTINUE
  1901.       RETURN
  1902.       END
  1903.       SUBROUTINE DSVDC(X,LDX,N,P,S,E,U,LDU,V,LDV,WORK,JOB,INFO)         
  1904. 00000010
  1905.       INTEGER LDX,N,P,LDU,LDV,JOB,INFO
  1906.       DOUBLE PRECISION X(LDX,1),S(1),E(1),U(LDU,1),V(LDV,1),WORK(1)
  1907. C
  1908. C
  1909. C     DSVDC IS A SUBROUTINE TO REDUCE A DOUBLE PRECISION NXP MATRIX X
  1910. C     BY ORTHOGONAL TRANSFORMATIONS U AND V TO DIAGONAL FORM.  THE
  1911. C     DIAGONAL ELEMENTS S(I) ARE THE SINGULAR VALUES OF X.  THE
  1912. C     COLUMNS OF U ARE THE CORRESPONDING LEFT SINGULAR VECTORS,
  1913. C     AND THE COLUMNS OF V THE RIGHT SINGULAR VECTORS.
  1914. C
  1915. C     ON ENTRY
  1916. C
  1917. C         X         DOUBLE PRECISION(LDX,P), WHERE LDX.GE.N.
  1918. C                   X CONTAINS THE MATRIX WHOSE SINGULAR VALUE
  1919. C                   DECOMPOSITION IS TO BE COMPUTED.  X IS
  1920. C                   DESTROYED BY DSVDC.
  1921. C
  1922. C         LDX       INTEGER.
  1923. C                   LDX IS THE LEADING DIMENSION OF THE ARRAY X.
  1924. C
  1925. C         N         INTEGER.
  1926. C                   N IS THE NUMBER OF ROWS OF THE MATRIX X.
  1927. C
  1928. C         P         INTEGER.
  1929. C                   P IS THE NUMBER OF COLUMNS OF THE MATRIX X.
  1930. C
  1931. C         LDU       INTEGER.
  1932. C                   LDU IS THE LEADING DIMENSION OF THE ARRAY U.
  1933. C                   (SEE BELOW).
  1934. C
  1935. C         LDV       INTEGER.
  1936. C                   LDV IS THE LEADING DIMENSION OF THE ARRAY V.
  1937. C                   (SEE BELOW).
  1938. C
  1939. C         WORK      DOUBLE PRECISION(N).
  1940. C                   WORK IS A SCRATCH ARRAY.
  1941. C
  1942. C         JOB       INTEGER.
  1943. C                   JOB CONTROLS THE COMPUTATION OF THE SINGULAR
  1944. C                   VECTORS.  IT HAS THE DECIMAL EXPANSION AB
  1945. C                   WITH THE FOLLOWING MEANING
  1946. C
  1947. C                        A.EQ.0    DO NOT COMPUTE THE LEFT SINGULAR
  1948. C                                  VECTORS.
  1949. C                        A.EQ.1    RETURN THE N LEFT SINGULAR VECTORS
  1950. C                                  IN U.
  1951. C                        A.GE.2    RETURN THE FIRST MIN(N,P) SINGULAR
  1952. C                                  VECTORS IN U.
  1953. C                        B.EQ.0    DO NOT COMPUTE THE RIGHT SINGULAR
  1954. C                                  VECTORS.
  1955. C                        B.EQ.1    RETURN THE RIGHT SINGULAR VECTORS
  1956. C                                  IN V.
  1957. C
  1958. C     ON RETURN
  1959. C
  1960. C         S         DOUBLE PRECISION(MM), WHERE MM=MIN(N+1,P).
  1961. C                   THE FIRST MIN(N,P) ENTRIES OF S CONTAIN THE
  1962. C                   SINGULAR VALUES OF X ARRANGED IN DESCENDING
  1963. C                   ORDER OF MAGNITUDE.
  1964. C
  1965. C         E         DOUBLE PRECISION(P)
  1966. C                   E ORDINARILY CONTAINS ZEROS.  HOWEVER SEE THE
  1967. C                   DISCUSSION OF INFO FOR EXCEPTIONS.
  1968. C
  1969. C         U         DOUBLE PRECISION(LDU,K), WHERE LDU.GE.N.  IF
  1970. C                                   JOBA.EQ.1 THEN K.EQ.N, IF JOBA.GE.2
  1971. C                                   THEN K.EQ.MIN(N,P).
  1972. C                   U CONTAINS THE MATRIX OF LEFT SINGULAR VECTORS.
  1973. C                   U IS NOT REFERENCED IF JOBA.EQ.0.  IF N.LE.P
  1974. C                   OR IF JOBA.EQ.2, THEN U MAY BE IDENTIFIED WITH X
  1975. C                   IN THE SUBROUTINE CALL.
  1976. C
  1977. C         V         DOUBLE PRECISION(LDV,P), WHERE LDV.GE.P.
  1978. C                   V CONTAINS THE MATRIX OF RIGHT SINGULAR VECTORS.
  1979. C                   V IS NOT REFERENCED IF JOB.EQ.0.  IF P.LE.N,
  1980. C                   THEN V MAY BE IDENTIFIED WITH X IN THE
  1981. C                   SUBROUTINE CALL.
  1982. C
  1983. C         INFO      INTEGER.
  1984. C                   THE SINGULAR VALUES (AND THEIR CORRESPONDING
  1985. C                   SINGULAR VECTORS) S(INFO+1),S(INFO+2),...,S(M)
  1986. C                   ARE CORRECT (HERE M=MIN(N,P)).  THUS IF
  1987. C                   INFO.EQ.0, ALL THE SINGULAR VALUES AND THEIR
  1988. C                   VECTORS ARE CORRECT.  IN ANY EVENT, THE MATRIX
  1989. C                   B = TRANS(U)*X*V IS THE BIDIAGONAL MATRIX
  1990. C                   WITH THE ELEMENTS OF S ON ITS DIAGONAL AND THE
  1991. C                   ELEMENTS OF E ON ITS SUPER-DIAGONAL (TRANS(U)
  1992. C                   IS THE TRANSPOSE OF U).  THUS THE SINGULAR
  1993. C                   VALUES OF X AND B ARE THE SAME.
  1994. C
  1995. C     LINPACK. THIS VERSION DATED 08/14/78 .
  1996. C     G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB.
  1997. C
  1998. C     DSVDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS.
  1999. C
  2000. C     EXTERNAL DROT
  2001. C     BLAS DAXPY,DDOT,DSCAL,DSWAP,DNRM2,DROTG
  2002. C     FORTRAN DABS,DMAX1,MAX0,MIN0,MOD,DSQRT
  2003. C
  2004. C     INTERNAL VARIABLES
  2005. C
  2006.       INTEGER I,ITER,J,JOBU,K,KASE,KK,L,LL,LLS,LM1,LP1,LS,LU,M,MAXIT,
  2007.      *        MM,MM1,MP1,NCT,NCTP1,NCU,NRT,NRTP1
  2008.       DOUBLE PRECISION DDOT,T
  2009.       DOUBLE PRECISION B,C,CS,EL,EMM1,F,G,DNRM2,SCALE,SHIFT,SL,SM,SN,
  2010.      *                 SMM1,T1,TEST,ZTEST
  2011.       LOGICAL WANTU,WANTV
  2012. C
  2013. C
  2014. C     SET THE MAXIMUM NUMBER OF ITERATIONS.
  2015. C
  2016.       MAXIT = 30
  2017. C
  2018. C     DETERMINE WHAT IS TO BE COMPUTED.
  2019. C
  2020.       WANTU = .FALSE.
  2021.       WANTV = .FALSE.
  2022.       JOBU = MOD(JOB,100)/10
  2023.       NCU = N
  2024.       IF (JOBU .GT. 1) NCU = MIN0(N,P)
  2025.       IF (JOBU .NE. 0) WANTU = .TRUE.
  2026.       IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE.
  2027. C
  2028. C     REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS
  2029. C     IN S AND THE SUPER-DIAGONAL ELEMENTS IN E.
  2030. C
  2031.       INFO = 0
  2032.       NCT = MIN0(N-1,P)
  2033.       NRT = MAX0(0,MIN0(P-2,N))
  2034.       LU = MAX0(NCT,NRT)
  2035.       IF (LU .LT. 1) GO TO 170
  2036.       DO 160 L = 1, LU
  2037.          LP1 = L + 1
  2038.          IF (L .GT. NCT) GO TO 20
  2039. C
  2040. C           COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND
  2041. C           PLACE THE L-TH DIAGONAL IN S(L).
  2042. C
  2043.             S(L) = DNRM2(N-L+1,X(L,L),1)
  2044.             IF (S(L) .EQ. 0.0D0) GO TO 10
  2045.                IF (X(L,L) .NE. 0.0D0) S(L) = DSIGN(S(L),X(L,L))
  2046.                CALL DSCAL(N-L+1,1.0D0/S(L),X(L,L),1)
  2047.                X(L,L) = 1.0D0 + X(L,L)
  2048.    10       CONTINUE
  2049.             S(L) = -S(L)
  2050.    20    CONTINUE
  2051.          IF (P .LT. LP1) GO TO 50
  2052.          DO 40 J = LP1, P
  2053.             IF (L .GT. NCT) GO TO 30
  2054.             IF (S(L) .EQ. 0.0D0) GO TO 30
  2055. C
  2056. C              APPLY THE TRANSFORMATION.
  2057. C
  2058.                T = -DDOT(N-L+1,X(L,L),1,X(L,J),1)/X(L,L)
  2059.                CALL DAXPY(N-L+1,T,X(L,L),1,X(L,J),1)
  2060.    30       CONTINUE
  2061. C
  2062. C           PLACE THE L-TH ROW OF X INTO  E FOR THE
  2063. C           SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION.
  2064. C
  2065.             E(J) = X(L,J)
  2066.    40    CONTINUE
  2067.    50    CONTINUE
  2068.          IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 70
  2069. C
  2070. C           PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK
  2071. C           MULTIPLICATION.
  2072. C
  2073.             DO 60 I = L, N
  2074.                U(I,L) = X(I,L)
  2075.    60       CONTINUE
  2076.    70    CONTINUE
  2077.          IF (L .GT. NRT) GO TO 150
  2078. C
  2079. C           COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE
  2080. C           L-TH SUPER-DIAGONAL IN E(L).
  2081. C
  2082.             E(L) = DNRM2(P-L,E(LP1),1)
  2083.             IF (E(L) .EQ. 0.0D0) GO TO 80
  2084.                IF (E(LP1) .NE. 0.0D0) E(L) = DSIGN(E(L),E(LP1))
  2085.                CALL DSCAL(P-L,1.0D0/E(L),E(LP1),1)
  2086.                E(LP1) = 1.0D0 + E(LP1)
  2087.    80       CONTINUE
  2088.             E(L) = -E(L)
  2089.             IF (LP1 .GT. N .OR. E(L) .EQ. 0.0D0) GO TO 120
  2090. C
  2091. C              APPLY THE TRANSFORMATION.
  2092. C
  2093.                DO 90 I = LP1, N
  2094.                   WORK(I) = 0.0D0
  2095.    90          CONTINUE
  2096.                DO 100 J = LP1, P
  2097.                   CALL DAXPY(N-L,E(J),X(LP1,J),1,WORK(LP1),1)
  2098.   100          CONTINUE
  2099.                DO 110 J = LP1, P
  2100.                   CALL DAXPY(N-L,-E(J)/E(LP1),WORK(LP1),1,X(LP1,J),1)
  2101.   110          CONTINUE
  2102.   120       CONTINUE
  2103.             IF (.NOT.WANTV) GO TO 140
  2104. C
  2105. C              PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT
  2106. C              BACK MULTIPLICATION.
  2107. C
  2108.                DO 130 I = LP1, P
  2109.                   V(I,L) = E(I)
  2110.   130          CONTINUE
  2111.   140       CONTINUE
  2112.   150    CONTINUE
  2113.   160 CONTINUE
  2114.   170 CONTINUE
  2115. C
  2116. C     SET UP THE FINAL BIDIAGONAL MATRIX OF ORDER M.
  2117. C
  2118.       M = MIN0(P,N+1)
  2119.       NCTP1 = NCT + 1
  2120.       NRTP1 = NRT + 1
  2121.       IF (NCT .LT. P) S(NCTP1) = X(NCTP1,NCTP1)
  2122.       IF (N .LT. M) S(M) = 0.0D0
  2123.       IF (NRTP1 .LT. M) E(NRTP1) = X(NRTP1,M)
  2124.       E(M) = 0.0D0
  2125. C
  2126. C     IF REQUIRED, GENERATE U.
  2127. C
  2128.       IF (.NOT.WANTU) GO TO 300
  2129.          IF (NCU .LT. NCTP1) GO TO 200
  2130.          DO 190 J = NCTP1, NCU
  2131.             DO 180 I = 1, N
  2132.                U(I,J) = 0.0D0
  2133.   180       CONTINUE
  2134.             U(J,J) = 1.0D0
  2135.   190    CONTINUE
  2136.   200    CONTINUE
  2137.          IF (NCT .LT. 1) GO TO 290
  2138.          DO 280 LL = 1, NCT
  2139.             L = NCT - LL + 1
  2140.             IF (S(L) .EQ. 0.0D0) GO TO 250
  2141.                LP1 = L + 1
  2142.                IF (NCU .LT. LP1) GO TO 220
  2143.                DO 210 J = LP1, NCU
  2144.                   T = -DDOT(N-L+1,U(L,L),1,U(L,J),1)/U(L,L)
  2145.                   CALL DAXPY(N-L+1,T,U(L,L),1,U(L,J),1)
  2146.   210          CONTINUE
  2147.   220          CONTINUE
  2148.                CALL DSCAL(N-L+1,-1.0D0,U(L,L),1)
  2149.                U(L,L) = 1.0D0 + U(L,L)
  2150.                LM1 = L - 1
  2151.                IF (LM1 .LT. 1) GO TO 240
  2152.                DO 230 I = 1, LM1
  2153.                   U(I,L) = 0.0D0
  2154.   230          CONTINUE
  2155.   240          CONTINUE
  2156.             GO TO 270
  2157.   250       CONTINUE
  2158.                DO 260 I = 1, N
  2159.                   U(I,L) = 0.0D0
  2160.   260          CONTINUE
  2161.                U(L,L) = 1.0D0
  2162.   270       CONTINUE
  2163.   280    CONTINUE
  2164.   290    CONTINUE
  2165.   300 CONTINUE
  2166. C
  2167. C     IF IT IS REQUIRED, GENERATE V.
  2168. C
  2169.       IF (.NOT.WANTV) GO TO 350
  2170.          DO 340 LL = 1, P
  2171.             L = P - LL + 1
  2172.             LP1 = L + 1
  2173.             IF (L .GT. NRT) GO TO 320
  2174.             IF (E(L) .EQ. 0.0D0) GO TO 320
  2175.                DO 310 J = LP1, P
  2176.                   T = -DDOT(P-L,V(LP1,L),1,V(LP1,J),1)/V(LP1,L)
  2177.                   CALL DAXPY(P-L,T,V(LP1,L),1,V(LP1,J),1)
  2178.   310          CONTINUE
  2179.   320       CONTINUE
  2180.             DO 330 I = 1, P
  2181.                V(I,L) = 0.0D0
  2182.   330       CONTINUE
  2183.             V(L,L) = 1.0D0
  2184.   340    CONTINUE
  2185.   350 CONTINUE
  2186. C
  2187. C     MAIN ITERATION LOOP FOR THE SINGULAR VALUES.
  2188. C
  2189.       MM = M
  2190.       ITER = 0
  2191.   360 CONTINUE
  2192. C
  2193. C        QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND.
  2194. C
  2195. C     ...EXIT
  2196.          IF (M .EQ. 0) GO TO 620
  2197. C
  2198. C        IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET
  2199. C        FLAG AND RETURN.
  2200. C
  2201.          IF (ITER .LT. MAXIT) GO TO 370
  2202.             INFO = M
  2203. C     ......EXIT
  2204.             GO TO 620
  2205.   370    CONTINUE
  2206. C
  2207. C        THIS SECTION OF THE PROGRAM INSPECTS FOR
  2208. C        NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS.  ON
  2209. C        COMPLETION THE VARIABLES KASE AND L ARE SET AS FOLLOWS.
  2210. C
  2211. C           KASE = 1     IF S(M) AND E(L-1) ARE NEGLIGIBLE AND L.LT.M
  2212. C           KASE = 2     IF S(L) IS NEGLIGIBLE AND L.LT.M
  2213. C           KASE = 3     IF E(L-1) IS NEGLIGIBLE, L.LT.M, AND
  2214. C                        S(L), ..., S(M) ARE NOT NEGLIGIBLE (QR STEP).
  2215. C           KASE = 4     IF E(M-1) IS NEGLIGIBLE (CONVERGENCE).
  2216. C
  2217.          DO 390 LL = 1, M
  2218.             L = M - LL
  2219. C        ...EXIT
  2220.             IF (L .EQ. 0) GO TO 400
  2221.             TEST = DABS(S(L)) + DABS(S(L+1))
  2222.             ZTEST = TEST + DABS(E(L))
  2223.             IF (ZTEST .NE. TEST) GO TO 380
  2224.                E(L) = 0.0D0
  2225. C        ......EXIT
  2226.                GO TO 400
  2227.   380       CONTINUE
  2228.   390    CONTINUE
  2229.   400    CONTINUE
  2230.          IF (L .NE. M - 1) GO TO 410
  2231.             KASE = 4
  2232.          GO TO 480
  2233.   410    CONTINUE
  2234.             LP1 = L + 1
  2235.             MP1 = M + 1
  2236.             DO 430 LLS = LP1, MP1
  2237.                LS = M - LLS + LP1
  2238. C           ...EXIT
  2239.                IF (LS .EQ. L) GO TO 440
  2240.                TEST = 0.0D0
  2241.                IF (LS .NE. M) TEST = TEST + DABS(E(LS))
  2242.                IF (LS .NE. L + 1) TEST = TEST + DABS(E(LS-1))
  2243.                ZTEST = TEST + DABS(S(LS))
  2244.                IF (ZTEST .NE. TEST) GO TO 420
  2245.                   S(LS) = 0.0D0
  2246. C           ......EXIT
  2247.                   GO TO 440
  2248.   420          CONTINUE
  2249.   430       CONTINUE
  2250.   440       CONTINUE
  2251.             IF (LS .NE. L) GO TO 450
  2252.                KASE = 3
  2253.             GO TO 470
  2254.   450       CONTINUE
  2255.             IF (LS .NE. M) GO TO 460
  2256.                KASE = 1
  2257.             GO TO 470
  2258.   460       CONTINUE
  2259.                KASE = 2
  2260.                L = LS
  2261.   470       CONTINUE
  2262.   480    CONTINUE
  2263.          L = L + 1
  2264. C
  2265. C        PERFORM THE TASK INDICATED BY KASE.
  2266. C
  2267.          GO TO (490,520,540,570), KASE
  2268. C
  2269. C        DEFLATE NEGLIGIBLE S(M).
  2270. C
  2271.   490    CONTINUE
  2272.             MM1 = M - 1
  2273.             F = E(M-1)
  2274.             E(M-1) = 0.0D0
  2275.             DO 510 KK = L, MM1
  2276.                K = MM1 - KK + L
  2277.                T1 = S(K)
  2278.                CALL DROTG(T1,F,CS,SN)
  2279.                S(K) = T1
  2280.                IF (K .EQ. L) GO TO 500
  2281.                   F = -SN*E(K-1)
  2282.                   E(K-1) = CS*E(K-1)
  2283.   500          CONTINUE
  2284.                IF (WANTV) CALL DROT(P,V(1,K),1,V(1,M),1,CS,SN)
  2285.   510       CONTINUE
  2286.          GO TO 610
  2287. C
  2288. C        SPLIT AT NEGLIGIBLE S(L).
  2289. C
  2290.   520    CONTINUE
  2291.             F = E(L-1)
  2292.             E(L-1) = 0.0D0
  2293.             DO 530 K = L, M
  2294.                T1 = S(K)
  2295.                CALL DROTG(T1,F,CS,SN)
  2296.                S(K) = T1
  2297.                F = -SN*E(K)
  2298.                E(K) = CS*E(K)
  2299.                IF (WANTU) CALL DROT(N,U(1,K),1,U(1,L-1),1,CS,SN)
  2300.   530       CONTINUE
  2301.          GO TO 610
  2302. C
  2303. C        PERFORM ONE QR STEP.
  2304. C
  2305.   540    CONTINUE
  2306. C
  2307. C           CALCULATE THE SHIFT.
  2308. C
  2309.             SCALE = DMAX1(DABS(S(M)),DABS(S(M-1)),DABS(E(M-1)),
  2310.      *                    DABS(S(L)),DABS(E(L)))
  2311.             SM = S(M)/SCALE
  2312.             SMM1 = S(M-1)/SCALE
  2313.             EMM1 = E(M-1)/SCALE
  2314.             SL = S(L)/SCALE
  2315.             EL = E(L)/SCALE
  2316.             B = ((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0
  2317.             C = (SM*EMM1)**2
  2318.             SHIFT = 0.0D0
  2319.             IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 550
  2320.                SHIFT = DSQRT(B**2+C)
  2321.                IF (B .LT. 0.0D0) SHIFT = -SHIFT
  2322.                SHIFT = C/(B + SHIFT)
  2323.   550       CONTINUE
  2324.             F = (SL + SM)*(SL - SM) - SHIFT
  2325.             G = SL*EL
  2326. C
  2327. C           CHASE ZEROS.
  2328. C
  2329.             MM1 = M - 1
  2330.             DO 560 K = L, MM1
  2331.                CALL DROTG(F,G,CS,SN)
  2332.                IF (K .NE. L) E(K-1) = F
  2333.                F = CS*S(K) + SN*E(K)
  2334.                E(K) = CS*E(K) - SN*S(K)
  2335.                G = SN*S(K+1)
  2336.                S(K+1) = CS*S(K+1)
  2337.                IF (WANTV) CALL DROT(P,V(1,K),1,V(1,K+1),1,CS,SN)
  2338.                CALL DROTG(F,G,CS,SN)
  2339.                S(K) = F
  2340.                F = CS*E(K) + SN*S(K+1)
  2341.                S(K+1) = -SN*E(K) + CS*S(K+1)
  2342.                G = SN*E(K+1)
  2343.                E(K+1) = CS*E(K+1)
  2344.                IF (WANTU .AND. K .LT. N)
  2345.      *            CALL DROT(N,U(1,K),1,U(1,K+1),1,CS,SN)
  2346.   560       CONTINUE
  2347.             E(M-1) = F
  2348.             ITER = ITER + 1
  2349.          GO TO 610
  2350. C
  2351. C        CONVERGENCE.
  2352. C
  2353.   570    CONTINUE
  2354. C
  2355. C           MAKE THE SINGULAR VALUE  POSITIVE.
  2356. C
  2357.             IF (S(L) .GE. 0.0D0) GO TO 580
  2358.                S(L) = -S(L)
  2359.                IF (WANTV) CALL DSCAL(P,-1.0D0,V(1,L),1)
  2360.   580       CONTINUE
  2361. C
  2362. C           ORDER THE SINGULAR VALUE.
  2363. C
  2364.   590       IF (L .EQ. MM) GO TO 600
  2365. C           ...EXIT
  2366.                IF (S(L) .GE. S(L+1)) GO TO 600
  2367.                T = S(L)
  2368.                S(L) = S(L+1)
  2369.                S(L+1) = T
  2370.                IF (WANTV .AND. L .LT. P)
  2371.      *            CALL DSWAP(P,V(1,L),1,V(1,L+1),1)
  2372.                IF (WANTU .AND. L .LT. N)
  2373.      *            CALL DSWAP(N,U(1,L),1,U(1,L+1),1)
  2374.                L = L + 1
  2375.             GO TO 590
  2376.   600       CONTINUE
  2377.             ITER = 0
  2378.             M = M - 1
  2379.   610    CONTINUE
  2380.       GO TO 360
  2381.   620 CONTINUE
  2382.       RETURN
  2383.       END
  2384. C
  2385. C     ------------------------------------------------------------------
  2386. C
  2387. C     TITLE:  RS
  2388. C
  2389. C
  2390. C     THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF
  2391. C     SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK)
  2392. C     TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED)
  2393. C     OF A REAL SYMMETRIC MATRIX.
  2394. C
  2395. C     ON INPUT:
  2396. C
  2397. C        NM  MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL
  2398. C        ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2399. C        DIMENSION STATEMENT;
  2400. C
  2401. C        N  IS THE ORDER OF THE MATRIX  A;
  2402. C
  2403. C        A  CONTAINS THE REAL SYMMETRIC MATRIX;
  2404. C
  2405. C        MATZ  IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF
  2406. C        ONLY EIGENVALUES ARE DESIRED;  OTHERWISE IT IS SET TO
  2407. C        ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS.
  2408. C
  2409. C     ON OUTPUT:
  2410. C
  2411. C        W  CONTAINS THE EIGENVALUES IN ASCENDING ORDER;
  2412. C
  2413. C        Z  CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO;
  2414. C
  2415. C        IERR  IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN
  2416. C        ERROR COMPLETION CODE DESCRIBED IN SECTION 2B OF THE
  2417. C        DOCUMENTATION.  THE NORMAL COMPLETION CODE IS ZERO;
  2418. C
  2419. C        FV1  AND  FV2  ARE TEMPORARY STORAGE ARRAYS.
  2420. C
  2421. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  2422. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  2423. C
  2424. C     ------------------------------------------------------------------
  2425. C
  2426.       SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR)
  2427. C
  2428.       INTEGER N,NM,IERR,MATZ
  2429.       DOUBLE PRECISION A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N)
  2430. C
  2431.       IF (N .LE. NM) GO TO 10
  2432.       IERR = 10 * N
  2433.       GO TO 50
  2434. C
  2435.    10 IF (MATZ .NE. 0) GO TO 20
  2436. C
  2437. C::::: FIND EIGENVALUES ONLY ::::::::::
  2438.       CALL  TRED1(NM,N,A,W,FV1,FV2)
  2439.       CALL  TQLRAT(N,W,FV2,IERR)
  2440.       GO TO 50
  2441. C
  2442. C::::: FIND BOTH EIGENVALUES AND EIGENVECTORS
  2443. C:::::
  2444.    20 CALL  TRED2(NM,N,A,W,FV1,Z)
  2445.       CALL  TQL2(NM,N,W,FV1,Z,IERR)
  2446.    50 RETURN
  2447. C
  2448. C::::: LAST CARD OF RS
  2449. C::
  2450.       END
  2451. C
  2452. C     ------------------------------------------------------------------
  2453. C
  2454.       SUBROUTINE TRED1(NM,N,A,D,E,E2)
  2455. C
  2456.       INTEGER I,J,K,L,N,II,NM,JP1
  2457.       DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N)
  2458.       DOUBLE PRECISION F,G,H,SCALE
  2459.       DOUBLE PRECISION DSQRT,DABS,DSIGN
  2460. C
  2461. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1,
  2462. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  2463. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  2464. C
  2465. C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX
  2466. C     TO A SYMMETRIC TRIDIAGONAL MATRIX USING
  2467. C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  2468. C
  2469. C     ON INPUT:
  2470. C
  2471. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2472. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2473. C          DIMENSION STATEMENT;
  2474. C
  2475. C        N IS THE ORDER OF THE MATRIX;
  2476. C
  2477. C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
  2478. C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
  2479. C
  2480. C     ON OUTPUT:
  2481. C
  2482. C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS-
  2483. C          FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER
  2484. C          TRIANGLE.  THE FULL UPPER TRIANGLE OF A IS UNALTERED;
  2485. C
  2486. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX;
  2487. C
  2488. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  2489. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO;
  2490. C
  2491. C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
  2492. C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
  2493. C
  2494. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  2495. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  2496. C
  2497. C     ------------------------------------------------------------------
  2498. C
  2499.       DO 100 I = 1, N
  2500.   100 D(I) = A(I,I)
  2501. C
  2502. C::::: FOR I=N STEP -1 UNTIL 1 DO --
  2503. C::::::
  2504.       DO 300 II = 1, N
  2505.          I = N + 1 - II
  2506.          L = I - 1
  2507.          H = 0.0D0
  2508.          SCALE = 0.0D0
  2509.          IF (L .LT. 1) GO TO 130
  2510. C
  2511. C::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED)
  2512. C::::
  2513.          DO 120 K = 1, L
  2514.   120    SCALE = SCALE + DABS(A(I,K))
  2515. C
  2516.          IF (SCALE .NE. 0.0D0) GO TO 140
  2517.   130    E(I) = 0.0D0
  2518.          E2(I) = 0.0D0
  2519.          GO TO 290
  2520. C
  2521.   140    DO 150 K = 1, L
  2522.             A(I,K) = A(I,K) / SCALE
  2523.             H = H + A(I,K) * A(I,K)
  2524.   150    CONTINUE
  2525. C
  2526.          E2(I) = SCALE * SCALE * H
  2527.          F = A(I,L)
  2528.          G = -DSIGN(DSQRT(H),F)
  2529.          E(I) = SCALE * G
  2530.          H = H - F * G
  2531.          A(I,L) = F - G
  2532.          IF (L .EQ. 1) GO TO 270
  2533.          F = 0.0D0
  2534. C
  2535.          DO 240 J = 1, L
  2536.             G = 0.0D0
  2537. C
  2538. C::::: FORM ELEMENT OF A*U
  2539. C::::::
  2540.             DO 180 K = 1, J
  2541.   180       G = G + A(J,K) * A(I,K)
  2542. C
  2543.             JP1 = J + 1
  2544.             IF (L .LT. JP1) GO TO 220
  2545. C
  2546.             DO 200 K = JP1, L
  2547.   200       G = G + A(K,J) * A(I,K)
  2548. C
  2549. C::::: FORM ELEMENT OF P
  2550. C::::
  2551.   220       E(J) = G / H
  2552.             F = F + E(J) * A(I,J)
  2553.   240    CONTINUE
  2554. C
  2555.          H = F / (H + H)
  2556. C
  2557. C::::: FORM REDUCED A
  2558. C:
  2559.          DO 260 J = 1, L
  2560.             F = A(I,J)
  2561.             G = E(J) - H * F
  2562.             E(J) = G
  2563. C
  2564.             DO 260 K = 1, J
  2565.                A(J,K) = A(J,K) - F * E(K) - G * A(I,K)
  2566.   260    CONTINUE
  2567. C
  2568.   270    DO 280 K = 1, L
  2569.   280    A(I,K) = SCALE * A(I,K)
  2570. C
  2571.   290    H = D(I)
  2572.          D(I) = A(I,I)
  2573.          A(I,I) = H
  2574.   300 CONTINUE
  2575. C
  2576.       RETURN
  2577. C
  2578. C::::: LAST CARD OF TRED1
  2579. C:::::
  2580.       END
  2581. C
  2582. C     ------------------------------------------------------------------
  2583. C
  2584.       SUBROUTINE TQLRAT(N,D,E2,IERR)
  2585. C
  2586.       INTEGER I,J,L,M,N,II,L1,MML,IERR
  2587.       DOUBLE PRECISION D(N),E2(N)
  2588.       DOUBLE PRECISION B,C,F,G,H,P,R,S,T,MACHEP
  2589.       DOUBLE PRECISION DSQRT,DABS,DSIGN
  2590. C
  2591. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
  2592. C     ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
  2593. C
  2594. C     THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
  2595. C     TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD.
  2596. C
  2597. C     ON INPUT:
  2598. C
  2599. C        N IS THE ORDER OF THE MATRIX;
  2600. C
  2601. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX;
  2602. C
  2603. C        E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE
  2604. C          INPUT MATRIX IN ITS LAST N-1 POSITIONS.  E2(1) IS ARBITRARY.
  2605. C
  2606. C      ON OUTPUT:
  2607. C
  2608. C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
  2609. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
  2610. C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
  2611. C          THE SMALLEST EIGENVALUES;
  2612. C
  2613. C        E2 HAS BEEN DESTROYED;
  2614. C
  2615. C        IERR IS SET TO
  2616. C          ZERO       FOR NORMAL RETURN,
  2617. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
  2618. C                     DETERMINED AFTER 30 ITERATIONS.
  2619. C
  2620. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  2621. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  2622. C
  2623. C     ------------------------------------------------------------------
  2624. C
  2625. C
  2626. C::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
  2627. C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
  2628. C                AN APPROXIMATION TO MACHEP IS COMPUTED BELOW.
  2629. C
  2630. C::::::
  2631.       MACHEP = 1.D0 / 2.D0
  2632.       DO 2 I=1,500
  2633.          T = 1.D0 + MACHEP
  2634.          IF (T .EQ. 1.D0) GO TO 3
  2635.          MACHEP = MACHEP / 2.D0
  2636.     2 CONTINUE
  2637. C
  2638.     3 IERR = 0
  2639.       IF (N .EQ. 1) GO TO 1001
  2640. C
  2641.       DO 100 I = 2, N
  2642.   100 E2(I-1) = E2(I)
  2643. C
  2644.       F = 0.0D0
  2645.       B = 0.0D0
  2646.       E2(N) = 0.0D0
  2647. C
  2648.       DO 290 L = 1, N
  2649.          J = 0
  2650.          H = MACHEP * (DABS(D(L)) + DSQRT(E2(L)))
  2651.          IF (B .GT. H) GO TO 105
  2652.          B = H
  2653.          C = B * B
  2654. C
  2655. C::::: LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT
  2656. C
  2657.   105    DO 110 M = L, N
  2658.             IF (E2(M) .LE. C) GO TO 120
  2659. C
  2660. C::::: E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
  2661. C                THROUGH THE BOTTOM OF THE LOOP
  2662. C:::::::
  2663.   110    CONTINUE
  2664. C
  2665.   120    IF (M .EQ. L) GO TO 210
  2666.   130    IF (J .EQ. 30) GO TO 1000
  2667.          J = J + 1
  2668. C
  2669. C::::: FORM SHIFT
  2670. C:::::::
  2671.          L1 = L + 1
  2672.          S = DSQRT(E2(L))
  2673.          G = D(L)
  2674.          P = (D(L1) - G) / (2.0D0 * S)
  2675.          R = DSQRT(P*P+1.0D0)
  2676.          D(L) = S / (P + DSIGN(R,P))
  2677.          H = G - D(L)
  2678. C
  2679.          DO 140 I = L1, N
  2680.   140    D(I) = D(I) - H
  2681. C
  2682.          F = F + H
  2683. C
  2684. C::::: RATIONAL QL TRANSFORMATION
  2685. C:::
  2686.          G = D(M)
  2687.          IF (G .EQ. 0.0D0) G = B
  2688.          H = G
  2689.          S = 0.0D0
  2690.          MML = M - L
  2691. C
  2692. C::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
  2693.          DO 200 II = 1, MML
  2694.             I = M - II
  2695.             P = G * H
  2696.             R = P + E2(I)
  2697.             E2(I+1) = S * R
  2698.             S = E2(I) / R
  2699.             D(I+1) = H + S * (H + D(I))
  2700.             G = D(I) - E2(I) / G
  2701.             IF (G .EQ. 0.0D0) G = B
  2702.             H = G * P / R
  2703.   200    CONTINUE
  2704. C
  2705.          E2(L) = S * G
  2706.          D(L) = H
  2707. C
  2708. C::::: GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST
  2709. C
  2710.          IF (H .EQ. 0.0D0) GO TO 210
  2711.          IF (DABS(E2(L)) .LE. DABS(C/H)) GO TO 210
  2712.          E2(L) = H * E2(L)
  2713.          IF (E2(L) .NE. 0.0D0) GO TO 130
  2714.   210    P = D(L) + F
  2715. C
  2716. C::::: ORDER EIGENVALUES
  2717. C::::
  2718.          IF (L .EQ. 1) GO TO 250
  2719. C
  2720. C::::: FOR I=L STEP -1 UNTIL 2 DO --
  2721. C::::::
  2722.          DO 230 II = 2, L
  2723.             I = L + 2 - II
  2724.             IF (P .GE. D(I-1)) GO TO 270
  2725.             D(I) = D(I-1)
  2726.   230    CONTINUE
  2727. C
  2728.   250    I = 1
  2729.   270    D(I) = P
  2730.   290 CONTINUE
  2731. C
  2732.       GO TO 1001
  2733. C
  2734. C::::: SET ERROR -- NO CONVERGENCE TO AN
  2735. C                EIGENVALUE AFTER 30 ITERATIONS
  2736. C:::::::
  2737.  1000 IERR = L
  2738.  1001 RETURN
  2739. C
  2740. C::::: LAST CARD OF TQLRAT
  2741. C::::::
  2742.       END
  2743. C
  2744. C     ------------------------------------------------------------------
  2745. C
  2746.       SUBROUTINE TRED2(NM,N,A,D,E,Z)
  2747. C
  2748.       INTEGER I,J,K,L,N,II,NM,JP1
  2749.       DOUBLE PRECISION A(NM,N),D(N),E(N),Z(NM,N)
  2750.       DOUBLE PRECISION F,G,H,HH,SCALE
  2751.       DOUBLE PRECISION DSQRT,DABS,DSIGN
  2752. C
  2753. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2,
  2754. C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
  2755. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
  2756. C
  2757. C     THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A
  2758. C     SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING
  2759. C     ORTHOGONAL SIMILARITY TRANSFORMATIONS.
  2760. C
  2761. C     ON INPUT:
  2762. C
  2763. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2764. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2765. C          DIMENSION STATEMENT;
  2766. C
  2767. C        N IS THE ORDER OF THE MATRIX;
  2768. C
  2769. C        A CONTAINS THE REAL SYMMETRIC INPUT MATRIX.  ONLY THE
  2770. C          LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
  2771. C
  2772. C     ON OUTPUT:
  2773. C
  2774. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX;
  2775. C
  2776. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
  2777. C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO;
  2778. C
  2779. C        Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX
  2780. C          PRODUCED IN THE REDUCTION;
  2781. C
  2782. C        A AND Z MAY COINCIDE.  IF DISTINCT, A IS UNALTERED.
  2783. C
  2784. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  2785. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  2786. C
  2787. C     ------------------------------------------------------------------
  2788. C
  2789.       DO 100 I = 1, N
  2790. C
  2791.          DO 100 J = 1, I
  2792.             Z(I,J) = A(I,J)
  2793.   100 CONTINUE
  2794. C
  2795.       IF (N .EQ. 1) GO TO 320
  2796. C
  2797. C::::: FOR I=N STEP -1 UNTIL 2 DO --
  2798. C::::::
  2799.       DO 300 II = 2, N
  2800.          I = N + 2 - II
  2801.          L = I - 1
  2802.          H = 0.0D0
  2803.          SCALE = 0.0D0
  2804.          IF (L .LT. 2) GO TO 130
  2805. C
  2806. C::::: SCALE ROW (ALGOL TOL THEN NOT NEEDED)
  2807. C::::
  2808.          DO 120 K = 1, L
  2809.   120    SCALE = SCALE + DABS(Z(I,K))
  2810. C
  2811.          IF (SCALE .NE. 0.0D0) GO TO 140
  2812.   130    E(I) = Z(I,L)
  2813.          GO TO 290
  2814. C
  2815.   140    DO 150 K = 1, L
  2816.             Z(I,K) = Z(I,K) / SCALE
  2817.             H = H + Z(I,K) * Z(I,K)
  2818.   150    CONTINUE
  2819. C
  2820.          F = Z(I,L)
  2821.          G = -DSIGN(DSQRT(H),F)
  2822.          E(I) = SCALE * G
  2823.          H = H - F * G
  2824.          Z(I,L) = F - G
  2825.          F = 0.0D0
  2826. C
  2827.          DO 240 J = 1, L
  2828.             Z(J,I) = Z(I,J) / H
  2829.             G = 0.0D0
  2830. C
  2831. C::::: FORM ELEMENT OF A*U
  2832. C::::::
  2833.             DO 180 K = 1, J
  2834.   180       G = G + Z(J,K) * Z(I,K)
  2835. C
  2836.             JP1 = J + 1
  2837.             IF (L .LT. JP1) GO TO 220
  2838. C
  2839.             DO 200 K = JP1, L
  2840.   200       G = G + Z(K,J) * Z(I,K)
  2841. C
  2842. C::::: FORM ELEMENT OF P
  2843. C::::
  2844.   220       E(J) = G / H
  2845.             F = F + E(J) * Z(I,J)
  2846.   240    CONTINUE
  2847. C
  2848.          HH = F / (H + H)
  2849. C
  2850. C::::: FORM REDUCED A
  2851. C:
  2852.          DO 260 J = 1, L
  2853.             F = Z(I,J)
  2854.             G = E(J) - HH * F
  2855.             E(J) = G
  2856. C
  2857.             DO 260 K = 1, J
  2858.                Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K)
  2859.   260    CONTINUE
  2860. C
  2861.   290    D(I) = H
  2862.   300 CONTINUE
  2863. C
  2864.   320 D(1) = 0.0D0
  2865.       E(1) = 0.0D0
  2866. C
  2867. C::::: ACCUMULATION OF TRANSFORMATION MATRICES
  2868. C::::::
  2869.       DO 500 I = 1, N
  2870.          L = I - 1
  2871.          IF (D(I) .EQ. 0.0D0) GO TO 380
  2872. C
  2873.          DO 360 J = 1, L
  2874.             G = 0.0D0
  2875. C
  2876.             DO 340 K = 1, L
  2877.   340       G = G + Z(I,K) * Z(K,J)
  2878. C
  2879.             DO 360 K = 1, L
  2880.                Z(K,J) = Z(K,J) - G * Z(K,I)
  2881.   360    CONTINUE
  2882. C
  2883.   380    D(I) = Z(I,I)
  2884.          Z(I,I) = 1.0D0
  2885.          IF (L .LT. 1) GO TO 500
  2886. C
  2887.          DO 400 J = 1, L
  2888.             Z(I,J) = 0.0D0
  2889.             Z(J,I) = 0.0D0
  2890.   400    CONTINUE
  2891. C
  2892.   500 CONTINUE
  2893. C
  2894.       RETURN
  2895. C
  2896. C::::: LAST CARD OF TRED2
  2897. C:::::
  2898.       END
  2899. C
  2900. C     ------------------------------------------------------------------
  2901. C
  2902.       SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
  2903. C
  2904.       INTEGER I,J,K,L,M,N,II,L1,NM,MML,IERR
  2905.       DOUBLE PRECISION D(N),E(N),Z(NM,N)
  2906.       DOUBLE PRECISION B,C,F,G,H,P,R,S,T,MACHEP
  2907.       DOUBLE PRECISION DSQRT,DABS,DSIGN
  2908. C
  2909. C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
  2910. C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
  2911. C     WILKINSON.
  2912. C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
  2913. C
  2914. C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
  2915. C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
  2916. C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
  2917. C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
  2918. C     FULL MATRIX TO TRIDIAGONAL FORM.
  2919. C
  2920. C     ON INPUT:
  2921. C
  2922. C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  2923. C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM
  2924. C          DIMENSION STATEMENT;
  2925. C
  2926. C        N IS THE ORDER OF THE MATRIX;
  2927. C
  2928. C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX;
  2929. C
  2930. C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  2931. C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY;
  2932. C
  2933. C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
  2934. C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
  2935. C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
  2936. C          THE IDENTITY MATRIX.
  2937. C
  2938. C      ON OUTPUT:
  2939. C
  2940. C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
  2941. C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
  2942. C          UNORDERED FOR INDICES 1,2,...,IERR-1;
  2943. C
  2944. C        E HAS BEEN DESTROYED;
  2945. C
  2946. C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
  2947. C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
  2948. C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
  2949. C          EIGENVALUES;
  2950. C
  2951. C        IERR IS SET TO
  2952. C          ZERO       FOR NORMAL RETURN,
  2953. C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
  2954. C                     DETERMINED AFTER 30 ITERATIONS.
  2955. C
  2956. C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
  2957. C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
  2958. C
  2959. C     ------------------------------------------------------------------
  2960. C
  2961. C
  2962. C::::: MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
  2963. C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
  2964. C                AN APPROXIMATION TO MACHEP IS COMPUTED BELOW.
  2965. C
  2966. C::::::
  2967.       MACHEP = 1.D0 / 2.D0
  2968.       DO 2 I=1,500
  2969.          T = 1.D0 + MACHEP
  2970.          IF (T .EQ. 1.D0) GO TO 3
  2971.          MACHEP = MACHEP / 2.D0
  2972.     2 CONTINUE
  2973. C
  2974.     3 IERR = 0
  2975.       IF (N .EQ. 1) GO TO 1001
  2976. C
  2977.       DO 100 I = 2, N
  2978.   100 E(I-1) = E(I)
  2979. C
  2980.       F = 0.0D0
  2981.       B = 0.0D0
  2982.       E(N) = 0.0D0
  2983. C
  2984.       DO 240 L = 1, N
  2985.          J = 0
  2986.          H = MACHEP * (DABS(D(L)) + DABS(E(L)))
  2987.          IF (B .LT. H) B = H
  2988. C
  2989. C::::: LOOK FOR SMALL SUB-DIAGONAL ELEMENT
  2990. C::
  2991.          DO 110 M = L, N
  2992.             IF (DABS(E(M)) .LE. B) GO TO 120
  2993. C
  2994. C::::: E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
  2995. C                THROUGH THE BOTTOM OF THE LOOP
  2996. C:::::::
  2997.   110    CONTINUE
  2998. C
  2999.   120    IF (M .EQ. L) GO TO 220
  3000.   130    IF (J .EQ. 30) GO TO 1000
  3001.          J = J + 1
  3002. C
  3003. C::::: FORM SHIFT
  3004. C:::::::
  3005.          L1 = L + 1
  3006.          G = D(L)
  3007.          P = (D(L1) - G) / (2.0D0 * E(L))
  3008.          R = DSQRT(P*P+1.0D0)
  3009.          D(L) = E(L) / (P + DSIGN(R,P))
  3010.          H = G - D(L)
  3011. C
  3012.          DO 140 I = L1, N
  3013.   140    D(I) = D(I) - H
  3014. C
  3015.          F = F + H
  3016. C
  3017. C::::: QL TRANSFORMATION
  3018. C::::
  3019.          P = D(M)
  3020.          C = 1.0D0
  3021.          S = 0.0D0
  3022.          MML = M - L
  3023. C
  3024. C::::: FOR I=M-1 STEP -1 UNTIL L DO -- ::::::::::
  3025.          DO 200 II = 1, MML
  3026.             I = M - II
  3027.             G = C * E(I)
  3028.             H = C * P
  3029.             IF (DABS(P) .LT. DABS(E(I))) GO TO 150
  3030.             C = E(I) / P
  3031.             R = DSQRT(C*C+1.0D0)
  3032.             E(I+1) = S * P * R
  3033.             S = C / R
  3034.             C = 1.0D0 / R
  3035.             GO TO 160
  3036.   150       C = P / E(I)
  3037.             R = DSQRT(C*C+1.0D0)
  3038.             E(I+1) = S * E(I) * R
  3039.             S = 1.0D0 / R
  3040.             C = C * S
  3041.   160       P = C * D(I) - S * G
  3042.             D(I+1) = H + S * (C * G + S * D(I))
  3043. C
  3044. C::::: FORM VECTOR ::::::::::
  3045.             DO 180 K = 1, N
  3046.                H = Z(K,I+1)
  3047.                Z(K,I+1) = S * Z(K,I) + C * H
  3048.                Z(K,I) = C * Z(K,I) - S * H
  3049.   180       CONTINUE
  3050. C
  3051.   200    CONTINUE
  3052. C
  3053.          E(L) = S * P
  3054.          D(L) = C * P
  3055.          IF (DABS(E(L)) .GT. B) GO TO 130
  3056.   220    D(L) = D(L) + F
  3057.   240 CONTINUE
  3058. C
  3059. C::::: ORDER EIGENVALUES AND EIGENVECTORS
  3060. C:
  3061.       DO 300 II = 2, N
  3062.          I = II - 1
  3063.          K = I
  3064.          P = D(I)
  3065. C
  3066.          DO 260 J = II, N
  3067.             IF (D(J) .GE. P) GO TO 260
  3068.             K = J
  3069.             P = D(J)
  3070.   260    CONTINUE
  3071. C
  3072.          IF (K .EQ. I) GO TO 300
  3073.          D(K) = D(I)
  3074.          D(I) = P
  3075. C
  3076.          DO 280 J = 1, N
  3077.             P = Z(J,I)
  3078.             Z(J,I) = Z(J,K)
  3079.             Z(J,K) = P
  3080.   280    CONTINUE
  3081. C
  3082.   300 CONTINUE
  3083. C
  3084.       GO TO 1001
  3085. C
  3086. C::::: SET ERROR -- NO CONVERGENCE TO AN
  3087. C                EIGENVALUE AFTER 30 ITERATIONS
  3088. C:::::::
  3089.  1000 IERR = L
  3090.  1001 RETURN
  3091. C
  3092. C::::: LAST CARD OF TQL2
  3093. C::::
  3094.       END
  3095.       SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
  3096. C
  3097. C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
  3098. C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
  3099. C     JACK DONGARRA, LINPACK, 3/11/78.
  3100. C
  3101.       DOUBLE PRECISION DX(1),DY(1),DA
  3102.       INTEGER I,INCX,INCY,IXIY,M,MP1,N
  3103. C
  3104.       IF(N.LE.0)RETURN
  3105.       IF (DA .EQ. 0.0D0) RETURN
  3106.       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
  3107. C
  3108. C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
  3109. C          NOT EQUAL TO 1
  3110. C
  3111.       IX = 1
  3112.       IY = 1
  3113.       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
  3114.       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
  3115.       DO 10 I = 1,N
  3116.         DY(IY) = DY(IY) + DA*DX(IX)
  3117.         IX = IX + INCX
  3118.         IY = IY + INCY
  3119.    10 CONTINUE
  3120.       RETURN
  3121. C
  3122. C        CODE FOR BOTH INCREMENTS EQUAL TO 1
  3123. C
  3124. C
  3125. C        CLEAN-UP LOOP
  3126. C
  3127.    20 M = MOD(N,4)
  3128.       IF( M .EQ. 0 ) GO TO 40
  3129.       DO 30 I = 1,M
  3130.         DY(I) = DY(I) + DA*DX(I)
  3131.    30 CONTINUE
  3132.       IF( N .LT. 4 ) RETURN
  3133.    40 MP1 = M + 1
  3134.       DO 50 I = MP1,N,4
  3135.         DY(I) = DY(I) + DA*DX(I)
  3136.         DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
  3137.         DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
  3138.         DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
  3139.    50 CONTINUE
  3140.       RETURN
  3141.       END
  3142.       SUBROUTINE  DCOPY(N,DX,INCX,DY,INCY)
  3143. C
  3144. C     COPIES A VECTOR, X, TO A VECTOR, Y.
  3145. C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
  3146. C     JACK DONGARRA, LINPACK, 3/11/78.
  3147. C
  3148.       DOUBLE PRECISION DX(1),DY(1)
  3149.       INTEGER I,INCX,INCY,IX,IY,M,MP1,N
  3150. C
  3151.       IF(N.LE.0)RETURN
  3152.       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
  3153. C
  3154. C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
  3155. C          NOT EQUAL TO 1
  3156. C
  3157.       IX = 1
  3158.       IY = 1
  3159.       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
  3160.       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
  3161.       DO 10 I = 1,N
  3162.         DY(IY) = DX(IX)
  3163.         IX = IX + INCX
  3164.         IY = IY + INCY
  3165.    10 CONTINUE
  3166.       RETURN
  3167. C
  3168. C        CODE FOR BOTH INCREMENTS EQUAL TO 1
  3169. C
  3170. C
  3171. C        CLEAN-UP LOOP
  3172. C
  3173.    20 M = MOD(N,7)
  3174.       IF( M .EQ. 0 ) GO TO 40
  3175.       DO 30 I = 1,M
  3176.         DY(I) = DX(I)
  3177.    30 CONTINUE
  3178.       IF( N .LT. 7 ) RETURN
  3179.    40 MP1 = M + 1
  3180.       DO 50 I = MP1,N,7
  3181.         DY(I) = DX(I)
  3182.         DY(I + 1) = DX(I + 1)
  3183.         DY(I + 2) = DX(I + 2)
  3184.         DY(I + 3) = DX(I + 3)
  3185.         DY(I + 4) = DX(I + 4)
  3186.         DY(I + 5) = DX(I + 5)
  3187.         DY(I + 6) = DX(I + 6)
  3188.    50 CONTINUE
  3189.       RETURN
  3190.       END
  3191.       DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY)
  3192. C
  3193. C     FORMS THE DOT PRODUCT OF TWO VECTORS.
  3194. C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
  3195. C     JACK DONGARRA, LINPACK, 3/11/78.
  3196. C
  3197.       DOUBLE PRECISION DX(1),DY(1),DTEMP
  3198.       INTEGER I,INCX,INCY,IX,IY,M,MP1,N
  3199. C
  3200.       DDOT = 0.0D0
  3201.       DTEMP = 0.0D0
  3202.       IF(N.LE.0)RETURN
  3203.       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
  3204. C
  3205. C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
  3206. C          NOT EQUAL TO 1
  3207. C
  3208.       IX = 1
  3209.       IY = 1
  3210.       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
  3211.       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
  3212.       DO 10 I = 1,N
  3213.         DTEMP = DTEMP + DX(IX)*DY(IY)
  3214.         IX = IX + INCX
  3215.         IY = IY + INCY
  3216.    10 CONTINUE
  3217.       DDOT = DTEMP
  3218.       RETURN
  3219. C
  3220. C        CODE FOR BOTH INCREMENTS EQUAL TO 1
  3221. C
  3222. C
  3223. C        CLEAN-UP LOOP
  3224. C
  3225.    20 M = MOD(N,5)
  3226.       IF( M .EQ. 0 ) GO TO 40
  3227.       DO 30 I = 1,M
  3228.         DTEMP = DTEMP + DX(I)*DY(I)
  3229.    30 CONTINUE
  3230.       IF( N .LT. 5 ) GO TO 60
  3231.    40 MP1 = M + 1
  3232.       DO 50 I = MP1,N,5
  3233.         DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) +
  3234.      *   DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
  3235.    50 CONTINUE
  3236.    60 DDOT = DTEMP
  3237.       RETURN
  3238.       END
  3239.       DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX)
  3240.       INTEGER          NEXT
  3241.       DOUBLE PRECISION   DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
  3242.       DATA   ZERO, ONE /0.0D0, 1.0D0/
  3243. C
  3244. C     EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
  3245. C     INCREMENT INCX .
  3246. C     IF    N .LE. 0 RETURN WITH RESULT = 0.
  3247. C     IF N .GE. 1 THEN INCX MUST BE .GE. 1
  3248. C
  3249. C           C.L.LAWSON, 1978 JAN 08
  3250. C
  3251. C     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
  3252. C     HOPEFULLY APPLICABLE TO ALL MACHINES.
  3253. C         CUTLO = MAXIMUM OF  DSQRT(U/EPS)  OVER ALL KNOWN MACHINES.
  3254. C         CUTHI = MINIMUM OF  DSQRT(V)      OVER ALL KNOWN MACHINES.
  3255. C     WHERE
  3256. C         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
  3257. C         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
  3258. C         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
  3259. C
  3260. C     BRIEF OUTLINE OF ALGORITHM..
  3261. C
  3262. C     PHASE 1    SCANS ZERO COMPONENTS.
  3263. C     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
  3264. C     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
  3265. C     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
  3266. C     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
  3267. C
  3268. C     VALUES FOR CUTLO AND CUTHI..
  3269. C     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
  3270. C     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
  3271. C     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
  3272. C                   UNIVAC AND DEC AT 2**(-103)
  3273. C                   THUS CUTLO = 2**(-51) = 4.44089E-16
  3274. C     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
  3275. C                   THUS CUTHI = 2**(63.5) = 1.30438E19
  3276. C     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
  3277. C                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
  3278. C     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D19
  3279. C     DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
  3280. C     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
  3281.       DATA CUTLO, CUTHI / 8.232D-11,  1.304D19 /
  3282. C
  3283.       IF(N .GT. 0) GO TO 10
  3284.          DNRM2  = ZERO
  3285.          GO TO 300
  3286. C
  3287.    10 ASSIGN 30 TO NEXT
  3288.       SUM = ZERO
  3289.       NN = N * INCX
  3290. C                                                 BEGIN MAIN LOOP
  3291.       I = 1
  3292.    20    GO TO NEXT,(30, 50, 70, 110)
  3293.    30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
  3294.       ASSIGN 50 TO NEXT
  3295.       XMAX = ZERO
  3296. C
  3297. C                        PHASE 1.  SUM IS ZERO
  3298. C
  3299.    50 IF( DX(I) .EQ. ZERO) GO TO 200
  3300.       IF( DABS(DX(I)) .GT. CUTLO) GO TO 85
  3301. C
  3302. C                                PREPARE FOR PHASE 2.
  3303.       ASSIGN 70 TO NEXT
  3304.       GO TO 105
  3305. C
  3306. C                                PREPARE FOR PHASE 4.
  3307. C
  3308.   100 I = J
  3309.       ASSIGN 110 TO NEXT
  3310.       SUM = (SUM / DX(I)) / DX(I)
  3311.   105 XMAX = DABS(DX(I))
  3312.       GO TO 115
  3313. C
  3314. C                   PHASE 2.  SUM IS SMALL.
  3315. C                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
  3316. C
  3317.    70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75
  3318. C
  3319. C                     COMMON CODE FOR PHASES 2 AND 4.
  3320. C                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
  3321. C
  3322.   110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115
  3323.          SUM = ONE + SUM * (XMAX / DX(I))**2
  3324.          XMAX = DABS(DX(I))
  3325.          GO TO 200
  3326. C
  3327.   115 SUM = SUM + (DX(I)/XMAX)**2
  3328.       GO TO 200
  3329. C
  3330. C
  3331. C                  PREPARE FOR PHASE 3.
  3332. C
  3333.    75 SUM = (SUM * XMAX) * XMAX
  3334. C
  3335. C
  3336. C     FOR REAL OR D.P. SET HITEST = CUTHI/N
  3337. C     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
  3338. C
  3339.    85 HITEST = CUTHI/FLOAT( N )
  3340. C
  3341. C                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
  3342. C
  3343.       DO 95 J =I,NN,INCX
  3344.       IF(DABS(DX(J)) .GE. HITEST) GO TO 100
  3345.    95    SUM = SUM + DX(J)**2
  3346.       DNRM2 = DSQRT( SUM )
  3347.       GO TO 300
  3348. C
  3349.   200 CONTINUE
  3350.       I = I + INCX
  3351.       IF ( I .LE. NN ) GO TO 20
  3352. C
  3353. C              END OF MAIN LOOP.
  3354. C
  3355. C              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
  3356. C
  3357.       DNRM2 = XMAX * DSQRT(SUM)
  3358.   300 CONTINUE
  3359.       RETURN
  3360.       END
  3361.       SUBROUTINE  DROT (N,DX,INCX,DY,INCY,C,S)                          
  3362. 00000010
  3363. C
  3364. C     APPLIES A PLANE ROTATION.
  3365. C     JACK DONGARRA, LINPACK, 3/11/78.
  3366. C
  3367.       DOUBLE PRECISION DX(1),DY(1),DTEMP,C,S
  3368.       INTEGER I,INCX,INCY,IX,IY,N
  3369. C
  3370.       IF(N.LE.0)RETURN
  3371.       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
  3372. C
  3373. C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
  3374. C         TO 1
  3375. C
  3376.       IX = 1
  3377.       IY = 1
  3378.       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
  3379.       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
  3380.       DO 10 I = 1,N
  3381.         DTEMP = C*DX(IX) + S*DY(IY)
  3382.         DY(IY) = C*DY(IY) - S*DX(IX)
  3383.         DX(IX) = DTEMP
  3384.         IX = IX + INCX
  3385.         IY = IY + INCY
  3386.    10 CONTINUE
  3387.       RETURN
  3388. C
  3389. C       CODE FOR BOTH INCREMENTS EQUAL TO 1
  3390. C
  3391.    20 DO 30 I = 1,N
  3392.         DTEMP = C*DX(I) + S*DY(I)
  3393.         DY(I) = C*DY(I) - S*DX(I)
  3394.         DX(I) = DTEMP
  3395.    30 CONTINUE
  3396.       RETURN
  3397.       END
  3398.       SUBROUTINE DROTG(DA,DB,C,S)                                       
  3399. 00000010
  3400. C
  3401. C     CONSTRUCT GIVENS PLANE ROTATION.
  3402. C     JACK DONGARRA, LINPACK, 3/11/78.
  3403. C
  3404.       DOUBLE PRECISION DA,DB,C,S,ROE,SCALE,R,Z
  3405. C
  3406.       ROE = DB
  3407.       IF( DABS(DA) .GT. DABS(DB) ) ROE = DA
  3408.       SCALE = DABS(DA) + DABS(DB)
  3409.       IF( SCALE .NE. 0.0D0 ) GO TO 10
  3410.          C = 1.0D0
  3411.          S = 0.0D0
  3412.          R = 0.0D0
  3413.          GO TO 20
  3414.    10 R = SCALE*DSQRT((DA/SCALE)**2 + (DB/SCALE)**2)
  3415.       R = DSIGN(1.0D0,ROE)*R
  3416.       C = DA/R
  3417.       S = DB/R
  3418.    20 Z = 1.0D0
  3419.       IF( DABS(DA) .GT. DABS(DB) ) Z = S
  3420.       IF( DABS(DB) .GE. DABS(DA) .AND. C .NE. 0.0D0 ) Z = 1.0D0/C
  3421.       DA = R
  3422.       DB = Z
  3423.       RETURN
  3424.       END
  3425.       SUBROUTINE  DSCAL(N,DA,DX,INCX)
  3426. C
  3427. C     SCALES A VECTOR BY A CONSTANT.
  3428. C     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
  3429. C     JACK DONGARRA, LINPACK, 3/11/78.
  3430. C
  3431.       DOUBLE PRECISION DA,DX(1)
  3432.       INTEGER I,INCX,M,MP1,N,NINCX
  3433. C
  3434.       IF(N.LE.0)RETURN
  3435.       IF(INCX.EQ.1)GO TO 20
  3436. C
  3437. C        CODE FOR INCREMENT NOT EQUAL TO 1
  3438. C
  3439.       NINCX = N*INCX
  3440.       DO 10 I = 1,NINCX,INCX
  3441.         DX(I) = DA*DX(I)
  3442.    10 CONTINUE
  3443.       RETURN
  3444. C
  3445. C        CODE FOR INCREMENT EQUAL TO 1
  3446. C
  3447. C
  3448. C        CLEAN-UP LOOP
  3449. C
  3450.    20 M = MOD(N,5)
  3451.       IF( M .EQ. 0 ) GO TO 40
  3452.       DO 30 I = 1,M
  3453.         DX(I) = DA*DX(I)
  3454.    30 CONTINUE
  3455.       IF( N .LT. 5 ) RETURN
  3456.    40 MP1 = M + 1
  3457.       DO 50 I = MP1,N,5
  3458.         DX(I) = DA*DX(I)
  3459.         DX(I + 1) = DA*DX(I + 1)
  3460.         DX(I + 2) = DA*DX(I + 2)
  3461.         DX(I + 3) = DA*DX(I + 3)
  3462.         DX(I + 4) = DA*DX(I + 4)
  3463.    50 CONTINUE
  3464.       RETURN
  3465.       END
  3466.       SUBROUTINE  DSWAP (N,DX,INCX,DY,INCY)
  3467. C
  3468. C     INTERCHANGES TWO VECTORS.
  3469. C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
  3470. C     JACK DONGARRA, LINPACK, 3/11/78.
  3471. C
  3472.       DOUBLE PRECISION DX(1),DY(1),DTEMP
  3473.       INTEGER I,INCX,INCY,IX,IY,M,MP1,N
  3474. C
  3475.       IF(N.LE.0)RETURN
  3476.       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
  3477. C
  3478. C       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
  3479. C         TO 1
  3480. C
  3481.       IX = 1
  3482.       IY = 1
  3483.       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
  3484.       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
  3485.       DO 10 I = 1,N
  3486.         DTEMP = DX(IX)
  3487.         DX(IX) = DY(IY)
  3488.         DY(IY) = DTEMP
  3489.         IX = IX + INCX
  3490.         IY = IY + INCY
  3491.    10 CONTINUE
  3492.       RETURN
  3493. C
  3494. C       CODE FOR BOTH INCREMENTS EQUAL TO 1
  3495. C
  3496. C
  3497. C       CLEAN-UP LOOP
  3498. C
  3499.    20 M = MOD(N,3)
  3500.       IF( M .EQ. 0 ) GO TO 40
  3501.       DO 30 I = 1,M
  3502.         DTEMP = DX(I)
  3503.         DX(I) = DY(I)
  3504.         DY(I) = DTEMP
  3505.    30 CONTINUE
  3506.       IF( N .LT. 3 ) RETURN
  3507.    40 MP1 = M + 1
  3508.       DO 50 I = MP1,N,3
  3509.         DTEMP = DX(I)
  3510.         DX(I) = DY(I)
  3511.         DY(I) = DTEMP
  3512.         DTEMP = DX(I + 1)
  3513.         DX(I + 1) = DY(I + 1)
  3514.         DY(I + 1) = DTEMP
  3515.         DTEMP = DX(I + 2)
  3516.         DX(I + 2) = DY(I + 2)
  3517.         DY(I + 2) = DTEMP
  3518.    50 CONTINUE
  3519.       RETURN
  3520.       END
  3521. LONGLEY DATA EXAMPLE FROM JASA 1967 62 819-841
  3522.  16  6 10 30  1  2
  3523.  83.0 234289 2356 1590 107608 1947 60323
  3524.  88.5 259426 2325 1456 108632 1948 61122
  3525.  88.2 258054 3682 1616 109773 1949 60171
  3526.  89.5 284599 3351 1650 110929 1950 61187
  3527.  96.2 328975 2099 3099 112075 1951 63221
  3528.  98.1 346999 1932 3594 113270 1952 63639
  3529.  99.0 365385 1870 3547 115094 1953 64989
  3530. 100.0 363112 3578 3350 116219 1954 63761
  3531. 101.2 397469 2904 3048 117388 1955 66019
  3532. 104.6 419180 2822 2857 118734 1956 67857
  3533. 108.4 442769 2936 2798 120445 1957 68169
  3534. 110.8 444546 4681 2637 121950 1958 66513
  3535. 112.6 482704 3813 2552 123366 1959 68655
  3536. 114.2 502601 3931 2514 125368 1960 69564
  3537. 115.7 518173 4806 2572 127852 1961 69331
  3538. 116.9 554894 4007 2827 130081 1962 70551
  3539.  16  6 21  0  3  2
  3540.  
  3541.