source/unres/src_MD/src/eigen.f

   1 C 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
   2 C 31 MAR 94 - MWS - ADD A VARIABLE TO END OF MACHSW COMMON
   3 C 26 JUN 93 - MWS - ETRED3: ADD RETURN FOR SPECIAL CASE N=1
   4 C  4 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
   5 C 30 AUG 91 - MWS - JACDIA: LIMIT ITERATIONS, USE EPSLON IN TEST.
   6 C 14 JUL 91 - MWS - JACOBI DIAGONALIZATION ALLOWS FOR LDVEC.NE.N
   7 C 29 JAN 91 - TLW - GLDIAG: CHANGED COMMON DIAGSW TO MACHSW
   8 C 29 OCT 90 - STE - FIX JACDIA UNDEFINED VARIABLE BUG
   9 C 14 SEP 90 - MK  - NEW JACOBI DIAGONALIZATION (KDIAG=3)
  10 C 27 MAR 88 - MWS - ALLOW FOR VECTOR ROUTINE IN GLDIAG
  11 C 11 AUG 87 - MWS - SANITIZE CONSTANTS IN EQLRAT
  12 C 15 FEB 87 - STE - FIX EINVIT SUB-MATRIX LOOP LIMIT
  13 C                   SCRATCH ARRAYS ARE N*8 REAL AND N INTEGER
  14 C  8 DEC 86 - STE - USE PERF INDEX FROM ESTPI1 TO JUDGE EINVIT FAILURE
  15 C 30 NOV 86 - STE - DELETE LIGENB, MAKE EVVRSP DEFAULT
  16 C                   (GIVEIS FAILS ON CRAY FOR BENCHMC AND BENCHCI)
  17 C  7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
  18 C 11 OCT 85 - STE - LIGENB,TQL2: USE DROT,DSWAP; TINVTB: SCALE VECTOR
  19 C                   BEFORE NORMALIZING; GENERIC FUNCTIONS
  20 C 24 FEB 84 - STE - INITIALIZE INDEX ARRAY FOR LIGENB IN GLDIAG
  21 C  1 DEC 83 - STE - CHANGE MACHEP FROM 2**-54 TO 2**-50
  22 C 28 SEP 82 - MWS - CONVERT TO IBM
  23 C
  24 C*MODULE EIGEN   *DECK EINVIT
  25       SUBROUTINE EINVIT(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
  26 C*
  27 C*    AUTHORS-
  28 C*       THIS IS A MODIFICATION OF TINVIT FROM EISPACK EDITION 3
  29 C*       DATED AUGUST 1983.
  30 C*       TINVIT IS A TRANSLATION OF THE INVERSE ITERATION TECHNIQUE
  31 C*       IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
  32 C*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
  33 C*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
  34 C*
  35 C*    PURPOSE -
  36 C*       THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
  37 C*       SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES.
  38 C*
  39 C*    METHOD -
  40 C*       INVERSE ITERATION.
  41 C*
  42 C*    ON ENTRY -
  43 C*       NM     - INTEGER
  44 C*                MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
  45 C*                ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
  46 C*                DIMENSION STATEMENT.
  47 C*       N      - INTEGER
  48 C*       D      - W.P. REAL (N)
  49 C*                CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
  50 C*       E      - W.P. REAL (N)
  51 C*                CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
  52 C*                IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
  53 C*       E2     - W.P. REAL (N)
  54 C*                CONTAINS THE SQUARES OF CORRESPONDING ELEMENTS OF E,
  55 C*                WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
  56 C*                E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
  57 C*                THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE
  58 C*                SUM OF THE MAGNITUDES OF D(I) AND D(I-1).  E2(1) MUST
  59 C*                CONTAIN 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER,
  60 C*                OR 2.0 IF THE EIGENVALUES ARE IN DESCENDING ORDER.
  61 C*                IF TQLRAT, BISECT, TRIDIB, OR IMTQLV
  62 C*                HAS BEEN USED TO FIND THE EIGENVALUES, THEIR
  63 C*                OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE.
  64 C*       M      - INTEGER
  65 C*                THE NUMBER OF SPECIFIED EIGENVECTORS.
  66 C*       W      - W.P. REAL (M)
  67 C*                CONTAINS THE M EIGENVALUES IN ASCENDING
  68 C*                OR DESCENDING ORDER.
  69 C*       IND    - INTEGER (M)
  70 C*                CONTAINS IN FIRST M POSITIONS THE SUBMATRIX INDICES
  71 C*                ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
  72 C*                1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX
  73 C*                FROM THE TOP, 2 FOR THOSE BELONGING TO THE SECOND
  74 C*                SUBMATRIX, ETC.
  75 C*       IERR   - INTEGER (LOGICAL UNIT NUMBER)
  76 C*                LOGICAL UNIT FOR ERROR MESSAGES
  77 C*
  78 C*    ON EXIT -
  79 C*       ALL INPUT ARRAYS ARE UNALTERED.
  80 C*       Z      - W.P. REAL (NM,M)
  81 C*                CONTAINS THE ASSOCIATED SET OF ORTHONORMAL
  82 C*                EIGENVECTORS. ANY VECTOR WHICH WHICH FAILS TO CONVERGE
  83 C*                IS LEFT AS IS (BUT NORMALIZED) WHEN ITERATING STOPPED.
  84 C*       IERR   - INTEGER
  85 C*                SET TO
  86 C*                ZERO    FOR NORMAL RETURN,
  87 C*                -R      IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
  88 C*                        EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS.
  89 C*                        (ONLY LAST FAILURE TO CONVERGE IS REPORTED)
  90 C*
  91 C*       RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
  92 C*
  93 C*       RV1    - W.P. REAL (N)
  94 C*                DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
  95 C*       RV2    - W.P. REAL (N)
  96 C*                SUPER(1)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
  97 C*       RV3    - W.P. REAL (N)
  98 C*                SUPER(2)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
  99 C*       RV4    - W.P. REAL (N)
 100 C*                ELEMENTS DEFINING L IN LU DECOMPOSITION
 101 C*       RV6    - W.P. REAL (N)
 102 C*                APPROXIMATE EIGENVECTOR
 103 C*
 104 C*    DIFFERENCES FROM EISPACK 3 -
 105 C*       EPS3 IS SCALED BY  EPSCAL  (ENHANCES CONVERGENCE, BUT
 106 C*          LOWERS ACCURACY)!
 107 C*       ONE MORE ITERATION (MINIMUM 2) IS PERFORMED AFTER CONVERGENCE
 108 C*          (ENHANCES ACCURACY)!
 109 C*       REPLACE LOOP WITH PYTHAG WITH SINGLE CALL TO DNRM2!
 110 C*       IF NOT CONVERGED, USE PERFORMANCE INDEX TO DECIDE ON ERROR
 111 C*          VALUE SETTING, BUT DO NOT STOP!
 112 C*       L.U. FOR ERROR MESSAGES PASSED THROUGH IERR
 113 C*       USE PARAMETER STATEMENTS AND GENERIC INTRINSIC FUNCTIONS
 114 C*       USE LEVEL 1 BLAS
 115 C*       USE IF-THEN-ELSE TO CLARIFY LOGIC
 116 C*       LOOP OVER SUBSPACES MADE INTO DO LOOP.
 117 C*       LOOP OVER INVERSE ITERATIONS MADE INTO DO LOOP
 118 C*       ZERO ONLY REQUIRED PORTIONS OF OUTPUT VECTOR
 119 C*
 120 C*    NOTE -
 121 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
 122 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
 123 C*
 124 C
 125       LOGICAL CONVGD,GOPARR,DSKWRK,MASWRK
 126 C
 127       INTEGER GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG
 128       INTEGER IND(M)
 129 C
 130       DOUBLE PRECISION D(N),E(N),E2(N),W(M),Z(NM,M)
 131       DOUBLE PRECISION RV1(N),RV2(N),RV3(N),RV4(N),RV6(N)
 132       DOUBLE PRECISION ANORM,EPS2,EPS3,EPS4,NORM,ORDER,RHO,U,UK,V
 133       DOUBLE PRECISION X0,X1,XU
 134       DOUBLE PRECISION EPSCAL,GRPTOL,HUNDRD,ONE,TEN,ZERO
 135       DOUBLE PRECISION EPSLON, ESTPI1, DASUM, DDOT, DNRM2
 136 C
 137       COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
 138 C
 139       PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00, GRPTOL = 0.001D+00)
 140       PARAMETER (EPSCAL = 0.5D+00, HUNDRD = 100.0D+00, TEN = 10.0D+00)
 141 C
 142   001 FORMAT(' EIGENVECTOR ROUTINE EINVIT DID NOT CONVERGE FOR VECTOR'
 143      *      ,I5,'.  NORM =',1P,E10.2,' PERFORMANCE INDEX =',E10.2/
 144      *      ' (AN ERROR HALT WILL OCCUR IF THE PI IS GREATER THAN 100)')
 145 C
 146 C-----------------------------------------------------------------------
 147 C
 148       LUEMSG = IERR
 149       IERR = 0
 150       X0 = ZERO
 151       UK = ZERO
 152       NORM = ZERO
 153       EPS2 = ZERO
 154       EPS3 = ZERO
 155       EPS4 = ZERO
 156       GROUP = 0
 157       TAG = 0
 158       ORDER = ONE - E2(1)
 159       Q = 0
 160       DO 930 SUBMAT = 1, N
 161          P = Q + 1
 162 C
 163 C        .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
 164 C
 165          DO 120 Q = P, N-1
 166             IF (E2(Q+1) .EQ. ZERO) GO TO 140
 167   120    CONTINUE
 168          Q = N
 169 C
 170 C        .......... FIND VECTORS BY INVERSE ITERATION ..........
 171 C
 172   140    CONTINUE
 173          TAG = TAG + 1
 174          ANORM = ZERO
 175          S = 0
 176 C
 177          DO 920 R = 1, M
 178             IF (IND(R) .NE. TAG) GO TO 920
 179             ITS = 1
 180             X1 = W(R)
 181             IF (S .NE. 0) GO TO 510
 182 C
 183 C        .......... CHECK FOR ISOLATED ROOT ..........
 184 C
 185             XU = ONE
 186             IF (P .EQ. Q) THEN
 187                RV6(P) = ONE
 188                CONVGD = .TRUE.
 189                GO TO 860
 190 C
 191             END IF
 192             NORM = ABS(D(P))
 193             DO 500 I = P+1, Q
 194                NORM = MAX( NORM, ABS(D(I)) + ABS(E(I)) )
 195   500       CONTINUE
 196 C
 197 C        .......... EPS2 IS THE CRITERION FOR GROUPING,
 198 C                   EPS3 REPLACES ZERO PIVOTS AND EQUAL
 199 C                   ROOTS ARE MODIFIED BY EPS3,
 200 C                   EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
 201 C
 202             EPS2 = GRPTOL * NORM
 203             EPS3 = EPSCAL * EPSLON(NORM)
 204             UK = Q - P + 1
 205             EPS4 = UK * EPS3
 206             UK = EPS4 / SQRT(UK)
 207             S = P
 208             GROUP = 0
 209             GO TO 520
 210 C
 211 C        .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
 212 C
 213   510       IF (ABS(X1-X0) .GE. EPS2) THEN
 214 C
 215 C                 ROOTS ARE SEPERATE
 216 C
 217                GROUP = 0
 218             ELSE
 219 C
 220 C                 ROOTS ARE CLOSE
 221 C
 222                GROUP = GROUP + 1
 223                IF (ORDER * (X1 - X0) .LE. EPS3) X1 = X0 + ORDER * EPS3
 224             END IF
 225 C
 226 C        .......... ELIMINATION WITH INTERCHANGES AND
 227 C                   INITIALIZATION OF VECTOR ..........
 228 C
 229   520       CONTINUE
 230 C
 231             U = D(P) - X1
 232             V = E(P+1)
 233             RV6(P) = UK
 234             DO 550 I = P+1, Q
 235                RV6(I) = UK
 236                IF (ABS(E(I)) .GT. ABS(U)) THEN
 237 C
 238 C                 EXCHANGE ROWS BEFORE ELIMINATION
 239 C
 240 C                  *** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
 241 C                      E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......
 242 C
 243                   XU = U / E(I)
 244                   RV4(I) = XU
 245                   RV1(I-1) = E(I)
 246                   RV2(I-1) = D(I) - X1
 247                   RV3(I-1) = E(I+1)
 248                   U = V - XU * RV2(I-1)
 249                   V = -XU * RV3(I-1)
 250 C
 251                ELSE
 252 C
 253 C                    STRAIGHT ELIMINATION
 254 C
 255                   XU = E(I) / U
 256                   RV4(I) = XU
 257                   RV1(I-1) = U
 258                   RV2(I-1) = V
 259                   RV3(I-1) = ZERO
 260                   U = D(I) - X1 - XU * V
 261                   V = E(I+1)
 262                END IF
 263   550       CONTINUE
 264 C
 265             IF (ABS(U) .LE. EPS3) U = EPS3
 266             RV1(Q) = U
 267             RV2(Q) = ZERO
 268             RV3(Q) = ZERO
 269 C
 270 C              DO INVERSE ITERATIONS
 271 C
 272             CONVGD = .FALSE.
 273             DO 800 ITS = 1, 5
 274                IF (ITS .EQ. 1) GO TO 600
 275 C
 276 C                    .......... FORWARD SUBSTITUTION ..........
 277 C
 278                   IF (NORM .EQ. ZERO) THEN
 279                      RV6(S) = EPS4
 280                      S = S + 1
 281                      IF (S .GT. Q) S = P
 282                   ELSE
 283                      XU = EPS4 / NORM
 284                      CALL DSCAL (Q-P+1, XU, RV6(P), 1)
 285                   END IF
 286 C
 287 C                     ... ELIMINATION OPERATIONS ON NEXT VECTOR
 288 C
 289                   DO 590 I = P+1, Q
 290                      U = RV6(I)
 291 C
 292 C                         IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
 293 C                         WAS PERFORMED EARLIER IN THE
 294 C                         TRIANGULARIZATION PROCESS ..........
 295 C
 296                      IF (RV1(I-1) .EQ. E(I)) THEN
 297                         U = RV6(I-1)
 298                         RV6(I-1) = RV6(I)
 299                      ELSE
 300                         U = RV6(I)
 301                      END IF
 302                      RV6(I) = U - RV4(I) * RV6(I-1)
 303   590             CONTINUE
 304   600          CONTINUE
 305 C
 306 C           .......... BACK SUBSTITUTION
 307 C
 308                RV6(Q) = RV6(Q) / RV1(Q)
 309                V = U
 310                U = RV6(Q)
 311                NORM = ABS(U)
 312                DO 620 I = Q-1, P, -1
 313                   RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
 314                   V = U
 315                   U = RV6(I)
 316                   NORM = NORM + ABS(U)
 317   620          CONTINUE
 318                IF (GROUP .EQ. 0) GO TO 700
 319 C
 320 C                 ....... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
 321 C                         MEMBERS OF GROUP ..........
 322 C
 323                   J = R
 324                   DO 680 JJ = 1, GROUP
 325   630                J = J - 1
 326                      IF (IND(J) .NE. TAG) GO TO 630
 327                      CALL DAXPY(Q-P+1, -DDOT(Q-P+1,RV6(P),1,Z(P,J),1),
 328      *                          Z(P,J),1,RV6(P),1)
 329   680             CONTINUE
 330                   NORM = DASUM(Q-P+1, RV6(P), 1)
 331   700          CONTINUE
 332 C
 333                IF (CONVGD) GO TO 840
 334                IF (NORM .GE. ONE) CONVGD = .TRUE.
 335   800       CONTINUE
 336 C
 337 C        .......... NORMALIZE SO THAT SUM OF SQUARES IS
 338 C                   1 AND EXPAND TO FULL ORDER ..........
 339 C
 340   840       CONTINUE
 341 C
 342             XU = ONE / DNRM2(Q-P+1,RV6(P),1)
 343 C
 344   860       CONTINUE
 345             DO 870 I = 1, P-1
 346                Z(I,R) = ZERO
 347   870       CONTINUE
 348             DO 890 I = P,Q
 349                Z(I,R) = RV6(I) * XU
 350   890       CONTINUE
 351             DO 900 I = Q+1, N
 352                Z(I,R) = ZERO
 353   900       CONTINUE
 354 C
 355             IF (.NOT.CONVGD) THEN
 356                RHO = ESTPI1(Q-P+1,X1,D(P),E(P),Z(P,R),ANORM)
 357                IF (RHO .GE. TEN .AND. LUEMSG .GT. 0 .AND. MASWRK)
 358      *             WRITE(LUEMSG,001) R,NORM,RHO
 359 C
 360 C               *** SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
 361 C
 362                IF (RHO .GT. HUNDRD) IERR = -R
 363             END IF
 364 C
 365             X0 = X1
 366   920    CONTINUE
 367 C
 368          IF (Q .EQ. N) GO TO 940
 369   930 CONTINUE
 370   940 CONTINUE
 371       RETURN
 372       END
 373 C*MODULE EIGEN   *DECK ELAUM
 374       SUBROUTINE ELAU(HINV,L,D,A,E)
 375 C
 376       DOUBLE PRECISION A(*)
 377       DOUBLE PRECISION D(L)
 378       DOUBLE PRECISION E(L)
 379       DOUBLE PRECISION F
 380       DOUBLE PRECISION G
 381       DOUBLE PRECISION HALF
 382       DOUBLE PRECISION HH
 383       DOUBLE PRECISION HINV
 384       DOUBLE PRECISION ZERO
 385 C
 386       PARAMETER (ZERO = 0.0D+00, HALF = 0.5D+00)
 387 C
 388       JL = L
 389       E(1) = A(1) * D(1)
 390       JK = 2
 391       DO 210 J = 2, JL
 392          F = D(J)
 393          G = ZERO
 394          JM1 = J - 1
 395 C
 396          DO 200 K = 1, JM1
 397             G = G + A(JK) * D(K)
 398             E(K) = E(K) + A(JK) * F
 399             JK = JK + 1
 400   200    CONTINUE
 401 C
 402          E(J) = G + A(JK) * F
 403          JK = JK + 1
 404   210 CONTINUE
 405 C
 406 C        .......... FORM P ..........
 407 C
 408       F = ZERO
 409       DO 245 J = 1, L
 410          E(J) = E(J) * HINV
 411          F = F + E(J) * D(J)
 412   245 CONTINUE
 413 C
 414 C     .......... FORM Q ..........
 415 C
 416       HH = F * HALF * HINV
 417       DO 250 J = 1, L
 418   250 E(J) = E(J) - HH * D(J)
 419 C
 420       RETURN
 421       END
 422 C*MODULE EIGEN   *DECK EPSLON
 423       DOUBLE PRECISION FUNCTION EPSLON (X)
 424 C*
 425 C*    AUTHORS -
 426 C*       THIS ROUTINE WAS TAKEN FROM EISPACK EDITION 3 DATED 4/6/83
 427 C*       THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE NOV 1986
 428 C*
 429 C*    PURPOSE -
 430 C*       ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
 431 C*
 432 C*    ON ENTRY -
 433 C*       X      - WORKING PRECISION REAL
 434 C*                VALUES TO FIND EPSLON FOR
 435 C*
 436 C*    ON EXIT -
 437 C*       EPSLON - WORKING PRECISION REAL
 438 C*                SMALLEST POSITIVE VALUE SUCH THAT X+EPSLON .NE. ZERO
 439 C*
 440 C*    QUALIFICATIONS -
 441 C*       THIS ROUTINE SHOULD PERFORM PROPERLY ON ALL SYSTEMS
 442 C*       SATISFYING THE FOLLOWING TWO ASSUMPTIONS,
 443 C*          1.  THE BASE USED IN REPRESENTING FLOATING POINT
 444 C*              NUMBERS IS NOT A POWER OF THREE.
 445 C*          2.  THE QUANTITY  A  IN STATEMENT 10 IS REPRESENTED TO
 446 C*              THE ACCURACY USED IN FLOATING POINT VARIABLES
 447 C*              THAT ARE STORED IN MEMORY.
 448 C*       THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO
 449 C*       FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING
 450 C*       ASSUMPTION 2.
 451 C*       UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT,
 452 C*              A  IS NOT EXACTLY EQUAL TO FOUR-THIRDS,
 453 C*              B  HAS A ZERO FOR ITS LAST BIT OR DIGIT,
 454 C*              C  IS NOT EXACTLY EQUAL TO ONE,
 455 C*              EPS  MEASURES THE SEPARATION OF 1.0 FROM
 456 C*                   THE NEXT LARGER FLOATING POINT NUMBER.
 457 C*       THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED
 458 C*       ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD.
 459 C*
 460 C*    DIFFERENCES FROM EISPACK 3 -
 461 C*       USE IS MADE OF PARAMETER STATEMENTS AND INTRINSIC FUNCTIONS
 462 C*       --NO EXECUTEABLE CODE CHANGES--
 463 C*
 464 C*    NOTE -
 465 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
 466 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
 467 C
 468       DOUBLE PRECISION A,B,C,EPS,X
 469       DOUBLE PRECISION ZERO, ONE, THREE, FOUR
 470 C
 471       PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, FOUR=4.0D+00)
 472 C
 473 C-----------------------------------------------------------------------
 474 C
 475       A = FOUR/THREE
 476    10 B = A - ONE
 477       C = B + B + B
 478       EPS = ABS(C - ONE)
 479       IF (EPS .EQ. ZERO) GO TO 10
 480       EPSLON = EPS*ABS(X)
 481       RETURN
 482       END
 483 C*MODULE EIGEN   *DECK EQLRAT
 484       SUBROUTINE EQLRAT(N,DIAG,E,E2IN,D,IND,IERR,E2)
 485 C*
 486 C*    AUTHORS -
 487 C*       THIS IS A MODIFICATION OF ROUTINE EQLRAT FROM EISPACK EDITION 3
 488 C*       DATED AUGUST 1983.
 489 C*       TQLRAT IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
 490 C*       ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
 491 C*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
 492 C*
 493 C*    PURPOSE -
 494 C*       THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
 495 C*       TRIDIAGONAL MATRIX
 496 C*
 497 C*    METHOD -
 498 C*       RATIONAL QL
 499 C*
 500 C*    ON ENTRY -
 501 C*       N      - INTEGER
 502 C*                THE ORDER OF THE MATRIX.
 503 C*       D      - W.P. REAL (N)
 504 C*                CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
 505 C*       E2     - W.P. REAL (N)
 506 C*                CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF
 507 C*                THE INPUT MATRIX IN ITS LAST N-1 POSITIONS.
 508 C*                E2(1) IS ARBITRARY.
 509 C*
 510 C*     ON EXIT -
 511 C*       D      - W.P. REAL (N)
 512 C*                CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
 513 C*                ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
 514 C*                ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
 515 C*                THE SMALLEST EIGENVALUES.
 516 C*       E2     - W.P. REAL (N)
 517 C*                DESTROYED.
 518 C*       IERR   - INTEGER
 519 C*                SET TO
 520 C*                ZERO       FOR NORMAL RETURN,
 521 C*                J          IF THE J-TH EIGENVALUE HAS NOT BEEN
 522 C*                           DETERMINED AFTER 30 ITERATIONS.
 523 C*
 524 C*    DIFFERENCES FROM EISPACK 3 -
 525 C*       G=G+B INSTEAD OF IF(G.EQ.0) G=B ; B=B/4
 526 C*       F77 BACKWARD LOOPS INSTEAD OF F66 CONSTRUCT
 527 C*       GENERIC INTRINSIC FUNCTIONS
 528 C*       ARRARY  IND  ADDED FOR USE BY EINVIT
 529 C*
 530 C*    NOTE -
 531 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
 532 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
 533 C
 534       INTEGER I,J,L,M,N,II,L1,IERR
 535       INTEGER IND(N)
 536 C
 537       DOUBLE PRECISION D(N),E(N),E2(N),DIAG(N),E2IN(N)
 538       DOUBLE PRECISION B,C,F,G,H,P,R,S,T,EPSLON
 539       DOUBLE PRECISION SCALE,ZERO,ONE
 540 C
 541       PARAMETER (ZERO = 0.0D+00, SCALE= 1.0D+00/64.0D+00, ONE = 1.0D+00)
 542 C
 543 C-----------------------------------------------------------------------
 544       IERR = 0
 545       D(1)=DIAG(1)
 546       IND(1) = 1
 547       K = 0
 548       ITAG = 0
 549       IF (N .EQ. 1) GO TO 1001
 550 C
 551       DO 100 I = 2, N
 552          D(I)=DIAG(I)
 553   100 E2(I-1) = E2IN(I)
 554 C
 555       F = ZERO
 556       T = ZERO
 557       B = EPSLON(ONE)
 558       C = B *B
 559       B = B * SCALE
 560       E2(N) = ZERO
 561 C
 562       DO 290 L = 1, N
 563          H = ABS(D(L)) + ABS(E(L))
 564          IF (T .GE. H) GO TO 105
 565             T = H
 566             B = EPSLON(T)
 567             C = B * B
 568             B = B * SCALE
 569   105    CONTINUE
 570 C     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
 571          M = L - 1
 572   110    M = M + 1
 573          IF (E2(M) .GT. C) GO TO 110
 574 C     .......... E2(N) IS ALWAYS ZERO, SO THERE IS AN EXIT
 575 C                FROM THE LOOP ..........
 576 C
 577          IF (M .LE. K) GO TO 125
 578             IF (M .NE. N) E2IN(M+1) = ZERO
 579             K = M
 580             ITAG = ITAG + 1
 581   125    CONTINUE
 582          IF (M .EQ. L) GO TO 210
 583 C
 584 C           ITERATE
 585 C
 586          DO 205 J = 1, 30
 587 C              .......... FORM SHIFT ..........
 588             L1 = L + 1
 589             S = SQRT(E2(L))
 590             G = D(L)
 591             P = (D(L1) - G) / (2.0D+00 * S)
 592             R = SQRT(P*P+1.0D+00)
 593             D(L) = S / (P + SIGN(R,P))
 594             H = G - D(L)
 595 C
 596             DO 140 I = L1, N
 597   140       D(I) = D(I) - H
 598 C
 599             F = F + H
 600 C              .......... RATIONAL QL TRANSFORMATION ..........
 601             G = D(M) + B
 602             H = G
 603             S = ZERO
 604             DO 200 I = M-1,L,-1
 605                P = G * H
 606                R = P + E2(I)
 607                E2(I+1) = S * R
 608                S = E2(I) / R
 609                D(I+1) = H + S * (H + D(I))
 610                G = D(I) - E2(I) / G   + B
 611                H = G * P / R
 612   200       CONTINUE
 613 C
 614             E2(L) = S * G
 615             D(L) = H
 616 C              .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST
 617             IF (H .EQ. ZERO) GO TO 210
 618             IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
 619             E2(L) = H * E2(L)
 620             IF (E2(L) .EQ. ZERO) GO TO 210
 621   205    CONTINUE
 622 C     .......... SET ERROR -- NO CONVERGENCE TO AN
 623 C                EIGENVALUE AFTER 30 ITERATIONS ..........
 624       IERR = L
 625       GO TO 1001
 626 C
 627 C           CONVERGED
 628 C
 629   210    P = D(L) + F
 630 C           .......... ORDER EIGENVALUES ..........
 631          I = 1
 632          IF (L .EQ. 1) GO TO 250
 633             IF (P .LT. D(1)) GO TO 230
 634                I = L
 635 C           .......... LOOP TO FIND ORDERED POSITION
 636   220          I = I - 1
 637                IF (P .LT. D(I)) GO TO 220
 638 C
 639                I = I + 1
 640                IF (I .EQ. L) GO TO 250
 641   230       CONTINUE
 642             DO 240 II = L, I+1, -1
 643                D(II) = D(II-1)
 644                IND(II) = IND(II-1)
 645   240       CONTINUE
 646 C
 647   250    CONTINUE
 648          D(I) = P
 649          IND(I) = ITAG
 650   290 CONTINUE
 651 C
 652  1001 RETURN
 653       END
 654 C*MODULE EIGEN   *DECK ESTPI1
 655       DOUBLE PRECISION FUNCTION ESTPI1 (N,EVAL,D,E,X,ANORM)
 656 C*
 657 C*    AUTHOR -
 658 C*       STEPHEN T. ELBERT (AMES LABORATORY-USDOE) DATE: 5 DEC 1986
 659 C*
 660 C*    PURPOSE -
 661 C*       EVALUATE SYMMETRIC TRIDIAGONAL MATRIX PERFORMANCE INDEX
 662 C*       *        *         *                  *           *
 663 C*       FOR 1 EIGENVECTOR
 664 C*           *
 665 C*
 666 C*    METHOD -
 667 C*       THIS ROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX A*X-X*EVAL
 668 C*       WHERE  A  IS A SYMMETRIC TRIDIAGONAL MATRIX STORED
 669 C*       IN THE DIAGONAL (D) AND SUB-DIAGONAL (E) VECTORS, EVAL IS THE
 670 C*       EIGENVALUE OF AN EIGENVECTOR OF  A,  NAMELY  X.
 671 C*       THIS NORM IS SCALED BY MACHINE ACCURACY FOR THE PROBLEM SIZE.
 672 C*       ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
 673 C*
 674 C*    ON ENTRY -
 675 C*       N      - INTEGER
 676 C*                THE ORDER OF THE MATRIX  A.
 677 C*       EVAL   - W.P. REAL
 678 C*                THE EIGENVALUE CORRESPONDING TO VECTOR  X.
 679 C*       D      - W.P. REAL (N)
 680 C*                THE DIAGONAL VECTOR OF  A.
 681 C*       E      - W.P. REAL (N)
 682 C*                THE SUB-DIAGONAL VECTOR OF  A.
 683 C*       X      - W.P. REAL (N)
 684 C*                AN EIGENVECTOR OF  A.
 685 C*       ANORM  - W.P. REAL
 686 C*                THE NORM OF  A  IF IT HAS BEEN PREVIOUSLY COMPUTED.
 687 C*
 688 C*    ON EXIT -
 689 C*       ANORM  - W.P. REAL
 690 C*                THE NORM OF  A, COMPUTED IF INITIALLY ZERO.
 691 C*       ESTPI1 - W.P. REAL
 692 C*          !!A*X-X*EVAL!! / (EPSLON(10*N)*!!A!!*!!X!!);
 693 C*          WHERE EPSLON(X) IS THE SMALLEST NUMBER SUCH THAT
 694 C*             X + EPSLON(X) .NE. X
 695 C*
 696 C*          ESTPI1 .LT. 1 == SATISFACTORY PERFORMANCE
 697 C*                 .GE. 1 AND .LE. 100 == MARGINAL PERFORMANCE
 698 C*                 .GT. 100 == POOR PERFORMANCE
 699 C*          (SEE LECT. NOTES IN COMP. SCI. VOL.6 PP 124-125)
 700 C
 701       DOUBLE PRECISION ANORM,EVAL,RNORM,SIZE,XNORM
 702       DOUBLE PRECISION D(N), E(N), X(N)
 703       DOUBLE PRECISION EPSLON, ONE, ZERO
 704 C
 705       PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00)
 706 C
 707 C-----------------------------------------------------------------------
 708 C
 709       ESTPI1 = ZERO
 710       IF( N .LE. 1 ) RETURN
 711       SIZE = 10 * N
 712       IF (ANORM .EQ. ZERO) THEN
 713 C
 714 C              COMPUTE NORM OF  A
 715 C
 716          ANORM = MAX( ABS(D(1)) + ABS(E(2))
 717      *               ,ABS(D(N)) + ABS(E(N)))
 718          DO 110 I = 2, N-1
 719             ANORM = MAX( ANORM, ABS(E(I))+ABS(D(I))+ABS(E(I+1)))
 720   110    CONTINUE
 721          IF(ANORM .EQ. ZERO) ANORM = ONE
 722       END IF
 723 C
 724 C           COMPUTE NORMS OF RESIDUAL AND EIGENVECTOR
 725 C
 726       XNORM = ABS(X(1)) + ABS(X(N))
 727       RNORM = ABS( (D(1)-EVAL)*X(1) + E(2)*X(2))
 728      *       +ABS( (D(N)-EVAL)*X(N) + E(N)*X(N-1))
 729       DO 120 I = 2, N-1
 730          XNORM = XNORM + ABS(X(I))
 731          RNORM = RNORM + ABS(E(I)*X(I-1) + (D(I)-EVAL)*X(I)
 732      *                       + E(I+1)*X(I+1))
 733   120 CONTINUE
 734 C
 735       ESTPI1 = RNORM / (EPSLON(SIZE)*ANORM*XNORM)
 736       RETURN
 737       END
 738 C*MODULE EIGEN   *DECK ETRBK3
 739       SUBROUTINE ETRBK3(NM,N,NV,A,M,Z)
 740 C*
 741 C*    AUTHORS-
 742 C*       THIS IS A MODIFICATION OF ROUTINE TRBAK3 FROM EISPACK EDITION 3
 743 C*       DATED AUGUST 1983.
 744 C*       EISPACK TRBAK3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
 745 C*       NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
 746 C*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
 747 C*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
 748 C*
 749 C*    PURPOSE -
 750 C*       THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
 751 C*       MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
 752 C*       SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  ETRED3.
 753 C*
 754 C*    METHOD -
 755 C*       THE CALCULATION IS CARRIED OUT BY FORMING THE MATRIX PRODUCT
 756 C*          Q*Z
 757 C*       WHERE  Q  IS A PRODUCT OF THE ORTHOGONAL SYMMETRIC MATRICES
 758 C*                Q = PROD(I)[1 - U(I)*.TRANSPOSE.U(I)*H(I)]
 759 C*       U  IS THE AUGMENTED SUB-DIAGONAL ROWS OF  A  AND
 760 C*       Z  IS THE SET OF EIGENVECTORS OF THE TRIDIAGONAL
 761 C*       MATRIX  F  WHICH WAS FORMED FROM THE ORIGINAL SYMMETRIC
 762 C*       MATRIX  C  BY THE SIMILARITY TRANSFORMATION
 763 C*                F = Q(TRANSPOSE) C Q
 764 C*       NOTE THAT ETRBK3 PRESERVES VECTOR EUCLIDEAN NORMS.
 765 C*
 766 C*
 767 C*    COMPLEXITY -
 768 C*       M*N**2
 769 C*
 770 C*    ON ENTRY-
 771 C*       NM     - INTEGER
 772 C*                MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
 773 C*                ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
 774 C*                DIMENSION STATEMENT.
 775 C*       N      - INTEGER
 776 C*                THE ORDER OF THE MATRIX  A.
 777 C*       NV     - INTEGER
 778 C*                MUST BE SET TO THE DIMENSION OF THE ARRAY  A  AS
 779 C*                DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT.
 780 C*       A      - W.P. REAL (NV)
 781 C*                CONTAINS INFORMATION ABOUT THE ORTHOGONAL
 782 C*                TRANSFORMATIONS USED IN THE REDUCTION BY  ETRED3  IN
 783 C*                ITS FIRST  NV = N*(N+1)/2 POSITIONS.
 784 C*       M      - INTEGER
 785 C*                THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
 786 C*       Z      - W.P REAL (NM,M)
 787 C*                CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
 788 C*                IN ITS FIRST M COLUMNS.
 789 C*
 790 C*    ON EXIT-
 791 C*       Z      - W.P. REAL (NM,M)
 792 C*                CONTAINS THE TRANSFORMED EIGENVECTORS
 793 C*                IN ITS FIRST M COLUMNS.
 794 C*
 795 C*    DIFFERENCES WITH EISPACK 3 -
 796 C*       THE TWO INNER LOOPS ARE REPLACED BY DDOT AND DAXPY.
 797 C*       MULTIPLICATION USED INSTEAD OF DIVISION TO FIND S.
 798 C*       OUTER LOOP RANGE CHANGED FROM 2,N TO 3,N.
 799 C*       ADDRESS POINTERS FOR  A  SIMPLIFIED.
 800 C*
 801 C*    NOTE -
 802 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
 803 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
 804 C
 805       INTEGER I,II,IM1,IZ,J,M,N,NM,NV
 806 C
 807       DOUBLE PRECISION A(NV),Z(NM,M)
 808       DOUBLE PRECISION H,S,DDOT,ZERO
 809 C
 810       PARAMETER (ZERO = 0.0D+00)
 811 C
 812 C-----------------------------------------------------------------------
 813 C
 814       IF (M .EQ. 0) RETURN
 815       IF (N .LE. 2) RETURN
 816 C
 817       II=3
 818       DO 140 I = 3, N
 819          IZ=II+1
 820          II=II+I
 821          H = A(II)
 822          IF (H .EQ. ZERO) GO TO 140
 823             IM1 = I - 1
 824             DO 130 J = 1, M
 825                S = -( DDOT(IM1,A(IZ),1,Z(1,J),1) * H) * H
 826                CALL DAXPY(IM1,S,A(IZ),1,Z(1,J),1)
 827   130       CONTINUE
 828   140 CONTINUE
 829       RETURN
 830       END
 831 C*MODULE EIGEN   *DECK ETRED3
 832       SUBROUTINE ETRED3(N,NV,A,D,E,E2)
 833 C*
 834 C*    AUTHORS -
 835 C*       THIS IS A MODIFICATION OF ROUTINE TRED3 FROM EISPACK EDITION 3
 836 C*       DATED AUGUST 1983.
 837 C*       EISPACK TRED3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
 838 C*       NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
 839 C*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
 840 C*       THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE JUN 1986
 841 C*
 842 C*    PURPOSE -
 843 C*       THIS ROUTINE REDUCES A REAL SYMMETRIC (PACKED) MATRIX, STORED
 844 C*       AS A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
 845 C*       USING ORTHOGONAL SIMILARITY TRANSFORMATIONS, PRESERVING THE
 846 C*       INFORMATION ABOUT THE TRANSFORMATIONS IN  A.
 847 C*
 848 C*    METHOD -
 849 C*       THE TRIDIAGONAL REDUCTION IS PERFORMED IN THE FOLLOWING WAY.
 850 C*       STARTING WITH J=N, THE ELEMENTS IN THE J-TH ROW TO THE
 851 C*       LEFT OF THE DIAGONAL ARE FIRST SCALED, TO AVOID POSSIBLE
 852 C*       UNDERFLOW IN THE TRANSFORMATION THAT MIGHT RESULT IN SEVERE
 853 C*       DEPARTURE FROM ORTHOGONALITY.  THE SUM OF SQUARES  SIGMA  OF
 854 C*       THESE SCALED ELEMENTS IS NEXT FORMED.  THEN, A VECTOR  U  AND
 855 C*       A SCALAR
 856 C*                      H = U(TRANSPOSE) * U / 2
 857 C*       DEFINE A REFLECTION OPERATOR
 858 C*                      P = I - U * U(TRANSPOSE) / H
 859 C*       WHICH IS ORTHOGONAL AND SYMMETRIC AND FOR WHICH THE
 860 C*       SIMILIARITY TRANSFORMATION  PAP  ELIMINATES THE ELEMENTS IN
 861 C*       THE J-TH ROW OF  A  TO THE LEFT OF THE SUBDIAGONAL AND THE
 862 C*       SYMMETRICAL ELEMENTS IN THE J-TH COLUMN.
 863 C*
 864 C*       THE NON-ZERO COMPONENTS OF  U  ARE THE ELEMENTS OF THE J-TH
 865 C*       ROW TO THE LEFT OF THE DIAGONAL WITH THE LAST OF THEM
 866 C*       AUGMENTED BY THE SQUARE ROOT OF  SIGMA  PREFIXED BY THE SIGN
 867 C*       OF THE SUBDIAGONAL ELEMENT.  BY STORING THE TRANSFORMED SUB-
 868 C*       DIAGONAL ELEMENT IN  E(J)  AND NOT OVERWRITING THE ROW
 869 C*       ELEMENTS ELIMINATED IN THE TRANSFORMATION, FULL INFORMATION
 870 C*       ABOUT  P  IS SAVE FOR LATER USE IN  ETRBK3.
 871 C*
 872 C*       THE TRANSFORMATION SETS  E2(J)  EQUAL TO  SIGMA  AND  E(J)
 873 C*       EQUAL TO THE SQUARE ROOT OF  SIGMA  PREFIXED BY THE SIGN
 874 C*       OF THE REPLACED SUBDIAGONAL ELEMENT.
 875 C*
 876 C*       THE ABOVE STEPS ARE REPEATED ON FURTHER ROWS OF THE
 877 C*       TRANSFORMED  A  IN REVERSE ORDER UNTIL  A  IS REDUCED TO TRI-
 878 C*       DIAGONAL FORM, THAT IS, REPEATED FOR  J = N-1,N-2,...,3.
 879 C*
 880 C*    COMPLEXITY -
 881 C*       2/3 N**3
 882 C*
 883 C*    ON ENTRY-
 884 C*       N      - INTEGER
 885 C*                THE ORDER OF THE MATRIX.
 886 C*       NV     - INTEGER
 887 C*                MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
 888 C*                AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT
 889 C*       A      - W.P. REAL (NV)
 890 C*                CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
 891 C*                INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
 892 C*                ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
 893 C*
 894 C*    ON EXIT-
 895 C*       A      - W.P. REAL (NV)
 896 C*                CONTAINS INFORMATION ABOUT THE ORTHOGONAL
 897 C*                TRANSFORMATIONS USED IN THE REDUCTION.
 898 C*       D      - W.P. REAL (N)
 899 C*                CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL
 900 C*                MATRIX.
 901 C*       E      - W.P. REAL (N)
 902 C*                CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
 903 C*                MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO
 904 C*       E2     - W.P. REAL (N)
 905 C*                CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF
 906 C*                E. MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
 907 C*
 908 C*    DIFFERENCES FROM EISPACK 3 -
 909 C*       OUTER LOOP CHANGED FROM II=1,N TO I=N,3,-1
 910 C*       PARAMETER STATEMENT AND GENERIC INTRINSIC FUNCTIONS USED
 911 C*       SCALE.NE.0 TEST NOW SPOTS TRI-DIAGONAL FORM
 912 C*       VALUES LESS THAN EPSLON CLEARED TO ZERO
 913 C*       USE BLAS(1)
 914 C*       U NOT COPIED TO D, LEFT IN A
 915 C*       E2 COMPUTED FROM E
 916 C*       INNER LOOPS SPLIT INTO ROUTINES ELAU AND FREDA
 917 C*       INVERSE OF H STORED INSTEAD OF H
 918 C*
 919 C*    NOTE -
 920 C*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
 921 C*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
 922 C
 923       INTEGER I,IIA,IZ0,L,N,NV
 924 C
 925       DOUBLE PRECISION A(NV),D(N),E(N),E2(N)
 926       DOUBLE PRECISION AIIMAX,F,G,H,HROOT,SCALE,SCALEI
 927       DOUBLE PRECISION DASUM, DNRM2
 928       DOUBLE PRECISION ONE, ZERO
 929 C
 930       PARAMETER (ZERO = 0.0D+00, ONE = 1.0D+00)
 931 C
 932 C-----------------------------------------------------------------------
 933 C
 934       IF (N .LE. 2) GO TO 310
 935       IZ0 = (N*N+N)/2
 936       AIIMAX = ABS(A(IZ0))
 937       DO 300 I = N, 3, -1
 938          L = I - 1
 939          IIA = IZ0
 940          IZ0 = IZ0 - I
 941          AIIMAX = MAX(AIIMAX, ABS(A(IIA)))
 942          SCALE = DASUM (L, A(IZ0+1), 1)
 943          IF(SCALE .EQ. ABS(A(IIA-1)) .OR. AIIMAX+SCALE .EQ. AIIMAX) THEN
 944 C
 945 C           THIS ROW IS ALREADY IN TRI-DIAGONAL FORM
 946 C
 947             D(I) = A(IIA)
 948             IF (AIIMAX+D(I) .EQ. AIIMAX) D(I) = ZERO
 949             E(I) = A(IIA-1)
 950             IF (AIIMAX+E(I) .EQ. AIIMAX) E(I) = ZERO
 951             E2(I) = E(I)*E(I)
 952             A(IIA) = ZERO
 953             GO TO 300
 954 C
 955          END IF
 956 C
 957          SCALEI = ONE / SCALE
 958          CALL DSCAL(L,SCALEI,A(IZ0+1),1)
 959          HROOT = DNRM2(L,A(IZ0+1),1)
 960 C
 961          F = A(IZ0+L)
 962          G = -SIGN(HROOT,F)
 963          E(I) = SCALE * G
 964          E2(I) = E(I)*E(I)
 965          H = HROOT*HROOT - F * G
 966          A(IZ0+L) = F - G
 967          D(I) = A(IIA)
 968          A(IIA) = ONE / SQRT(H)
 969 C           .......... FORM P THEN Q IN E(1:L) ..........
 970          CALL ELAU(ONE/H,L,A(IZ0+1),A,E)
 971 C           .......... FORM REDUCED A ..........
 972          CALL FREDA(L,A(IZ0+1),A,E)
 973 C
 974   300 CONTINUE
 975   310 CONTINUE
 976       E(1) = ZERO
 977       E2(1)= ZERO
 978       D(1) = A(1)
 979       IF(N.EQ.1) RETURN
 980 C
 981       E(2) = A(2)
 982       E2(2)= A(2)*A(2)
 983       D(2) = A(3)
 984       RETURN
 985       END
 986 C*MODULE EIGEN   *DECK EVVRSP
 987       SUBROUTINE EVVRSP(MSGFL,N,NVECT,LENA,NV,A,B,IND,ROOT,
 988      *                  VECT,IORDER,IERR)
 989 C*
 990 C*    AUTHOR:  S. T. ELBERT, AMES LABORATORY-USDOE, JUNE 1985
 991 C*
 992 C*    PURPOSE -
 993 C*       FINDS   (ALL) EIGENVALUES    AND    (SOME OR ALL) EIGENVECTORS
 994 C*                     *    *                                   *
 995 C*       OF A REAL SYMMETRIC PACKED MATRIX.
 996 C*            *    *         *
 997 C*
 998 C*    METHOD -
 999 C*       THE METHOD AS PRESENTED IN THIS ROUTINE CONSISTS OF FOUR STEPS:
1000 C*       FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE
1001 C*       HOUSEHOLDER TECHNIQUE (ORTHOGONAL SIMILARITY TRANSFORMATIONS).
1002 C*       SECOND, THE ROOTS ARE LOCATED USING THE RATIONAL QL METHOD.
1003 C*       THIRD, THE VECTORS OF THE TRIDIAGONAL FORM ARE EVALUATED BY THE
1004 C*       INVERSE ITERATION TECHNIQUE.  VECTORS FOR DEGENERATE OR NEAR-
1005 C*       DEGENERATE ROOTS ARE FORCED TO BE ORTHOGONAL.
1006 C*       FOURTH, THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE
1007 C*       ORIGINAL ARRAY.
1008 C*
1009 C*       THESE ROUTINES ARE MODIFICATIONS OF THE EISPACK 3
1010 C*       ROUTINES TRED3, TQLRAT, TINVIT AND TRBAK3
1011 C*
1012 C*       FOR FURTHER DETAILS, SEE EISPACK USERS GUIDE, B. T. SMITH
1013 C*       ET AL, SPRINGER-VERLAG, LECTURE NOTES IN COMPUTER SCIENCE,
1014 C*       VOL. 6, 2-ND EDITION, 1976.  ANOTHER GOOD REFERENCE IS
1015 C*       THE SYMMETRIC EIGENVALUE PROBLEM BY B. N. PARLETT
1016 C*       PUBLISHED BY PRENTICE-HALL, INC., ENGLEWOOD CLIFFS, N.J. (1980)
1017 C*
1018 C*    ON ENTRY -
1019 C*       MSGFL  - INTEGER (LOGICAL UNIT NO.)
1020 C*                FILE WHERE ERROR MESSAGES WILL BE PRINTED.
1021 C*                IF MSGFL IS 0, ERROR MESSAGES WILL BE PRINTED ON LU 6.
1022 C*                IF MSGFL IS NEGATIVE, NO ERROR MESSAGES PRINTED.
1023 C*       N      - INTEGER
1024 C*                ORDER OF MATRIX A.
1025 C*       NVECT  - INTEGER
1026 C*                NUMBER OF VECTORS DESIRED.  0 .LE. NVECT .LE. N.
1027 C*       LENA   - INTEGER
1028 C*                DIMENSION OF  A  IN CALLING ROUTINE.  MUST NOT BE LESS
1029 C*                THAN (N*N+N)/2.
1030 C*       NV     - INTEGER
1031 C*                ROW DIMENSION OF VECT IN CALLING ROUTINE.   N .LE. NV.
1032 C*       A      - WORKING PRECISION REAL (LENA)
1033 C*                INPUT MATRIX, ROWS OF THE LOWER TRIANGLE PACKED INTO
1034 C*                LINEAR ARRAY OF DIMENSION N*(N+1)/2.  THE PACKED ORDER
1035 C*                IS A(1,1), A(2,1), A(2,2), A(3,1), A(3,2), ...
1036 C*       B      - WORKING PRECISION REAL (N,8)
1037 C*                SCRATCH ARRAY, 8*N ELEMENTS
1038 C*       IND    - INTEGER (N)
1039 C*                SCRATCH ARRAY OF LENGTH N.
1040 C*       IORDER - INTEGER
1041 C*                ROOT ORDERING FLAG.
1042 C*                = 0, ROOTS WILL BE PUT IN ASCENDING ORDER.
1043 C*                = 2, ROOTS WILL BE PUT IN DESCENDING ORDER.
1044 C*
1045 C*    ON EXIT -
1046 C*       A      - DESTORYED.  NOW HOLDS REFLECTION OPERATORS.
1047 C*       ROOT   - WORKING PRECISION REAL (N)
1048 C*                ALL EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
1049 C*                  IF IORDER = 0, ROOT(1) .LE. ... .LE. ROOT(N)
1050 C*                  IF IORDER = 2, ROOT(1) .GE. ... .GE. ROOT(N)
1051 C*       VECT   - WORKING PRECISION REAL (NV,NVECT)
1052 C*                EIGENVECTORS FOR ROOT(1), ..., ROOT(NVECT).
1053 C*       IERR   - INTEGER
1054 C*                = 0 IF NO ERROR DETECTED,
1055 C*                = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
1056 C*                = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
1057 C*                (FAILURES SHOULD BE VERY RARE.  CONTACT C. MOLER.)
1058 C*
1059 C
1060       LOGICAL GOPARR,DSKWRK,MASWRK
1061 C
1062       DOUBLE PRECISION A(LENA)
1063       DOUBLE PRECISION B(N,8)
1064       DOUBLE PRECISION ROOT(N)
1065       DOUBLE PRECISION T
1066       DOUBLE PRECISION VECT(NV,*)
1067 C
1068       INTEGER IND(N)
1069 C
1070       COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1071 C
1072   900 FORMAT(26H0*** EVVRSP PARAMETERS ***/
1073      +       14H ***      N = ,I8,4H ***/
1074      +       14H ***  NVECT = ,I8,4H ***/
1075      +       14H ***   LENA = ,I8,4H ***/
1076      +       14H ***     NV = ,I8,4H ***/
1077      +       14H *** IORDER = ,I8,4H ***/
1078      +       14H ***   IERR = ,I8,4H ***)
1079   901 FORMAT(37H VALUE OF LENA IS LESS THAN (N*N+N)/2)
1080   902 FORMAT(39H EQLRAT HAS FAILED TO CONVERGE FOR ROOT,I5)
1081   903 FORMAT(18H NV IS LESS THAN N)
1082   904 FORMAT(41H EINVIT HAS FAILED TO CONVERGE FOR VECTOR,I5)
1083   905 FORMAT(51H VALUE OF IORDER MUST BE 0 (SMALLEST ROOT FIRST) OR
1084      *      ,23H 2 (LARGEST ROOT FIRST))
1085   906 FORMAT(' VALUE OF N IS LESS THAN OR EQUAL ZERO')
1086 C
1087 C-----------------------------------------------------------------------
1088 C
1089       LMSGFL=MSGFL
1090       IF (MSGFL .EQ. 0) LMSGFL=6
1091       IERR = N - 1
1092       IF (N .LE. 0) GO TO 800
1093       IERR = N + 1
1094       IF ( (N*N+N)/2 .GT. LENA) GO TO 810
1095 C
1096 C        REDUCE REAL SYMMETRIC MATRIX A TO TRIDIAGONAL FORM
1097 C
1098       CALL ETRED3(N,LENA,A,B(1,1),B(1,2),B(1,3))
1099 C
1100 C        FIND ALL EIGENVALUES OF TRIDIAGONAL MATRIX
1101 C
1102       CALL EQLRAT(N,B(1,1),B(1,2),B(1,3),ROOT,IND,IERR,B(1,4))
1103       IF (IERR .NE. 0) GO TO 820
1104 C
1105 C         CHECK THE DESIRED ORDER OF THE EIGENVALUES
1106 C
1107       B(1,3) = IORDER
1108       IF (IORDER .EQ. 0) GO TO 300
1109          IF (IORDER .NE. 2) GO TO 850
1110 C
1111 C         ORDER ROOTS IN DESCENDING ORDER (LARGEST FIRST)...
1112 C        TURN ROOT AND IND ARRAYS END FOR END
1113 C
1114          DO 210 I = 1, N/2
1115             J = N+1-I
1116             T = ROOT(I)
1117             ROOT(I) = ROOT(J)
1118             ROOT(J) = T
1119             L = IND(I)
1120             IND(I) = IND(J)
1121             IND(J) = L
1122   210    CONTINUE
1123 C
1124 C           FIND I AND J MARKING THE START AND END OF A SEQUENCE
1125 C           OF DEGENERATE ROOTS
1126 C
1127          I=0
1128   220    CONTINUE
1129             I = I+1
1130             IF (I .GT. N) GO TO 300
1131             DO 230 J=I,N
1132                IF (ROOT(J) .NE. ROOT(I)) GO TO 240
1133   230       CONTINUE
1134             J = N+1
1135   240       CONTINUE
1136             J = J-1
1137             IF (J .EQ. I) GO TO 220
1138 C
1139 C                    TURN AROUND IND BETWEEN I AND J
1140 C
1141             JSV = J
1142             KLIM = (J-I+1)/2
1143             DO 250 K=1,KLIM
1144                L = IND(J)
1145                IND(J) = IND(I)
1146                IND(I) = L
1147                I = I+1
1148                J = J-1
1149   250       CONTINUE
1150             I = JSV
1151          GO TO 220
1152 C
1153   300 CONTINUE
1154 C
1155       IF (NVECT .LE. 0) RETURN
1156       IF (NV .LT. N) GO TO 830
1157 C
1158 C        FIND EIGENVECTORS OF TRI-DIAGONAL MATRIX VIA INVERSE ITERATION
1159 C
1160       IERR = LMSGFL
1161       CALL EINVIT(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,IND,
1162      +            VECT,IERR,B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
1163       IF (IERR .NE. 0) GO TO 840
1164 C
1165 C        FIND EIGENVECTORS OF SYMMETRIC MATRIX VIA BACK TRANSFORMATION
1166 C
1167   400 CONTINUE
1168       CALL ETRBK3(NV,N,LENA,A,NVECT,VECT)
1169       RETURN
1170 C
1171 C        ERROR MESSAGE SECTION
1172 C
1173   800 IF (LMSGFL .LT. 0) RETURN
1174       IF (MASWRK) WRITE(LMSGFL,906)
1175       GO TO 890
1176 C
1177   810 IF (LMSGFL .LT. 0) RETURN
1178       IF (MASWRK) WRITE(LMSGFL,901)
1179       GO TO 890
1180 C
1181   820 IF (LMSGFL .LT. 0) RETURN
1182       IF (MASWRK) WRITE(LMSGFL,902) IERR
1183       GO TO 890
1184 C
1185   830 IF (LMSGFL .LT. 0) RETURN
1186       IF (MASWRK) WRITE(LMSGFL,903)
1187       GO TO 890
1188 C
1189   840 CONTINUE
1190       IF ((LMSGFL .GT. 0).AND.MASWRK) WRITE(LMSGFL,904) -IERR
1191       GO TO 400
1192 C
1193   850 IERR=-1
1194       IF (LMSGFL .LT. 0) RETURN
1195       IF (MASWRK) WRITE(LMSGFL,905)
1196       GO TO 890
1197 C
1198   890 CONTINUE
1199       IF (MASWRK) WRITE(LMSGFL,900) N,NVECT,LENA,NV,IORDER,IERR
1200       RETURN
1201       END
1202 C*MODULE EIGEN   *DECK FREDA
1203       SUBROUTINE FREDA(L,D,A,E)
1204 C
1205       DOUBLE PRECISION A(*)
1206       DOUBLE PRECISION D(L)
1207       DOUBLE PRECISION E(L)
1208       DOUBLE PRECISION F
1209       DOUBLE PRECISION G
1210 C
1211       JK = 1
1212 C
1213 C     .......... FORM REDUCED A ..........
1214 C
1215       DO 280 J = 1, L
1216          F = D(J)
1217          G = E(J)
1218 C
1219          DO 260 K = 1, J
1220             A(JK) = A(JK) - F * E(K) - G * D(K)
1221             JK = JK + 1
1222   260    CONTINUE
1223 C
1224   280 CONTINUE
1225       RETURN
1226       END
1227 C*MODULE EIGEN   *DECK GIVEIS
1228       SUBROUTINE GIVEIS(N,NVECT,NV,A,B,INDB,ROOT,VECT,IERR)
1229       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1230       DIMENSION A(*),B(N,8),INDB(N),ROOT(N),VECT(NV,NVECT)
1231 C
1232 C     EISPACK-BASED SUBSTITUTE FOR QCPE ROUTINE GIVENS.
1233 C     FINDS ALL EIGENVALUES AND SOME EIGENVECTORS OF A REAL SYMMETRIC
1234 C     MATRIX.   AUTHOR.. C. MOLER AND D. SPANGLER, N.R.C.C., 4/1/79.
1235 C
1236 C     INPUT..
1237 C     N     = ORDER OF MATRIX .
1238 C     NVECT = NUMBER OF VECTORS DESIRED.  0 .LE. NVECT .LE. N .
1239 C     NV    = LEADING DIMENSION OF VECT .
1240 C     A     = INPUT MATRIX, COLUMNS OF THE UPPER TRIANGLE PACKED INTO
1241 C             LINEAR ARRAY OF DIMENSION N*(N+1)/2 .
1242 C     B     = SCRATCH ARRAY, 8*N ELEMENTS (NOTE THIS IS MORE THAN
1243 C             PREVIOUS VERSIONS OF GIVENS.)
1244 C    IND    = INDEX ARRAY OF N ELEMENTS
1245 C
1246 C     OUTPUT..
1247 C     A       DESTROYED .
1248 C     ROOT  = ALL EIGENVALUES, ROOT(1) .LE. ... .LE. ROOT(N) .
1249 C             (FOR OTHER ORDERINGS, SEE BELOW.)
1250 C     VECT  = EIGENVECTORS FOR ROOT(1),..., ROOT(NVECT) .
1251 C     IERR  = 0 IF NO ERROR DETECTED,
1252 C           = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
1253 C           = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
1254 C             (FAILURES SHOULD BE VERY RARE.  CONTACT MOLER.)
1255 C
1256 C     CALLS MODIFIED EISPACK ROUTINES TRED3B, IMTQLV, TINVTB, AND
1257 C     TRBK3B.  THE ROUTINES TRED3B, TINVTB, AND TRBK3B.
1258 C     THE ORIGINAL EISPACK ROUTINES TRED3, TINVIT, AND TRBAK3
1259 C     WERE MODIFIED BY THE INTRODUCTION OF TWO ROUTINES FROM THE
1260 C     BLAS LIBRARY - DDOT AND DAXPY.
1261 C
1262 C         IF TINVIT FAILS TO CONVERGE, TQL2 IS CALLED
1263 C
1264 C         SEE EISPACK USERS GUIDE, B. T. SMITH ET AL, SPRINGER-VERLAG
1265 C     LECTURE NOTES IN COMPUTER SCIENCE, VOL. 6, 2-ND EDITION, 1976 .
1266 C     NOTE THAT IMTQLV AND TINVTB HAVE INTERNAL MACHINE
1267 C     DEPENDENT CONSTANTS.
1268 C
1269       DATA ONE, ZERO /1.0D+00, 0.0D+00/
1270       CALL TRED3B(N,(N*N+N)/2,A,B(1,1),B(1,2),B(1,3))
1271       CALL IMTQLV(N,B(1,1),B(1,2),B(1,3),ROOT,INDB,IERR,B(1,4))
1272       IF (IERR .NE. 0) RETURN
1273 C
1274 C     TO REORDER ROOTS...
1275 C     K = N/2
1276 C     B(1,3) = 2.0D+00
1277 C     DO 50 I = 1, K
1278 C        J = N+1-I
1279 C        T = ROOT(I)
1280 C        ROOT(I) = ROOT(J)
1281 C        ROOT(J) = T
1282 C 50  CONTINUE
1283 C
1284       IF (NVECT .LE. 0) RETURN
1285       CALL TINVTB(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,INDB,VECT,IERR,
1286      +     B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
1287       IF (IERR .EQ. 0) GO TO 160
1288 C
1289 C      IF INVERSE ITERATION GIVES AN ERROR IN DETERMINING THE
1290 C      EIGENVECTORS, TRY THE QL ALGORITHM IF ALL THE EIGENVECTORS
1291 C      ARE DESIRED.
1292 C
1293       IF (NVECT .NE. N) RETURN
1294       DO 120 I = 1, NVECT
1295       DO 100 J = 1, N
1296       VECT(I,J) = ZERO
1297   100 CONTINUE
1298       VECT(I,I) = ONE
1299   120 CONTINUE
1300       CALL TQL2 (NV,N,B(1,1),B(1,2),VECT,IERR)
1301       DO 140 I = 1, NVECT
1302       ROOT(I) = B(I,1)
1303   140 CONTINUE
1304       IF (IERR .NE. 0) RETURN
1305   160 CALL TRBK3B(NV,N,(N*N+N)/2,A,NVECT,VECT)
1306       RETURN
1307       END
1308 C*MODULE EIGEN   *DECK GLDIAG
1309       SUBROUTINE GLDIAG(LDVECT,NVECT,N,H,WRK,EIG,VECTOR,IERR,IWRK)
1310 C
1311       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1312 C
1313       LOGICAL GOPARR,DSKWRK,MASWRK
1314 C
1315       DIMENSION H(*),WRK(N,8),EIG(N),VECTOR(LDVECT,NVECT),IWRK(N)
1316 C
1317       COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
1318       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
1319       COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1320 C
1321 C     ----- GENERAL ROUTINE TO DIAGONALIZE A SYMMETRIC MATRIX -----
1322 C     IF KDIAG = 0, USE A ROUTINE FROM THE VECTOR LIBRARY,
1323 C                   IF AVAILABLE (SEE THE SUBROUTINE 'GLDIAG'
1324 C                   IN VECTOR.SRC), OR EVVRSP OTHERWISE
1325 C              = 1, USE EVVRSP
1326 C              = 2, USE GIVEIS
1327 C              = 3, USE JACOBI
1328 C
1329 C           N      = DIMENSION (ORDER) OF MATRIX TO BE SOLVED
1330 C           LDVECT = LEADING DIMENSION OF VECTOR
1331 C           NVECT  = NUMBER OF VECTORS DESIRED
1332 C           H      = MATRIX TO BE DIAGONALIZED
1333 C           WRK    = N*8 W.P. REAL WORDS OF SCRATCH SPACE
1334 C           EIG    = EIGENVECTORS (OUTPUT)
1335 C           VECTOR = EIGENVECTORS (OUTPUT)
1336 C           IERR   = ERROR FLAG (OUTPUT)
1337 C           IWRK   = N INTEGER WORDS OF SCRATCH SPACE
1338 C
1339       IERR = 0
1340 C
1341 C         ----- USE STEVE ELBERT'S ROUTINE -----
1342 C
1343       IF(KDIAG.LE.1  .OR.  KDIAG.GT.3) THEN
1344          LENH = (N*N+N)/2
1345          KORDER =0
1346          CALL EVVRSP(IW,N,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR
1347      *              ,KORDER,IERR)
1348       END IF
1349 C
1350 C         ----- USE MODIFIED EISPAK ROUTINE -----
1351 C
1352       IF(KDIAG.EQ.2)
1353      *   CALL GIVEIS(N,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
1354 C
1355 C         ----- USE JACOBI ROTATION ROUTINE -----
1356 C
1357       IF(KDIAG.EQ.3) THEN
1358          IF(NVECT.EQ.N) THEN
1359             CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,N)
1360          ELSE
1361             IF (MASWRK) WRITE(IW,9000) N,NVECT,LDVECT
1362             CALL ABRT
1363          END IF
1364       END IF
1365       RETURN
1366 C
1367  9000 FORMAT(1X,'IN -GLDIAG-, N,NVECT,LDVECT=',3I8/
1368      *       1X,'THE JACOBI CODE CANNOT COPE WITH N.NE.NVECT!'/
1369      *       1X,'SO THIS RUN DOES NOT PERMIT KDIAG=3.')
1370       END
1371 C*MODULE EIGEN   *DECK IMTQLV
1372       SUBROUTINE IMTQLV(N,D,E,E2,W,IND,IERR,RV1)
1373       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1374       INTEGER TAG
1375       DOUBLE PRECISION MACHEP
1376       DIMENSION D(N),E(N),E2(N),W(N),RV1(N),IND(N)
1377 C
1378 C     THIS ROUTINE IS A VARIANT OF  IMTQL1  WHICH IS A TRANSLATION OF
1379 C     ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND
1380 C     WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
1381 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
1382 C
1383 C     THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
1384 C     MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM
1385 C     THEIR CORRESPONDING SUBMATRIX INDICES.
1386 C
1387 C     ON INPUT-
1388 C
1389 C        N IS THE ORDER OF THE MATRIX,
1390 C
1391 C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1392 C
1393 C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1394 C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
1395 C
1396 C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
1397 C          E2(1) IS ARBITRARY.
1398 C
1399 C     ON OUTPUT-
1400 C
1401 C        D AND E ARE UNALTERED,
1402 C
1403 C        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
1404 C          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
1405 C          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
1406 C          E2(1) IS ALSO SET TO ZERO,
1407 C
1408 C        W CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
1409 C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
1410 C          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
1411 C          THE SMALLEST EIGENVALUES,
1412 C
1413 C        IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE
1414 C          CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES
1415 C          BELONGING TO THE FIRST SUBMATRIX FROM THE TOP,
1416 C          2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.,
1417 C
1418 C        IERR IS SET TO
1419 C          ZERO       FOR NORMAL RETURN,
1420 C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
1421 C                     DETERMINED AFTER 30 ITERATIONS,
1422 C
1423 C        RV1 IS A TEMPORARY STORAGE ARRAY.
1424 C
1425 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1426 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1427 C
1428 C     ------------------------------------------------------------------
1429 C
1430 C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1431 C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1432 C
1433 C                **********
1434       MACHEP = 2.0D+00**(-50)
1435 C
1436       IERR = 0
1437       K = 0
1438       TAG = 0
1439 C
1440       DO 100 I = 1, N
1441       W(I) = D(I)
1442       IF (I .NE. 1) RV1(I-1) = E(I)
1443   100 CONTINUE
1444 C
1445       E2(1) = 0.0D+00
1446       RV1(N) = 0.0D+00
1447 C
1448       DO 360 L = 1, N
1449       J = 0
1450 C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
1451   120 DO 140 M = L, N
1452       IF (M .EQ. N) GO TO 160
1453       IF (ABS(RV1(M)) .LE. MACHEP * (ABS(W(M)) + ABS(W(M+1)))) GO TO
1454      +     160
1455 C     ********** GUARD AGAINST UNDERFLOWED ELEMENT OF E2 **********
1456       IF (E2(M+1) .EQ. 0.0D+00) GO TO 180
1457   140 CONTINUE
1458 C
1459   160 IF (M .LE. K) GO TO 200
1460       IF (M .NE. N) E2(M+1) = 0.0D+00
1461   180 K = M
1462       TAG = TAG + 1
1463   200 P = W(L)
1464       IF (M .EQ. L) GO TO 280
1465       IF (J .EQ. 30) GO TO 380
1466       J = J + 1
1467 C     ********** FORM SHIFT **********
1468       G = (W(L+1) - P) / (2.0D+00 * RV1(L))
1469       R = SQRT(G*G+1.0D+00)
1470       G = W(M) - P + RV1(L) / (G + SIGN(R,G))
1471       S = 1.0D+00
1472       C = 1.0D+00
1473       P = 0.0D+00
1474       MML = M - L
1475 C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
1476       DO 260 II = 1, MML
1477       I = M - II
1478       F = S * RV1(I)
1479       B = C * RV1(I)
1480       IF (ABS(F) .LT. ABS(G)) GO TO 220
1481       C = G / F
1482       R = SQRT(C*C+1.0D+00)
1483       RV1(I+1) = F * R
1484       S = 1.0D+00 / R
1485       C = C * S
1486       GO TO 240
1487   220 S = F / G
1488       R = SQRT(S*S+1.0D+00)
1489       RV1(I+1) = G * R
1490       C = 1.0D+00 / R
1491       S = S * C
1492   240 G = W(I+1) - P
1493       R = (W(I) - G) * S + 2.0D+00 * C * B
1494       P = S * R
1495       W(I+1) = G + P
1496       G = C * R - B
1497   260 CONTINUE
1498 C
1499       W(L) = W(L) - P
1500       RV1(L) = G
1501       RV1(M) = 0.0D+00
1502       GO TO 120
1503 C     ********** ORDER EIGENVALUES **********
1504   280 IF (L .EQ. 1) GO TO 320
1505 C     ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
1506       DO 300 II = 2, L
1507       I = L + 2 - II
1508       IF (P .GE. W(I-1)) GO TO 340
1509       W(I) = W(I-1)
1510       IND(I) = IND(I-1)
1511   300 CONTINUE
1512 C
1513   320 I = 1
1514   340 W(I) = P
1515       IND(I) = TAG
1516   360 CONTINUE
1517 C
1518       GO TO 400
1519 C     ********** SET ERROR -- NO CONVERGENCE TO AN
1520 C                EIGENVALUE AFTER 30 ITERATIONS **********
1521   380 IERR = L
1522   400 RETURN
1523 C     ********** LAST CARD OF IMTQLV **********
1524       END
1525 C*MODULE EIGEN   *DECK JACDG
1526       SUBROUTINE JACDG(A,VEC,EIG,JBIG,BIG,LDVEC,N)
1527 C
1528       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1529 C
1530       DIMENSION A(*),VEC(LDVEC,N),EIG(N),JBIG(N),BIG(N)
1531 C
1532       PARAMETER (ONE=1.0D+00)
1533 C
1534 C     ----- JACOBI DIAGONALIZATION OF SYMMETRIC MATRIX -----
1535 C     SYMMETRIC MATRIX -A- OF DIMENSION -N- IS DESTROYED ON EXIT.
1536 C     ALL EIGENVECTORS ARE FOUND, SO -VEC- MUST BE SQUARE,
1537 C     UNLESS SOMEONE TAKES THE TROUBLE TO LOOK AT -NMAX- BELOW.
1538 C     -BIG- AND -JBIG- ARE SCRATCH WORK ARRAYS.
1539 C
1540       CALL VCLR(VEC,1,LDVEC*N)
1541       DO 20 I = 1,N
1542         VEC(I,I) = ONE
1543    20 CONTINUE
1544 C
1545       NB1 = N
1546       NB2 = (NB1*NB1+NB1)/2
1547       NMIN = 1
1548       NMAX = NB1
1549 C
1550       CALL JACDIA(A,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1551 C
1552       DO 30 I=1,N
1553         EIG(I) = A((I*I+I)/2)
1554    30 CONTINUE
1555 C
1556       CALL JACORD(VEC,EIG,NB1,LDVEC)
1557       RETURN
1558       END
1559 C*MODULE EIGEN   *DECK JACDIA
1560       SUBROUTINE JACDIA(F,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1561       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1562       LOGICAL GOPARR,DSKWRK,MASWRK
1563       DIMENSION F(NB2),VEC(LDVEC,NB1),BIG(NB1),JBIG(NB1)
1564 C
1565       COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1566 C
1567       PARAMETER (ROOT2=0.707106781186548D+00 )
1568       PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, D1050=1.05D+00,
1569      *           D1500=1.5D+00, D3875=3.875D+00,
1570      *           D0500=0.5D+00, D1375=1.375D+00, D0250=0.25D+00 )
1571       PARAMETER (C2=1.0D-12, C3=4.0D-16,
1572      *           C4=2.0D-16, C5=8.0D-09, C6=3.0D-06 )
1573 C
1574 C      F IS THE MATRIX TO BE DIAGONALIZED, F IS STORED TRIANGULAR
1575 C      VEC IS THE ARRAY OF EIGENVECTORS, DIMENSION NB1*NB1
1576 C      BIG AND JBIG ARE TEMPORARY SCRATCH AREAS OF DIMENSION NB1
1577 C      THE ROTATIONS AMONG THE FIRST NMIN BASIS FUNCTIONS ARE NOT
1578 C      ACCOUNTED FOR.
1579 C      THE ROTATIONS AMONG THE LAST NB1-NMAX BASIS FUNCTIONS ARE NOT
1580 C      ACCOUNTED FOR.
1581 C
1582       IEAA=0
1583       IEAB=0
1584       TT=ZERO
1585       EPS = 64.0D+00*EPSLON(ONE)
1586 C
1587 C      LOOP OVER COLUMNS (K) OF TRIANGULAR MATRIX TO DETERMINE
1588 C      LARGEST OFF-DIAGONAL ELEMENTS IN ROW(I).
1589 C
1590       DO 20 I=1,NB1
1591          BIG(I)=ZERO
1592          JBIG(I)=0
1593          IF(I.LT.NMIN  .OR.  I.EQ.1) GO TO 20
1594          II = (I*I-I)/2
1595          J=MIN(I-1,NMAX)
1596          DO 10 K=1,J
1597             IF(ABS(BIG(I)).GE.ABS(F(II+K))) GO TO 10
1598             BIG(I)=F(II+K)
1599             JBIG(I)=K
1600    10    CONTINUE
1601    20 CONTINUE
1602 C
1603 C     ----- 2X2 JACOBI ITERATIONS BEGIN HERE -----
1604 C
1605       MAXIT=MAX(NB2*20,500)
1606       ITER=0
1607    30 CONTINUE
1608       ITER=ITER+1
1609 C
1610 C      FIND SMALLEST DIAGONAL ELEMENT
1611 C
1612       SD=D1050
1613       JJ=0
1614       DO 40 J=1,NB1
1615          JJ=JJ+J
1616          SD= MIN(SD,ABS(F(JJ)))
1617    40 CONTINUE
1618       TEST = MAX(EPS, C2*MAX(SD,C6))
1619 C
1620 C      FIND LARGEST OFF-DIAGONAL ELEMENT
1621 C
1622       T=ZERO
1623       I1=MAX(2,NMIN)
1624       IB = I1
1625       DO 50 I=I1,NB1
1626          IF(T.GE.ABS(BIG(I))) GO TO 50
1627          T= ABS(BIG(I))
1628          IB=I
1629    50 CONTINUE
1630 C
1631 C      TEST FOR CONVERGENCE, THEN DETERMINE ROTATION.
1632 C
1633       IF(T.LT.TEST) RETURN
1634 C                   ******
1635 C
1636       IF(ITER.GT.MAXIT) THEN
1637          IF (MASWRK) THEN
1638             WRITE(6,*) 'JACOBI DIAGONALIZATION FAILS, DIMENSION=',NB1
1639             WRITE(6,9020) ITER,T,TEST,SD
1640          ENDIF
1641          CALL ABRT
1642          STOP
1643       END IF
1644 C
1645       IA=JBIG(IB)
1646       IAA=IA*(IA-1)/2
1647       IBB=IB*(IB-1)/2
1648       DIF=F(IAA+IA)-F(IBB+IB)
1649       IF(ABS(DIF).GT.C3*T) GO TO 70
1650       SX=ROOT2
1651       CX=ROOT2
1652       GO TO 110
1653    70 T2X2=BIG(IB)/DIF
1654       T2X25=T2X2*T2X2
1655       IF(T2X25 . GT . C4) GO TO 80
1656       CX=ONE
1657       SX=T2X2
1658       GO TO 110
1659    80 IF(T2X25 . GT . C5) GO TO 90
1660       SX=T2X2*(ONE-D1500*T2X25)
1661       CX=ONE-D0500*T2X25
1662       GO TO 110
1663    90 IF(T2X25 . GT . C6) GO TO 100
1664       CX=ONE+T2X25*(T2X25*D1375 - D0500)
1665       SX= T2X2*(ONE + T2X25*(T2X25*D3875 - D1500))
1666       GO TO 110
1667   100 T=D0250  / SQRT(D0250   + T2X25)
1668       CX= SQRT(D0500   + T)
1669       SX= SIGN( SQRT(D0500   - T),T2X2)
1670   110 IEAR=IAA+1
1671       IEBR=IBB+1
1672 C
1673       DO 230 IR=1,NB1
1674          T=F(IEAR)*SX
1675          F(IEAR)=F(IEAR)*CX+F(IEBR)*SX
1676          F(IEBR)=T-F(IEBR)*CX
1677          IF(IR-IA) 220,120,130
1678   120    TT=F(IEBR)
1679          IEAA=IEAR
1680          IEAB=IEBR
1681          F(IEBR)=BIG(IB)
1682          IEAR=IEAR+IR-1
1683          IF(JBIG(IR)) 200,220,200
1684   130    T=F(IEAR)
1685          IT=IA
1686          IEAR=IEAR+IR-1
1687          IF(IR-IB) 180,150,160
1688   150    F(IEAA)=F(IEAA)*CX+F(IEAB)*SX
1689          F(IEAB)=TT*CX+F(IEBR)*SX
1690          F(IEBR)=TT*SX-F(IEBR)*CX
1691          IEBR=IEBR+IR-1
1692          GO TO 200
1693   160    IF(  ABS(T) . GE .  ABS(F(IEBR))) GO TO 170
1694          IF(IB.GT.NMAX) GO TO 170
1695          T=F(IEBR)
1696          IT=IB
1697   170    IEBR=IEBR+IR-1
1698   180    IF(  ABS(T) . LT .  ABS(BIG(IR))) GO TO 190
1699          BIG(IR) = T
1700          JBIG(IR) = IT
1701          GO TO 220
1702   190    IF(IA . NE . JBIG(IR) . AND . IB . NE . JBIG(IR))  GO TO 220
1703   200    KQ=IEAR-IR-IA+1
1704          BIG(IR)=ZERO
1705          IR1=MIN(IR-1,NMAX)
1706          DO 210 I=1,IR1
1707             K=KQ+I
1708             IF(ABS(BIG(IR)) . GE . ABS(F(K)))  GO TO 210
1709             BIG(IR) = F(K)
1710             JBIG(IR)=I
1711   210    CONTINUE
1712   220    IEAR=IEAR+1
1713   230    IEBR=IEBR+1
1714 C
1715       DO 240 I=1,NB1
1716          T1=VEC(I,IA)*CX + VEC(I,IB)*SX
1717          T2=VEC(I,IA)*SX - VEC(I,IB)*CX
1718          VEC(I,IA)=T1
1719          VEC(I,IB)=T2
1720   240 CONTINUE
1721       GO TO 30
1722 C
1723  9020 FORMAT(1X,'ITER=',I6,' T,TEST,SD=',1P,3E20.10)
1724       END
1725 C*MODULE EIGEN   *DECK JACORD
1726       SUBROUTINE JACORD(VEC,EIG,N,LDVEC)
1727       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1728       DIMENSION VEC(LDVEC,N),EIG(N)
1729 C
1730 C     ---- SORT EIGENDATA INTO ASCENDING ORDER -----
1731 C
1732       DO 290 I = 1, N
1733          JJ = I
1734          DO 270 J = I, N
1735             IF (EIG(J) .LT. EIG(JJ)) JJ = J
1736   270    CONTINUE
1737          IF (JJ .EQ. I) GO TO 290
1738          T = EIG(JJ)
1739          EIG(JJ) = EIG(I)
1740          EIG(I) = T
1741          DO 280 J = 1, N
1742             T = VEC(J,JJ)
1743             VEC(J,JJ) = VEC(J,I)
1744             VEC(J,I) = T
1745   280    CONTINUE
1746   290 CONTINUE
1747       RETURN
1748       END
1749 C*MODULE EIGEN   *DECK TINVTB
1750       SUBROUTINE TINVTB(NM,N,D,E,E2,M,W,IND,Z,
1751      *                  IERR,RV1,RV2,RV3,RV4,RV6)
1752       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1753       DIMENSION D(N),E(N),E2(N),W(M),Z(NM,M),
1754      *          RV1(N),RV2(N),RV3(N),RV4(N),RV6(N),IND(M)
1755       DOUBLE PRECISION MACHEP,NORM
1756       INTEGER P,Q,R,S,TAG,GROUP
1757 C     ------------------------------------------------------------------
1758 C
1759 C     THIS ROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH-
1760 C     NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
1761 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
1762 C
1763 C     THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
1764 C     SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
1765 C     USING INVERSE ITERATION.
1766 C
1767 C     ON INPUT-
1768 C
1769 C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1770 C          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
1771 C          DIMENSION STATEMENT,
1772 C
1773 C        N IS THE ORDER OF THE MATRIX,
1774 C
1775 C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1776 C
1777 C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1778 C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
1779 C
1780 C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E,
1781 C          WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
1782 C          E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
1783 C          THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM
1784 C          OF THE MAGNITUDES OF D(I) AND D(I-1).  E2(1) MUST CONTAIN
1785 C          0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0
1786 C          IF THE EIGENVALUES ARE IN DESCENDING ORDER.  IF  BISECT,
1787 C          TRIDIB, OR  IMTQLV  HAS BEEN USED TO FIND THE EIGENVALUES,
1788 C          THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE,
1789 C
1790 C        M IS THE NUMBER OF SPECIFIED EIGENVALUES,
1791 C
1792 C        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER,
1793 C
1794 C        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
1795 C          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
1796 C          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
1797 C          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.
1798 C
1799 C     ON OUTPUT-
1800 C
1801 C        ALL INPUT ARRAYS ARE UNALTERED,
1802 C
1803 C        Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
1804 C          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO,
1805 C
1806 C        IERR IS SET TO
1807 C          ZERO       FOR NORMAL RETURN,
1808 C          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
1809 C                     EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS,
1810 C
1811 C        RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
1812 C
1813 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1814 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1815 C
1816 C     ------------------------------------------------------------------
1817 C
1818 C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1819 C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1820 C
1821 C                **********
1822       MACHEP = 2.0D+00**(-50)
1823 C
1824       IERR = 0
1825       IF (M .EQ. 0) GO TO 680
1826       TAG = 0
1827       ORDER = 1.0D+00 - E2(1)
1828       XU = 0.0D+00
1829       UK = 0.0D+00
1830       X0 = 0.0D+00
1831       U  = 0.0D+00
1832       EPS2 = 0.0D+00
1833       EPS3 = 0.0D+00
1834       EPS4 = 0.0D+00
1835       GROUP = 0
1836       Q = 0
1837 C     ********** ESTABLISH AND PROCESS NEXT SUBMATRIX **********
1838   100 P = Q + 1
1839       IP = P + 1
1840 C
1841       DO 120 Q = P, N
1842       IF (Q .EQ. N) GO TO 140
1843       IF (E2(Q+1) .EQ. 0.0D+00) GO TO 140
1844   120 CONTINUE
1845 C     ********** FIND VECTORS BY INVERSE ITERATION **********
1846   140 TAG = TAG + 1
1847       IQMP = Q - P + 1
1848       S = 0
1849 C
1850       DO 660 R = 1, M
1851       IF (IND(R) .NE. TAG) GO TO 660
1852       ITS = 1
1853       X1 = W(R)
1854       IF (S .NE. 0) GO TO 220
1855 C     ********** CHECK FOR ISOLATED ROOT **********
1856       XU = 1.0D+00
1857       IF (P .NE. Q) GO TO 160
1858       RV6(P) = 1.0D+00
1859       GO TO 600
1860   160 NORM = ABS(D(P))
1861 C
1862       DO 180 I = IP, Q
1863   180 NORM = NORM + ABS(D(I)) + ABS(E(I))
1864 C     ********** EPS2 IS THE CRITERION FOR GROUPING,
1865 C                EPS3 REPLACES ZERO PIVOTS AND EQUAL
1866 C                ROOTS ARE MODIFIED BY EPS3,
1867 C                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW **********
1868       EPS2 = 1.0D-03 * NORM
1869       EPS3 = MACHEP * NORM
1870       UK = IQMP
1871       EPS4 = UK * EPS3
1872       UK = EPS4 / SQRT(UK)
1873       S = P
1874   200 GROUP = 0
1875       GO TO 240
1876 C     ********** LOOK FOR CLOSE OR COINCIDENT ROOTS **********
1877   220 IF (ABS(X1-X0) .GE. EPS2) GO TO 200
1878       GROUP = GROUP + 1
1879       IF (ORDER * (X1 - X0) .LE. 0.0D+00) X1 = X0 + ORDER * EPS3
1880 C     ********** ELIMINATION WITH INTERCHANGES AND
1881 C                INITIALIZATION OF VECTOR **********
1882   240 V = 0.0D+00
1883 C
1884       DO 300 I = P, Q
1885       RV6(I) = UK
1886       IF (I .EQ. P) GO TO 280
1887       IF (ABS(E(I)) .LT. ABS(U)) GO TO 260
1888 C     ********** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
1889 C                E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY **********
1890       XU = U / E(I)
1891       RV4(I) = XU
1892       RV1(I-1) = E(I)
1893       RV2(I-1) = D(I) - X1
1894       RV3(I-1) = 0.0D+00
1895       IF (I .NE. Q) RV3(I-1) = E(I+1)
1896       U = V - XU * RV2(I-1)
1897       V = -XU * RV3(I-1)
1898       GO TO 300
1899   260 XU = E(I) / U
1900       RV4(I) = XU
1901       RV1(I-1) = U
1902       RV2(I-1) = V
1903       RV3(I-1) = 0.0D+00
1904   280 U = D(I) - X1 - XU * V
1905       IF (I .NE. Q) V = E(I+1)
1906   300 CONTINUE
1907 C
1908       IF (U .EQ. 0.0D+00) U = EPS3
1909       RV1(Q) = U
1910       RV2(Q) = 0.0D+00
1911       RV3(Q) = 0.0D+00
1912 C     ********** BACK SUBSTITUTION
1913 C                FOR I=Q STEP -1 UNTIL P DO -- **********
1914   320 DO 340 II = P, Q
1915       I = P + Q - II
1916       RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
1917       V = U
1918       U = RV6(I)
1919   340 CONTINUE
1920 C     ********** ORTHOGONALIZE WITH RESPECT TO PREVIOUS
1921 C                MEMBERS OF GROUP **********
1922       IF (GROUP .EQ. 0) GO TO 400
1923       J = R
1924 C
1925       DO 380 JJ = 1, GROUP
1926   360 J = J - 1
1927       IF (IND(J) .NE. TAG) GO TO 360
1928       XU = DDOT(IQMP,RV6(P),1,Z(P,J),1)
1929 C
1930       CALL DAXPY(IQMP,-XU,Z(P,J),1,RV6(P),1)
1931 C
1932   380 CONTINUE
1933 C
1934   400 NORM = 0.0D+00
1935 C
1936       DO 420 I = P, Q
1937   420 NORM = NORM + ABS(RV6(I))
1938 C
1939       IF (NORM .GE. 1.0D+00) GO TO 560
1940 C     ********** FORWARD SUBSTITUTION **********
1941       IF (ITS .EQ. 5) GO TO 540
1942       IF (NORM .NE. 0.0D+00) GO TO 440
1943       RV6(S) = EPS4
1944       S = S + 1
1945       IF (S .GT. Q) S = P
1946       GO TO 480
1947   440 XU = EPS4 / NORM
1948 C
1949       DO 460 I = P, Q
1950   460 RV6(I) = RV6(I) * XU
1951 C     ********** ELIMINATION OPERATIONS ON NEXT VECTOR
1952 C                ITERATE **********
1953   480 DO 520 I = IP, Q
1954       U = RV6(I)
1955 C     ********** IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
1956 C                WAS PERFORMED EARLIER IN THE
1957 C                TRIANGULARIZATION PROCESS **********
1958       IF (RV1(I-1) .NE. E(I)) GO TO 500
1959       U = RV6(I-1)
1960       RV6(I-1) = RV6(I)
1961   500 RV6(I) = U - RV4(I) * RV6(I-1)
1962   520 CONTINUE
1963 C
1964       ITS = ITS + 1
1965       GO TO 320
1966 C     ********** SET ERROR -- NON-CONVERGED EIGENVECTOR **********
1967   540 IERR = -R
1968       XU = 0.0D+00
1969       GO TO 600
1970 C     ********** NORMALIZE SO THAT SUM OF SQUARES IS
1971 C                1 AND EXPAND TO FULL ORDER **********
1972   560 U = 0.0D+00
1973 C
1974       DO 580 I = P, Q
1975       RV6(I) = RV6(I) / NORM
1976   580 U = U + RV6(I)**2
1977 C
1978       XU = 1.0D+00 / SQRT(U)
1979 C
1980   600 DO 620 I = 1, N
1981   620 Z(I,R) = 0.0D+00
1982 C
1983       DO 640 I = P, Q
1984   640 Z(I,R) = RV6(I) * XU
1985 C
1986       X0 = X1
1987   660 CONTINUE
1988 C
1989       IF (Q .LT. N) GO TO 100
1990   680 RETURN
1991 C     ********** LAST CARD OF TINVIT **********
1992       END
1993 C*MODULE EIGEN   *DECK TQL2
1994 C
1995 C     ------------------------------------------------------------------
1996 C
1997       SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
1998       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1999       DOUBLE PRECISION MACHEP
2000       DIMENSION D(N),E(N),Z(NM,N)
2001 C
2002 C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
2003 C     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
2004 C     WILKINSON.
2005 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
2006 C
2007 C     THIS ROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
2008 C     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
2009 C     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
2010 C     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
2011 C     FULL MATRIX TO TRIDIAGONAL FORM.
2012 C
2013 C     ON INPUT-
2014 C
2015 C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2016 C          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2017 C          DIMENSION STATEMENT,
2018 C
2019 C        N IS THE ORDER OF THE MATRIX,
2020 C
2021 C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2022 C
2023 C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2024 C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
2025 C
2026 C        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
2027 C          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
2028 C          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
2029 C          THE IDENTITY MATRIX.
2030 C
2031 C      ON OUTPUT-
2032 C
2033 C        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
2034 C          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
2035 C          UNORDERED FOR INDICES 1,2,...,IERR-1,
2036 C
2037 C        E HAS BEEN DESTROYED,
2038 C
2039 C        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
2040 C          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
2041 C          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
2042 C          EIGENVALUES,
2043 C
2044 C        IERR IS SET TO
2045 C          ZERO       FOR NORMAL RETURN,
2046 C          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
2047 C                     DETERMINED AFTER 30 ITERATIONS.
2048 C
2049 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2050 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2051 C
2052 C     ------------------------------------------------------------------
2053 C
2054 C     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2055 C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2056 C
2057 C                **********
2058       MACHEP = 2.0D+00**(-50)
2059 C
2060       IERR = 0
2061       IF (N .EQ. 1) GO TO 400
2062 C
2063       DO 100 I = 2, N
2064   100 E(I-1) = E(I)
2065 C
2066       F = 0.0D+00
2067       B = 0.0D+00
2068       E(N) = 0.0D+00
2069 C
2070       DO 300 L = 1, N
2071       J = 0
2072       H = MACHEP * (ABS(D(L)) + ABS(E(L)))
2073       IF (B .LT. H) B = H
2074 C     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
2075       DO 120 M = L, N
2076       IF (ABS(E(M)) .LE. B) GO TO 140
2077 C     ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
2078 C                THROUGH THE BOTTOM OF THE LOOP **********
2079   120 CONTINUE
2080 C
2081   140 IF (M .EQ. L) GO TO 280
2082   160 IF (J .EQ. 30) GO TO 380
2083       J = J + 1
2084 C     ********** FORM SHIFT **********
2085       L1 = L + 1
2086       G = D(L)
2087       P = (D(L1) - G) / (2.0D+00 * E(L))
2088       R = SQRT(P*P+1.0D+00)
2089       D(L) = E(L) / (P + SIGN(R,P))
2090       H = G - D(L)
2091 C
2092       DO 180 I = L1, N
2093   180 D(I) = D(I) - H
2094 C
2095       F = F + H
2096 C     ********** QL TRANSFORMATION **********
2097       P = D(M)
2098       C = 1.0D+00
2099       S = 0.0D+00
2100       MML = M - L
2101 C     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
2102       DO 260 II = 1, MML
2103       I = M - II
2104       G = C * E(I)
2105       H = C * P
2106       IF (ABS(P) .LT. ABS(E(I))) GO TO 200
2107       C = E(I) / P
2108       R = SQRT(C*C+1.0D+00)
2109       E(I+1) = S * P * R
2110       S = C / R
2111       C = 1.0D+00 / R
2112       GO TO 220
2113   200 C = P / E(I)
2114       R = SQRT(C*C+1.0D+00)
2115       E(I+1) = S * E(I) * R
2116       S = 1.0D+00 / R
2117       C = C * S
2118   220 P = C * D(I) - S * G
2119       D(I+1) = H + S * (C * G + S * D(I))
2120 C     ********** FORM VECTOR **********
2121       CALL DROT(N,Z(1,I+1),1,Z(1,I),1,C,S)
2122 C
2123   260 CONTINUE
2124 C
2125       E(L) = S * P
2126       D(L) = C * P
2127       IF (ABS(E(L)) .GT. B) GO TO 160
2128   280 D(L) = D(L) + F
2129   300 CONTINUE
2130 C     ********** ORDER EIGENVALUES AND EIGENVECTORS **********
2131       DO 360 II = 2, N
2132       I = II - 1
2133       K = I
2134       P = D(I)
2135 C
2136       DO 320 J = II, N
2137       IF (D(J) .GE. P) GO TO 320
2138       K = J
2139       P = D(J)
2140   320 CONTINUE
2141 C
2142       IF (K .EQ. I) GO TO 360
2143       D(K) = D(I)
2144       D(I) = P
2145 C
2146       CALL DSWAP(N,Z(1,I),1,Z(1,K),1)
2147 C
2148   360 CONTINUE
2149 C
2150       GO TO 400
2151 C     ********** SET ERROR -- NO CONVERGENCE TO AN
2152 C                EIGENVALUE AFTER 30 ITERATIONS **********
2153   380 IERR = L
2154   400 RETURN
2155 C     ********** LAST CARD OF TQL2 **********
2156       END
2157 C*MODULE EIGEN   *DECK TRBK3B
2158 C
2159 C     ------------------------------------------------------------------
2160 C
2161       SUBROUTINE TRBK3B(NM,N,NV,A,M,Z)
2162       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2163       DIMENSION A(NV),Z(NM,M)
2164 C
2165 C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
2166 C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
2167 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
2168 C
2169 C     THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
2170 C     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
2171 C     SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  TRED3B.
2172 C
2173 C     ON INPUT-
2174 C
2175 C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2176 C          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2177 C          DIMENSION STATEMENT,
2178 C
2179 C        N IS THE ORDER OF THE MATRIX,
2180 C
2181 C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2182 C          AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2183 C
2184 C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS
2185 C          USED IN THE REDUCTION BY  TRED3B IN ITS FIRST
2186 C          N*(N+1)/2 POSITIONS,
2187 C
2188 C        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
2189 C
2190 C        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
2191 C          IN ITS FIRST M COLUMNS.
2192 C
2193 C     ON OUTPUT-
2194 C
2195 C        Z CONTAINS THE TRANSFORMED EIGENVECTORS
2196 C          IN ITS FIRST M COLUMNS.
2197 C
2198 C     NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
2199 C
2200 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2201 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2202 C
2203 C     ------------------------------------------------------------------
2204 C
2205       IF (M .EQ. 0) GO TO 140
2206       IF (N .EQ. 1) GO TO 140
2207 C
2208       DO 120 I = 2, N
2209       L = I - 1
2210       IZ = (I * L) / 2
2211       IK = IZ + I
2212       H = A(IK)
2213       IF (H .EQ. 0.0D+00) GO TO 120
2214 C
2215       DO 100 J = 1, M
2216       S = -DDOT(L,A(IZ+1),1,Z(1,J),1)
2217 C
2218 C     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
2219       S = (S / H) / H
2220 C
2221       CALL DAXPY(L,S,A(IZ+1),1,Z(1,J),1)
2222 C
2223   100 CONTINUE
2224 C
2225   120 CONTINUE
2226 C
2227   140 RETURN
2228 C     ********** LAST CARD OF TRBAK3 **********
2229       END
2230 C*MODULE EIGEN   *DECK TRED3B
2231 C
2232 C     ------------------------------------------------------------------
2233 C
2234       SUBROUTINE TRED3B(N,NV,A,D,E,E2)
2235       IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2236       DIMENSION A(NV),D(N),E(N),E2(N)
2237 C
2238 C     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
2239 C     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
2240 C     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
2241 C
2242 C     THIS ROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS
2243 C     A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
2244 C     USING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
2245 C
2246 C     ON INPUT-
2247 C
2248 C        N IS THE ORDER OF THE MATRIX,
2249 C
2250 C        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2251 C          AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2252 C
2253 C        A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
2254 C          INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
2255 C          ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
2256 C
2257 C     ON OUTPUT-
2258 C
2259 C        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
2260 C          TRANSFORMATIONS USED IN THE REDUCTION,
2261 C
2262 C        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
2263 C
2264 C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
2265 C          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO,
2266 C
2267 C        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
2268 C          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
2269 C
2270 C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2271 C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2272 C
2273 C     ------------------------------------------------------------------
2274 C
2275 C     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
2276       DO 300 II = 1, N
2277       I = N + 1 - II
2278       L = I - 1
2279       IZ = (I * L) / 2
2280       H = 0.0D+00
2281       SCALE = 0.0D+00
2282       IF (L .LT. 1) GO TO 120
2283 C     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
2284       DO 100 K = 1, L
2285       IZ = IZ + 1
2286       D(K) = A(IZ)
2287       SCALE = SCALE + ABS(D(K))
2288   100 CONTINUE
2289 C
2290       IF (SCALE .NE. 0.0D+00) GO TO 140
2291   120 E(I) = 0.0D+00
2292       E2(I) = 0.0D+00
2293       GO TO 280
2294 C
2295   140 DO 160 K = 1, L
2296       D(K) = D(K) / SCALE
2297       H = H + D(K) * D(K)
2298   160 CONTINUE
2299 C
2300       E2(I) = SCALE * SCALE * H
2301       F = D(L)
2302       G = -SIGN(SQRT(H),F)
2303       E(I) = SCALE * G
2304       H = H - F * G
2305       D(L) = F - G
2306       A(IZ) = SCALE * D(L)
2307       IF (L .EQ. 1) GO TO 280
2308       F = 0.0D+00
2309 C
2310       JK = 1
2311       DO 220 J = 1, L
2312       JM1 = J - 1
2313       DT = D(J)
2314       G = 0.0D+00
2315 C     ********** FORM ELEMENT OF A*U **********
2316       IF (JM1 .EQ. 0) GO TO 200
2317       DO 180 K = 1, JM1
2318       E(K) = E(K) + DT * A(JK)
2319       G = G + D(K) * A(JK)
2320       JK = JK + 1
2321   180 CONTINUE
2322   200 E(J) = G + A(JK) * DT
2323       JK = JK + 1
2324 C     ********** FORM ELEMENT OF P **********
2325   220 CONTINUE
2326       F = 0.0D+00
2327       DO 240 J = 1, L
2328       E(J) = E(J) / H
2329       F = F + E(J) * D(J)
2330   240 CONTINUE
2331 C
2332       HH = F / (H + H)
2333       JK = 0
2334 C     ********** FORM REDUCED A **********
2335       DO 260 J = 1, L
2336       F = D(J)
2337       G = E(J) - HH * F
2338       E(J) = G
2339 C
2340       DO 260 K = 1, J
2341       JK = JK + 1
2342       A(JK) = A(JK) - F * E(K) - G * D(K)
2343   260 CONTINUE
2344 C
2345   280 D(I) = A(IZ+1)
2346       A(IZ+1) = SCALE * SQRT(H)
2347   300 CONTINUE
2348 C
2349       RETURN
2350 C     ********** LAST CARD OF TRED3 **********
2351       END