source/unres/md_calc.f90

   1       module MD_calc
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use MD_data, only:D_ban,IP
   5       use geometry_data
   6 !      use prng ! prng.f90 or prng_32.f90
   7       implicit none
   8 !
   9 !-----------------------------------------------------------------------------
  10       contains
  11 !-----------------------------------------------------------------------------
  12 ! add.f
  13 !-----------------------------------------------------------------------------
  14       subroutine ABRT
  15       STOP 'IN ABRT'
  16       end subroutine ABRT
  17 !-----------------------------------------------------------------------------
  18 !*MODULE MTHLIB  *DECK VCLR
  19       subroutine VCLR(A,INCA,N)
  20 !
  21 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
  22 !
  23       real(kind=8),DIMENSION(*) :: A
  24 !
  25       real(kind=8),PARAMETER :: ZERO=0.0D+00
  26       integer :: INCA,N
  27       integer :: l,la
  28 !
  29 !     ----- ZERO OUT VECTOR -A-, USING INCREMENT -INCA- -----
  30 !
  31       IF (INCA .NE. 1) GO TO 200
  32       DO 110 L=1,N
  33          A(L) = ZERO
  34   110 CONTINUE
  35       RETURN
  36 !
  37   200 CONTINUE
  38       LA=1-INCA
  39       DO 210 L=1,N
  40          LA=LA+INCA
  41          A(LA) = ZERO
  42   210 CONTINUE
  43       return
  44       end subroutine VCLR
  45 !-----------------------------------------------------------------------------
  46 ! banach.f
  47 !-----------------------------------------------------------------------------
  48       subroutine BANACH(N,NMAX,A,X,osob)
  49 !**********************
  50 !     Banachiewicz
  51 !      implicit real*8 (a-h,o-z)
  52 !      include 'DIMENSIONS'
  53       integer :: N,NMAX
  54       real(kind=8),DIMENSION(NMAX,NMAX) :: A
  55       real(kind=8),DIMENSION(NMAX) :: X
  56 !el      real(kind=8),DIMENSION(6*nres) :: D    !(MAXRES6) maxres6=6*maxres
  57 !el      COMMON /BANII/ D
  58       logical :: osob
  59       real(kind=8) :: xx,aij,aijd
  60       integer :: i,j,k,jjjj
  61
  62 !el      allocate(D_ban(6*nres))
  63
  64       osob=.false.
  65       if (dabs(a(1,1)).lt.1.0d-15) then
  66         osob=.true.
  67         return
  68       endif
  69       D_ban(1)=1./A(1,1)
  70       DO 80 I=2,N
  71       A(I,1)=A(1,I)
  72       DO 81 J=2,I-1
  73       XX=A(J,I)
  74       DO 82 K=1,J-1
  75       XX=XX-A(I,K)*A(J,K)
  76    82 CONTINUE
  77       A(I,J)=XX
  78    81 CONTINUE
  79       XX=A(I,I)
  80       JJJJ=I-1
  81       DO 83 J=1,JJJJ
  82       AIJ=A(I,J)
  83       AIJD=AIJ*D_ban(J)
  84       A(I,J)=AIJD
  85       XX=XX-AIJ*AIJD
  86    83 CONTINUE
  87       if (dabs(xx).lt.1.0d-15) then
  88         osob=.true.
  89         return
  90       endif
  91       D_ban(I)=1./XX
  92    80 CONTINUE
  93 !
  94       CALL BANAII(N,NMAX,A,X)
  95       return
  96       end subroutine BANACH
  97 !-----------------------------------------------------------------------------
  98       subroutine BANAII(N,NMAX,A,X)
  99 !************************
 100 !      implicit real*8 (a-h,o-z)
 101 !      include 'DIMENSIONS'
 102       integer :: N,NMAX
 103       real(kind=8),DIMENSION(NMAX,NMAX) :: A
 104       real(kind=8),DIMENSION(NMAX) :: X
 105 !el      real(kind=8),DIMENSION(6*nres) :: D    !(MAXRES6) maxres6=6*maxres
 106 !el      COMMON /BANII/ D   --->   D_ban
 107       real(kind=8) :: Z
 108       integer :: i,j,jjjj
 109       DO 90 I=1,N
 110       Z=X(I)
 111       JJJJ=I-1
 112       DO 91 J=JJJJ,1,-1
 113       Z=Z-A(I,J)*X(J)
 114    91 CONTINUE
 115       X(I)=Z
 116    90 CONTINUE
 117       DO 92 I=N,1,-1
 118       Z=X(I)*D_ban(I)
 119       JJJJ=I+1
 120       DO 93 J=JJJJ,N
 121       Z=Z-A(J,I)*X(J)
 122    93 CONTINUE
 123       X(I)=Z
 124    92 CONTINUE
 125       return
 126       end subroutine BANAII
 127 !-----------------------------------------------------------------------------
 128       subroutine MATINVERT(N,NMAX,A,A1,osob)
 129
 130 !      implicit real*8 (a-h,o-z)
 131 !      include 'DIMENSIONS'
 132       integer :: N,NMAX
 133       real(kind=8),DIMENSION(NMAX,NMAX) :: A,A1
 134 !el      real(kind=8),DIMENSION(6*nres) :: D    !(MAXRES6) maxres6=6*maxres
 135 !el      COMMON /BANII/ D
 136       real(kind=8),DIMENSION(NMAX) :: X
 137       logical :: osob
 138       integer :: i,j
 139       DO I=1,N
 140         X(I)=0.0
 141       ENDDO
 142       X(1)=1.0
 143       CALL BANACH(N,NMAX,A,X,osob)
 144       if (osob) return
 145       DO I=1,N
 146         A1(I,1)=X(I)
 147       ENDDO
 148       DO I=2,N
 149         DO J=1,N
 150           X(J)=0.0
 151         ENDDO
 152         X(I)=1.0
 153         CALL BANAII(N,NMAX,A,X)
 154         DO J=1,N
 155           A1(J,I)=X(J)
 156         ENDDO
 157       ENDDO
 158       return
 159       end subroutine MATINVERT
 160 !-----------------------------------------------------------------------------
 161 ! bond_move.f
 162 !-----------------------------------------------------------------------------
 163       subroutine bond_move(nbond,nstart,psi,lprint,error)
 164
 165       use mcm_data, only:print_mc
 166       use geometry, only:alpha,beta,refsys,matmult
 167 ! Move NBOND fragment starting from the CA(nstart) by angle PSI.
 168 !      implicit real*8 (a-h,o-z)
 169 !      include 'DIMENSIONS'
 170       integer :: nbond,nstart
 171       real(kind=8) :: psi
 172       logical :: fail,error,lprint
 173 !      include 'COMMON.GEO'
 174 !      include 'COMMON.CHAIN'
 175 !      include 'COMMON.VAR'
 176 !      include 'COMMON.IOUNITS'
 177 !      include 'COMMON.MCM'
 178       real(kind=8),dimension(3) :: x,e1,e2,e3
 179       real(kind=8),dimension(3,3) :: e,rot,trans
 180       real(kind=8) :: cospsi,sinpsi,rij
 181       integer :: i,j,nend,i2,i3,i4,k
 182       error=.false.
 183       nend=nstart+nbond
 184       if (print_mc.gt.2) then
 185       write (iout,*) 'nstart=',nstart,' nend=',nend,' nbond=',nbond
 186       write (iout,*) 'psi=',psi
 187       write (iout,'(a)') 'Original coordinates of the fragment'
 188       do i=nstart,nend
 189         write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
 190       enddo
 191       endif
 192       if (nstart.lt.1 .or. nend .gt.nres .or. nbond.lt.2 .or. &
 193        nbond.ge.nres-1) then
 194         write (iout,'(a)') 'Bad data in BOND_MOVE.'
 195         error=.true.
 196         return
 197       endif
 198 ! Generate the reference system.
 199       i2=nend
 200       i3=nstart
 201       i4=nstart+1
 202       call refsys(i2,i3,i4,e1,e2,e3,error)
 203 ! Return, if couldn't define the reference system.
 204       if (error) return
 205 ! Compute the transformation matrix.
 206       cospsi=dcos(psi)
 207       sinpsi=dsin(psi)
 208       rot(1,1)=1.0D0
 209       rot(1,2)=0.0D0
 210       rot(1,3)=0.0D0
 211       rot(2,1)=0.0D0
 212       rot(2,2)=cospsi
 213       rot(2,3)=-sinpsi
 214       rot(3,1)=0.0D0
 215       rot(3,2)=sinpsi
 216       rot(3,3)=cospsi
 217       do i=1,3
 218         e(1,i)=e1(i)
 219         e(2,i)=e2(i)
 220         e(3,i)=e3(i)
 221       enddo
 222
 223       if (print_mc.gt.2) then
 224       write (iout,'(a)') 'Reference system and matrix r:'
 225       do i=1,3
 226         write(iout,'(i5,2(3f10.5,5x))')i,(e(i,j),j=1,3),(rot(i,j),j=1,3)
 227       enddo
 228       endif
 229
 230       call matmult(rot,e,trans)
 231       do i=1,3
 232         do j=1,3
 233           e(i,1)=e1(i)
 234           e(i,2)=e2(i)
 235           e(i,3)=e3(i)
 236         enddo
 237       enddo
 238       call matmult(e,trans,trans)
 239
 240       if (lprint) then
 241       write (iout,'(a)') 'The trans matrix:'
 242       do i=1,3
 243         write (iout,'(i5,3f10.5)') i,(trans(i,j),j=1,3)
 244       enddo
 245       endif
 246
 247       do i=nstart,nend
 248         do j=1,3
 249           rij=c(j,nstart)
 250           do k=1,3
 251             rij=rij+trans(j,k)*(c(k,i)-c(k,nstart))
 252           enddo
 253           x(j)=rij
 254         enddo
 255         do j=1,3
 256           c(j,i)=x(j)
 257         enddo
 258       enddo
 259
 260       if (lprint) then
 261       write (iout,'(a)') 'Rotated coordinates of the fragment'
 262       do i=nstart,nend
 263         write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
 264       enddo
 265       endif
 266
 267 !     call int_from_cart(.false.,lprint)
 268       if (nstart.gt.1) then
 269         theta(nstart+1)=alpha(nstart-1,nstart,nstart+1)
 270         phi(nstart+2)=beta(nstart-1,nstart,nstart+1,nstart+2)
 271         if (nstart.gt.2) phi(nstart+1)= &
 272                       beta(nstart-2,nstart-1,nstart,nstart+1)
 273       endif
 274       if (nend.lt.nres) then
 275         theta(nend+1)=alpha(nend-1,nend,nend+1)
 276         phi(nend+1)=beta(nend-2,nend-1,nend,nend+1)
 277         if (nend.lt.nres-1) phi(nend+2)= &
 278                       beta(nend-1,nend,nend+1,nend+2)
 279       endif
 280       if (print_mc.gt.2) then
 281       write (iout,'(/a,i3,a,i3,a/)') &
 282        'Moved internal coordinates of the ',nstart,'-',nend,&
 283        ' fragment:'
 284       do i=nstart+1,nstart+2
 285         write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
 286       enddo
 287       do i=nend+1,nend+2
 288         write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
 289       enddo
 290       endif
 291       return
 292       end subroutine bond_move
 293 !-----------------------------------------------------------------------------
 294 ! eigen.f
 295 !-----------------------------------------------------------------------------
 296 ! 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
 297 ! 31 MAR 94 - MWS - ADD A VARIABLE TO END OF MACHSW COMMON
 298 ! 26 JUN 93 - MWS - ETRED3: ADD RETURN FOR SPECIAL CASE N=1
 299 !  4 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
 300 ! 30 AUG 91 - MWS - JACDIA: LIMIT ITERATIONS, USE EPSLON IN TEST.
 301 ! 14 JUL 91 - MWS - JACOBI DIAGONALIZATION ALLOWS FOR LDVEC.NE.N
 302 ! 29 JAN 91 - TLW - GLDIAG: CHANGED COMMON DIAGSW TO MACHSW
 303 ! 29 OCT 90 - STE - FIX JACDIA UNDEFINED VARIABLE BUG
 304 ! 14 SEP 90 - MK  - NEW JACOBI DIAGONALIZATION (KDIAG=3)
 305 ! 27 MAR 88 - MWS - ALLOW FOR VECTOR ROUTINE IN GLDIAG
 306 ! 11 AUG 87 - MWS - SANITIZE CONSTANTS IN EQLRAT
 307 ! 15 FEB 87 - STE - FIX EINVIT SUB-MATRIX LOOP LIMIT
 308 !                   SCRATCH ARRAYS ARE N*8 REAL AND N INTEGER
 309 !  8 DEC 86 - STE - USE PERF INDEX FROM ESTPI1 TO JUDGE EINVIT FAILURE
 310 ! 30 NOV 86 - STE - DELETE LIGENB, MAKE EVVRSP DEFAULT
 311 !                   (GIVEIS FAILS ON CRAY FOR BENCHMC AND BENCHCI)
 312 !  7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
 313 ! 11 OCT 85 - STE - LIGENB,TQL2: USE DROT,DSWAP; TINVTB: SCALE VECTOR
 314 !                   BEFORE NORMALIZING; GENERIC FUNCTIONS
 315 ! 24 FEB 84 - STE - INITIALIZE INDEX ARRAY FOR LIGENB IN GLDIAG
 316 !  1 DEC 83 - STE - CHANGE MACHEP FROM 2**-54 TO 2**-50
 317 ! 28 SEP 82 - MWS - CONVERT TO IBM
 318 !
 319 !*MODULE EIGEN   *DECK EINVIT
 320       subroutine EINVIT(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
 321 !*
 322 !*    AUTHORS-
 323 !*       THIS IS A MODIFICATION OF TINVIT FROM EISPACK EDITION 3
 324 !*       DATED AUGUST 1983.
 325 !*       TINVIT IS A TRANSLATION OF THE INVERSE ITERATION TECHNIQUE
 326 !*       IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
 327 !*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
 328 !*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
 329 !*
 330 !*    PURPOSE -
 331 !*       THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
 332 !*       SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES.
 333 !*
 334 !*    METHOD -
 335 !*       INVERSE ITERATION.
 336 !*
 337 !*    ON ENTRY -
 338 !*       NM     - INTEGER
 339 !*                MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
 340 !*                ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
 341 !*                DIMENSION STATEMENT.
 342 !*       N      - INTEGER
 343 !*       D      - W.P. REAL (N)
 344 !*                CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
 345 !*       E      - W.P. REAL (N)
 346 !*                CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
 347 !*                IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
 348 !*       E2     - W.P. REAL (N)
 349 !*                CONTAINS THE SQUARES OF CORRESPONDING ELEMENTS OF E,
 350 !*                WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
 351 !*                E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
 352 !*                THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE
 353 !*                SUM OF THE MAGNITUDES OF D(I) AND D(I-1).  E2(1) MUST
 354 !*                CONTAIN 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER,
 355 !*                OR 2.0 IF THE EIGENVALUES ARE IN DESCENDING ORDER.
 356 !*                IF TQLRAT, BISECT, TRIDIB, OR IMTQLV
 357 !*                HAS BEEN USED TO FIND THE EIGENVALUES, THEIR
 358 !*                OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE.
 359 !*       M      - INTEGER
 360 !*                THE NUMBER OF SPECIFIED EIGENVECTORS.
 361 !*       W      - W.P. REAL (M)
 362 !*                CONTAINS THE M EIGENVALUES IN ASCENDING
 363 !*                OR DESCENDING ORDER.
 364 !*       IND    - INTEGER (M)
 365 !*                CONTAINS IN FIRST M POSITIONS THE SUBMATRIX INDICES
 366 !*                ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
 367 !*                1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX
 368 !*                FROM THE TOP, 2 FOR THOSE BELONGING TO THE SECOND
 369 !*                SUBMATRIX, ETC.
 370 !*       IERR   - INTEGER (LOGICAL UNIT NUMBER)
 371 !*                LOGICAL UNIT FOR ERROR MESSAGES
 372 !*
 373 !*    ON EXIT -
 374 !*       ALL INPUT ARRAYS ARE UNALTERED.
 375 !*       Z      - W.P. REAL (NM,M)
 376 !*                CONTAINS THE ASSOCIATED SET OF ORTHONORMAL
 377 !*                EIGENVECTORS. ANY VECTOR WHICH WHICH FAILS TO CONVERGE
 378 !*                IS LEFT AS IS (BUT NORMALIZED) WHEN ITERATING STOPPED.
 379 !*       IERR   - INTEGER
 380 !*                SET TO
 381 !*                ZERO    FOR NORMAL RETURN,
 382 !*                -R      IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
 383 !*                        EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS.
 384 !*                        (ONLY LAST FAILURE TO CONVERGE IS REPORTED)
 385 !*
 386 !*       RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
 387 !*
 388 !*       RV1    - W.P. REAL (N)
 389 !*                DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
 390 !*       RV2    - W.P. REAL (N)
 391 !*                SUPER(1)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
 392 !*       RV3    - W.P. REAL (N)
 393 !*                SUPER(2)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
 394 !*       RV4    - W.P. REAL (N)
 395 !*                ELEMENTS DEFINING L IN LU DECOMPOSITION
 396 !*       RV6    - W.P. REAL (N)
 397 !*                APPROXIMATE EIGENVECTOR
 398 !*
 399 !*    DIFFERENCES FROM EISPACK 3 -
 400 !*       EPS3 IS SCALED BY  EPSCAL  (ENHANCES CONVERGENCE, BUT
 401 !*          LOWERS ACCURACY)!
 402 !*       ONE MORE ITERATION (MINIMUM 2) IS PERFORMED AFTER CONVERGENCE
 403 !*          (ENHANCES ACCURACY)!
 404 !*       REPLACE LOOP WITH PYTHAG WITH SINGLE CALL TO DNRM2!
 405 !*       IF NOT CONVERGED, USE PERFORMANCE INDEX TO DECIDE ON ERROR
 406 !*          VALUE SETTING, BUT DO NOT STOP!
 407 !*       L.U. FOR ERROR MESSAGES PASSED THROUGH IERR
 408 !*       USE PARAMETER STATEMENTS AND GENERIC INTRINSIC FUNCTIONS
 409 !*       USE LEVEL 1 BLAS
 410 !*       USE IF-THEN-ELSE TO CLARIFY LOGIC
 411 !*       LOOP OVER SUBSPACES MADE INTO DO LOOP.
 412 !*       LOOP OVER INVERSE ITERATIONS MADE INTO DO LOOP
 413 !*       ZERO ONLY REQUIRED PORTIONS OF OUTPUT VECTOR
 414 !*
 415 !*    NOTE -
 416 !*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
 417 !*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
 418 !*
 419 !
 420       use comm_par
 421       LOGICAL :: CONVGD !el,GOPARR,DSKWRK,MASWRK
 422 !
 423       INTEGER :: GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG
 424       INTEGER :: IND(M)
 425 !
 426       real(kind=8),dimension(N) :: D,E2
 427       real(kind=8) :: E(*)!el E(L)
 428       real(kind=8) :: W(M),Z(NM,M)
 429       real(kind=8),dimension(N) :: RV1,RV2,RV3,RV4,RV6
 430       real(kind=8) :: ANORM,EPS2,EPS3,EPS4,NORM,ORDER,RHO,U,UK,V
 431       real(kind=8) :: X0,X1,XU
 432 !      real(kind=8) :: ESTPI1 !, DASUM, DDOT, DNRM2 EPSLON,
 433 !
 434 !el      integer :: ME,MASTER,NPROC,IBTYP,IPTIM
 435 !el      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
 436 !
 437       real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00, GRPTOL = 0.001D+00
 438       real(kind=8),PARAMETER :: EPSCAL = 0.5D+00, HUNDRD = 100.0D+00, TEN = 10.0D+00
 439 !
 440   001 FORMAT(' EIGENVECTOR ROUTINE EINVIT DID NOT CONVERGE FOR VECTOR' &
 441             ,I5,'.  NORM =',1P,E10.2,' PERFORMANCE INDEX =',E10.2/ &
 442             ' (AN ERROR HALT WILL OCCUR IF THE PI IS GREATER THAN 100)')
 443       integer :: LUEMSG
 444 !
 445 !-----------------------------------------------------------------------
 446 !
 447       LUEMSG = IERR
 448       IERR = 0
 449       X0 = ZERO
 450       UK = ZERO
 451       NORM = ZERO
 452       EPS2 = ZERO
 453       EPS3 = ZERO
 454       EPS4 = ZERO
 455       GROUP = 0
 456       TAG = 0
 457       ORDER = ONE - E2(1)
 458       Q = 0
 459       DO 930 SUBMAT = 1, N
 460          P = Q + 1
 461 !
 462 !        .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
 463 !
 464          DO 120 Q = P, N-1
 465             IF (E2(Q+1) .EQ. ZERO) GO TO 140
 466   120    CONTINUE
 467          Q = N
 468 !
 469 !        .......... FIND VECTORS BY INVERSE ITERATION ..........
 470 !
 471   140    CONTINUE
 472          TAG = TAG + 1
 473          ANORM = ZERO
 474          S = 0
 475 !
 476          DO 920 R = 1, M
 477             IF (IND(R) .NE. TAG) GO TO 920
 478             ITS = 1
 479             X1 = W(R)
 480             IF (S .NE. 0) GO TO 510
 481 !
 482 !        .......... CHECK FOR ISOLATED ROOT ..........
 483 !
 484             XU = ONE
 485             IF (P .EQ. Q) THEN
 486                RV6(P) = ONE
 487                CONVGD = .TRUE.
 488                GO TO 860
 489 !
 490             END IF
 491             NORM = ABS(D(P))
 492             DO 500 I = P+1, Q
 493                NORM = MAX( NORM, ABS(D(I)) + ABS(E(I)) )
 494   500       CONTINUE
 495 !
 496 !        .......... EPS2 IS THE CRITERION FOR GROUPING,
 497 !                   EPS3 REPLACES ZERO PIVOTS AND EQUAL
 498 !                   ROOTS ARE MODIFIED BY EPS3,
 499 !                   EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
 500 !
 501             EPS2 = GRPTOL * NORM
 502             EPS3 = EPSCAL * EPSLON(NORM)
 503             UK = Q - P + 1
 504             EPS4 = UK * EPS3
 505             UK = EPS4 / SQRT(UK)
 506             S = P
 507             GROUP = 0
 508             GO TO 520
 509 !
 510 !        .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
 511 !
 512   510       IF (ABS(X1-X0) .GE. EPS2) THEN
 513 !
 514 !                 ROOTS ARE SEPERATE
 515 !
 516                GROUP = 0
 517             ELSE
 518 !
 519 !                 ROOTS ARE CLOSE
 520 !
 521                GROUP = GROUP + 1
 522                IF (ORDER * (X1 - X0) .LE. EPS3) X1 = X0 + ORDER * EPS3
 523             END IF
 524 !
 525 !        .......... ELIMINATION WITH INTERCHANGES AND
 526 !                   INITIALIZATION OF VECTOR ..........
 527 !
 528   520       CONTINUE
 529 !
 530             U = D(P) - X1
 531             V = E(P+1)
 532             RV6(P) = UK
 533             DO 550 I = P+1, Q
 534                RV6(I) = UK
 535                IF (ABS(E(I)) .GT. ABS(U)) THEN
 536 !
 537 !                 EXCHANGE ROWS BEFORE ELIMINATION
 538 !
 539 !                  *** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
 540 !                      E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......
 541 !
 542                   XU = U / E(I)
 543                   RV4(I) = XU
 544                   RV1(I-1) = E(I)
 545                   RV2(I-1) = D(I) - X1
 546                   RV3(I-1) = E(I+1)
 547                   U = V - XU * RV2(I-1)
 548                   V = -XU * RV3(I-1)
 549 !
 550                ELSE
 551 !
 552 !                    STRAIGHT ELIMINATION
 553 !
 554                   XU = E(I) / U
 555                   RV4(I) = XU
 556                   RV1(I-1) = U
 557                   RV2(I-1) = V
 558                   RV3(I-1) = ZERO
 559                   U = D(I) - X1 - XU * V
 560                   V = E(I+1)
 561                END IF
 562   550       CONTINUE
 563 !
 564             IF (ABS(U) .LE. EPS3) U = EPS3
 565             RV1(Q) = U
 566             RV2(Q) = ZERO
 567             RV3(Q) = ZERO
 568 !
 569 !              DO INVERSE ITERATIONS
 570 !
 571             CONVGD = .FALSE.
 572             DO 800 ITS = 1, 5
 573                IF (ITS .EQ. 1) GO TO 600
 574 !
 575 !                    .......... FORWARD SUBSTITUTION ..........
 576 !
 577                   IF (NORM .EQ. ZERO) THEN
 578                      RV6(S) = EPS4
 579                      S = S + 1
 580                      IF (S .GT. Q) S = P
 581                   ELSE
 582                      XU = EPS4 / NORM
 583                      CALL DSCAL (Q-P+1, XU, RV6(P), 1)
 584                   END IF
 585 !
 586 !                     ... ELIMINATION OPERATIONS ON NEXT VECTOR
 587 !
 588                   DO 590 I = P+1, Q
 589                      U = RV6(I)
 590 !
 591 !                         IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
 592 !                         WAS PERFORMED EARLIER IN THE
 593 !                         TRIANGULARIZATION PROCESS ..........
 594 !
 595                      IF (RV1(I-1) .EQ. E(I)) THEN
 596                         U = RV6(I-1)
 597                         RV6(I-1) = RV6(I)
 598                      ELSE
 599                         U = RV6(I)
 600                      END IF
 601                      RV6(I) = U - RV4(I) * RV6(I-1)
 602   590             CONTINUE
 603   600          CONTINUE
 604 !
 605 !           .......... BACK SUBSTITUTION
 606 !
 607                RV6(Q) = RV6(Q) / RV1(Q)
 608                V = U
 609                U = RV6(Q)
 610                NORM = ABS(U)
 611                DO 620 I = Q-1, P, -1
 612                   RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
 613                   V = U
 614                   U = RV6(I)
 615                   NORM = NORM + ABS(U)
 616   620          CONTINUE
 617                IF (GROUP .EQ. 0) GO TO 700
 618 !
 619 !                 ....... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
 620 !                         MEMBERS OF GROUP ..........
 621 !
 622                   J = R
 623                   DO 680 JJ = 1, GROUP
 624   630                J = J - 1
 625                      IF (IND(J) .NE. TAG) GO TO 630
 626                      CALL DAXPY(Q-P+1, -DDOT(Q-P+1,RV6(P),1,Z(P,J),1),&
 627                                 Z(P,J),1,RV6(P),1)
 628   680             CONTINUE
 629                   NORM = DASUM(Q-P+1, RV6(P), 1)
 630   700          CONTINUE
 631 !
 632                IF (CONVGD) GO TO 840
 633                IF (NORM .GE. ONE) CONVGD = .TRUE.
 634   800       CONTINUE
 635 !
 636 !        .......... NORMALIZE SO THAT SUM OF SQUARES IS
 637 !                   1 AND EXPAND TO FULL ORDER ..........
 638 !
 639   840       CONTINUE
 640 !
 641             XU = ONE / DNRM2(Q-P+1,RV6(P),1)
 642 !
 643   860       CONTINUE
 644             DO 870 I = 1, P-1
 645                Z(I,R) = ZERO
 646   870       CONTINUE
 647             DO 890 I = P,Q
 648                Z(I,R) = RV6(I) * XU
 649   890       CONTINUE
 650             DO 900 I = Q+1, N
 651                Z(I,R) = ZERO
 652   900       CONTINUE
 653 !
 654             IF (.NOT.CONVGD) THEN
 655                RHO = ESTPI1(Q-P+1,X1,D(P),E(P),Z(P,R),ANORM)
 656                IF (RHO .GE. TEN .AND. LUEMSG .GT. 0 .AND. MASWRK) &
 657                    WRITE(LUEMSG,001) R,NORM,RHO
 658 !
 659 !               *** SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
 660 !
 661                IF (RHO .GT. HUNDRD) IERR = -R
 662             END IF
 663 !
 664             X0 = X1
 665   920    CONTINUE
 666 !
 667          IF (Q .EQ. N) GO TO 940
 668   930 CONTINUE
 669   940 CONTINUE
 670       return
 671       end subroutine EINVIT
 672 !-----------------------------------------------------------------------------
 673 !*MODULE EIGEN   *DECK ELAUM
 674       subroutine ELAU(HINV,L,D,A,E)
 675 !
 676       integer :: L,JL,JK,J,JM1,K
 677       real(kind=8) :: A(*)
 678       real(kind=8) :: D(L)
 679 !el      real(kind=8) :: E(L)
 680       real(kind=8) :: E(*)!el E(L)
 681       real(kind=8) :: F
 682       real(kind=8) :: G
 683       real(kind=8) :: HH
 684       real(kind=8) :: HINV
 685 !
 686       real(kind=8),PARAMETER :: ZERO = 0.0D+00, HALF = 0.5D+00
 687 !
 688       JL = L
 689       E(1) = A(1) * D(1)
 690       JK = 2
 691       DO 210 J = 2, JL
 692          F = D(J)
 693          G = ZERO
 694          JM1 = J - 1
 695 !
 696          DO 200 K = 1, JM1
 697             G = G + A(JK) * D(K)
 698             E(K) = E(K) + A(JK) * F
 699             JK = JK + 1
 700   200    CONTINUE
 701 !
 702          E(J) = G + A(JK) * F
 703          JK = JK + 1
 704   210 CONTINUE
 705 !
 706 !        .......... FORM P ..........
 707 !
 708       F = ZERO
 709       DO 245 J = 1, L
 710          E(J) = E(J) * HINV
 711          F = F + E(J) * D(J)
 712   245 CONTINUE
 713 !
 714 !     .......... FORM Q ..........
 715 !
 716       HH = F * HALF * HINV
 717       DO 250 J = 1, L
 718   250 E(J) = E(J) - HH * D(J)
 719 !
 720       return
 721       end subroutine ELAU
 722 !-----------------------------------------------------------------------------
 723 !*MODULE EIGEN   *DECK EPSLON
 724       real(kind=8) function EPSLON(X)
 725 !*
 726 !*    AUTHORS -
 727 !*       THIS ROUTINE WAS TAKEN FROM EISPACK EDITION 3 DATED 4/6/83
 728 !*       THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE NOV 1986
 729 !*
 730 !*    PURPOSE -
 731 !*       ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
 732 !*
 733 !*    ON ENTRY -
 734 !*       X      - WORKING PRECISION REAL
 735 !*                VALUES TO FIND EPSLON FOR
 736 !*
 737 !*    ON EXIT -
 738 !*       EPSLON - WORKING PRECISION REAL
 739 !*                SMALLEST POSITIVE VALUE SUCH THAT X+EPSLON .NE. ZERO
 740 !*
 741 !*    QUALIFICATIONS -
 742 !*       THIS ROUTINE SHOULD PERFORM PROPERLY ON ALL SYSTEMS
 743 !*       SATISFYING THE FOLLOWING TWO ASSUMPTIONS,
 744 !*          1.  THE BASE USED IN REPRESENTING FLOATING POINT
 745 !*              NUMBERS IS NOT A POWER OF THREE.
 746 !*          2.  THE QUANTITY  A  IN STATEMENT 10 IS REPRESENTED TO
 747 !*              THE ACCURACY USED IN FLOATING POINT VARIABLES
 748 !*              THAT ARE STORED IN MEMORY.
 749 !*       THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO
 750 !*       FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING
 751 !*       ASSUMPTION 2.
 752 !*       UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT,
 753 !*              A  IS NOT EXACTLY EQUAL TO FOUR-THIRDS,
 754 !*              B  HAS A ZERO FOR ITS LAST BIT OR DIGIT,
 755 !*              C  IS NOT EXACTLY EQUAL TO ONE,
 756 !*              EPS  MEASURES THE SEPARATION OF 1.0 FROM
 757 !*                   THE NEXT LARGER FLOATING POINT NUMBER.
 758 !*       THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED
 759 !*       ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD.
 760 !*
 761 !*    DIFFERENCES FROM EISPACK 3 -
 762 !*       USE IS MADE OF PARAMETER STATEMENTS AND INTRINSIC FUNCTIONS
 763 !*       --NO EXECUTEABLE CODE CHANGES--
 764 !*
 765 !*    NOTE -
 766 !*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
 767 !*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
 768 !
 769       real(kind=8) :: A,B,C,EPS,X
 770 !
 771       real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, FOUR=4.0D+00
 772 !
 773 !-----------------------------------------------------------------------
 774 !
 775       A = FOUR/THREE
 776    10 B = A - ONE
 777       C = B + B + B
 778       EPS = ABS(C - ONE)
 779       IF (EPS .EQ. ZERO) GO TO 10
 780       EPSLON = EPS*ABS(X)
 781       return
 782       end function EPSLON
 783 !-----------------------------------------------------------------------------
 784 !*MODULE EIGEN   *DECK EQLRAT
 785       subroutine EQLRAT(N,DIAG,E,E2IN,D,IND,IERR,E2)
 786 !*
 787 !*    AUTHORS -
 788 !*       THIS IS A MODIFICATION OF ROUTINE EQLRAT FROM EISPACK EDITION 3
 789 !*       DATED AUGUST 1983.
 790 !*       TQLRAT IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
 791 !*       ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
 792 !*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
 793 !*
 794 !*    PURPOSE -
 795 !*       THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
 796 !*       TRIDIAGONAL MATRIX
 797 !*
 798 !*    METHOD -
 799 !*       RATIONAL QL
 800 !*
 801 !*    ON ENTRY -
 802 !*       N      - INTEGER
 803 !*                THE ORDER OF THE MATRIX.
 804 !*       D      - W.P. REAL (N)
 805 !*                CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
 806 !*       E2     - W.P. REAL (N)
 807 !*                CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF
 808 !*                THE INPUT MATRIX IN ITS LAST N-1 POSITIONS.
 809 !*                E2(1) IS ARBITRARY.
 810 !*
 811 !*     ON EXIT -
 812 !*       D      - W.P. REAL (N)
 813 !*                CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
 814 !*                ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
 815 !*                ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
 816 !*                THE SMALLEST EIGENVALUES.
 817 !*       E2     - W.P. REAL (N)
 818 !*                DESTROYED.
 819 !*       IERR   - INTEGER
 820 !*                SET TO
 821 !*                ZERO       FOR NORMAL RETURN,
 822 !*                J          IF THE J-TH EIGENVALUE HAS NOT BEEN
 823 !*                           DETERMINED AFTER 30 ITERATIONS.
 824 !*
 825 !*    DIFFERENCES FROM EISPACK 3 -
 826 !*       G=G+B INSTEAD OF IF(G.EQ.0) G=B ; B=B/4
 827 !*       F77 BACKWARD LOOPS INSTEAD OF F66 CONSTRUCT
 828 !*       GENERIC INTRINSIC FUNCTIONS
 829 !*       ARRARY  IND  ADDED FOR USE BY EINVIT
 830 !*
 831 !*    NOTE -
 832 !*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
 833 !*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
 834 !
 835       INTEGER :: I,J,L,M,N,II,L1,IERR
 836       INTEGER,dimension(N) :: IND
 837 !
 838       real(kind=8),dimension(N) :: D,E,E2,DIAG,E2IN
 839       real(kind=8) :: B,C,F,G,H,P,R,S,T !,EPSLON
 840 !
 841       real(kind=8),PARAMETER :: ZERO = 0.0D+00, SCALE= 1.0D+00/64.0D+00, ONE = 1.0D+00
 842 !
 843       integer :: K,ITAG
 844 !-----------------------------------------------------------------------
 845       IERR = 0
 846       D(1)=DIAG(1)
 847       IND(1) = 1
 848       K = 0
 849       ITAG = 0
 850       IF (N .EQ. 1) GO TO 1001
 851 !
 852       DO 100 I = 2, N
 853          D(I)=DIAG(I)
 854   100 E2(I-1) = E2IN(I)
 855 !
 856       F = ZERO
 857       T = ZERO
 858       B = EPSLON(ONE)
 859       C = B *B
 860       B = B * SCALE
 861       E2(N) = ZERO
 862 !
 863       DO 290 L = 1, N
 864          H = ABS(D(L)) + ABS(E(L))
 865          IF (T .GE. H) GO TO 105
 866             T = H
 867             B = EPSLON(T)
 868             C = B * B
 869             B = B * SCALE
 870   105    CONTINUE
 871 !     .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
 872          M = L - 1
 873   110    M = M + 1
 874          IF (E2(M) .GT. C) GO TO 110
 875 !     .......... E2(N) IS ALWAYS ZERO, SO THERE IS AN EXIT
 876 !                FROM THE LOOP ..........
 877 !
 878          IF (M .LE. K) GO TO 125
 879             IF (M .NE. N) E2IN(M+1) = ZERO
 880             K = M
 881             ITAG = ITAG + 1
 882   125    CONTINUE
 883          IF (M .EQ. L) GO TO 210
 884 !
 885 !           ITERATE
 886 !
 887          DO 205 J = 1, 30
 888 !              .......... FORM SHIFT ..........
 889             L1 = L + 1
 890             S = SQRT(E2(L))
 891             G = D(L)
 892             P = (D(L1) - G) / (2.0D+00 * S)
 893             R = SQRT(P*P+1.0D+00)
 894             D(L) = S / (P + SIGN(R,P))
 895             H = G - D(L)
 896 !
 897             DO 140 I = L1, N
 898   140       D(I) = D(I) - H
 899 !
 900             F = F + H
 901 !              .......... RATIONAL QL TRANSFORMATION ..........
 902             G = D(M) + B
 903             H = G
 904             S = ZERO
 905             DO 200 I = M-1,L,-1
 906                P = G * H
 907                R = P + E2(I)
 908                E2(I+1) = S * R
 909                S = E2(I) / R
 910                D(I+1) = H + S * (H + D(I))
 911                G = D(I) - E2(I) / G   + B
 912                H = G * P / R
 913   200       CONTINUE
 914 !
 915             E2(L) = S * G
 916             D(L) = H
 917 !              .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST
 918             IF (H .EQ. ZERO) GO TO 210
 919             IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
 920             E2(L) = H * E2(L)
 921             IF (E2(L) .EQ. ZERO) GO TO 210
 922   205    CONTINUE
 923 !     .......... SET ERROR -- NO CONVERGENCE TO AN
 924 !                EIGENVALUE AFTER 30 ITERATIONS ..........
 925       IERR = L
 926       GO TO 1001
 927 !
 928 !           CONVERGED
 929 !
 930   210    P = D(L) + F
 931 !           .......... ORDER EIGENVALUES ..........
 932          I = 1
 933          IF (L .EQ. 1) GO TO 250
 934             IF (P .LT. D(1)) GO TO 230
 935                I = L
 936 !           .......... LOOP TO FIND ORDERED POSITION
 937   220          I = I - 1
 938                IF (P .LT. D(I)) GO TO 220
 939 !
 940                I = I + 1
 941                IF (I .EQ. L) GO TO 250
 942   230       CONTINUE
 943             DO 240 II = L, I+1, -1
 944                D(II) = D(II-1)
 945                IND(II) = IND(II-1)
 946   240       CONTINUE
 947 !
 948   250    CONTINUE
 949          D(I) = P
 950          IND(I) = ITAG
 951   290 CONTINUE
 952 !
 953  1001 return
 954       end subroutine EQLRAT
 955 !-----------------------------------------------------------------------------
 956 !*MODULE EIGEN   *DECK ESTPI1
 957       real(kind=8) function ESTPI1(N,EVAL,D,E,X,ANORM)
 958 !*
 959 !*    AUTHOR -
 960 !*       STEPHEN T. ELBERT (AMES LABORATORY-USDOE) DATE: 5 DEC 1986
 961 !*
 962 !*    PURPOSE -
 963 !*       EVALUATE SYMMETRIC TRIDIAGONAL MATRIX PERFORMANCE INDEX
 964 !*       *        *         *                  *           *
 965 !*       FOR 1 EIGENVECTOR
 966 !*           *
 967 !*
 968 !*    METHOD -
 969 !*       THIS ROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX A*X-X*EVAL
 970 !*       WHERE  A  IS A SYMMETRIC TRIDIAGONAL MATRIX STORED
 971 !*       IN THE DIAGONAL (D) AND SUB-DIAGONAL (E) VECTORS, EVAL IS THE
 972 !*       EIGENVALUE OF AN EIGENVECTOR OF  A,  NAMELY  X.
 973 !*       THIS NORM IS SCALED BY MACHINE ACCURACY FOR THE PROBLEM SIZE.
 974 !*       ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
 975 !*
 976 !*    ON ENTRY -
 977 !*       N      - INTEGER
 978 !*                THE ORDER OF THE MATRIX  A.
 979 !*       EVAL   - W.P. REAL
 980 !*                THE EIGENVALUE CORRESPONDING TO VECTOR  X.
 981 !*       D      - W.P. REAL (N)
 982 !*                THE DIAGONAL VECTOR OF  A.
 983 !*       E      - W.P. REAL (N)
 984 !*                THE SUB-DIAGONAL VECTOR OF  A.
 985 !*       X      - W.P. REAL (N)
 986 !*                AN EIGENVECTOR OF  A.
 987 !*       ANORM  - W.P. REAL
 988 !*                THE NORM OF  A  IF IT HAS BEEN PREVIOUSLY COMPUTED.
 989 !*
 990 !*    ON EXIT -
 991 !*       ANORM  - W.P. REAL
 992 !*                THE NORM OF  A, COMPUTED IF INITIALLY ZERO.
 993 !*       ESTPI1 - W.P. REAL
 994 !*          !!A*X-X*EVAL!! / (EPSLON(10*N)*!!A!!*!!X!!);
 995 !*          WHERE EPSLON(X) IS THE SMALLEST NUMBER SUCH THAT
 996 !*             X + EPSLON(X) .NE. X
 997 !*
 998 !*          ESTPI1 .LT. 1 == SATISFACTORY PERFORMANCE
 999 !*                 .GE. 1 AND .LE. 100 == MARGINAL PERFORMANCE
1000 !*                 .GT. 100 == POOR PERFORMANCE
1001 !*          (SEE LECT. NOTES IN COMP. SCI. VOL.6 PP 124-125)
1002 !
1003       integer :: N,I
1004       real(kind=8) :: ANORM,EVAL,RNORM,SIZE,XNORM
1005       real(kind=8),dimension(N) :: D,X
1006       real(kind=8) :: E(*)!el E(L)
1007 !      real(kind=8) :: EPSLON
1008 !
1009       real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00
1010 !
1011 !-----------------------------------------------------------------------
1012 !
1013       ESTPI1 = ZERO
1014       IF( N .LE. 1 ) RETURN
1015       SIZE = 10 * N
1016       IF (ANORM .EQ. ZERO) THEN
1017 !
1018 !              COMPUTE NORM OF  A
1019 !
1020          ANORM = MAX( ABS(D(1)) + ABS(E(2)), &
1021                       ABS(D(N)) + ABS(E(N)))
1022          DO 110 I = 2, N-1
1023             ANORM = MAX( ANORM, ABS(E(I))+ABS(D(I))+ABS(E(I+1)))
1024   110    CONTINUE
1025          IF(ANORM .EQ. ZERO) ANORM = ONE
1026       END IF
1027 !
1028 !           COMPUTE NORMS OF RESIDUAL AND EIGENVECTOR
1029 !
1030       XNORM = ABS(X(1)) + ABS(X(N))
1031       RNORM = ABS( (D(1)-EVAL)*X(1) + E(2)*X(2)) &
1032              +ABS( (D(N)-EVAL)*X(N) + E(N)*X(N-1))
1033       DO 120 I = 2, N-1
1034          XNORM = XNORM + ABS(X(I))
1035          RNORM = RNORM + ABS(E(I)*X(I-1) + (D(I)-EVAL)*X(I) &
1036                              + E(I+1)*X(I+1))
1037   120 CONTINUE
1038 !
1039       ESTPI1 = RNORM / (EPSLON(SIZE)*ANORM*XNORM)
1040       return
1041       end function ESTPI1
1042 !-----------------------------------------------------------------------------
1043 !*MODULE EIGEN   *DECK ETRBK3
1044       subroutine ETRBK3(NM,N,NV,A,M,Z)
1045 !*
1046 !*    AUTHORS-
1047 !*       THIS IS A MODIFICATION OF ROUTINE TRBAK3 FROM EISPACK EDITION 3
1048 !*       DATED AUGUST 1983.
1049 !*       EISPACK TRBAK3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
1050 !*       NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
1051 !*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
1052 !*       THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
1053 !*
1054 !*    PURPOSE -
1055 !*       THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
1056 !*       MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
1057 !*       SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  ETRED3.
1058 !*
1059 !*    METHOD -
1060 !*       THE CALCULATION IS CARRIED OUT BY FORMING THE MATRIX PRODUCT
1061 !*          Q*Z
1062 !*       WHERE  Q  IS A PRODUCT OF THE ORTHOGONAL SYMMETRIC MATRICES
1063 !*                Q = PROD(I)[1 - U(I)*.TRANSPOSE.U(I)*H(I)]
1064 !*       U  IS THE AUGMENTED SUB-DIAGONAL ROWS OF  A  AND
1065 !*       Z  IS THE SET OF EIGENVECTORS OF THE TRIDIAGONAL
1066 !*       MATRIX  F  WHICH WAS FORMED FROM THE ORIGINAL SYMMETRIC
1067 !*       MATRIX  C  BY THE SIMILARITY TRANSFORMATION
1068 !*                F = Q(TRANSPOSE) C Q
1069 !*       NOTE THAT ETRBK3 PRESERVES VECTOR EUCLIDEAN NORMS.
1070 !*
1071 !*
1072 !*    COMPLEXITY -
1073 !*       M*N**2
1074 !*
1075 !*    ON ENTRY-
1076 !*       NM     - INTEGER
1077 !*                MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1078 !*                ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
1079 !*                DIMENSION STATEMENT.
1080 !*       N      - INTEGER
1081 !*                THE ORDER OF THE MATRIX  A.
1082 !*       NV     - INTEGER
1083 !*                MUST BE SET TO THE DIMENSION OF THE ARRAY  A  AS
1084 !*                DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT.
1085 !*       A      - W.P. REAL (NV)
1086 !*                CONTAINS INFORMATION ABOUT THE ORTHOGONAL
1087 !*                TRANSFORMATIONS USED IN THE REDUCTION BY  ETRED3  IN
1088 !*                ITS FIRST  NV = N*(N+1)/2 POSITIONS.
1089 !*       M      - INTEGER
1090 !*                THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
1091 !*       Z      - W.P REAL (NM,M)
1092 !*                CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
1093 !*                IN ITS FIRST M COLUMNS.
1094 !*
1095 !*    ON EXIT-
1096 !*       Z      - W.P. REAL (NM,M)
1097 !*                CONTAINS THE TRANSFORMED EIGENVECTORS
1098 !*                IN ITS FIRST M COLUMNS.
1099 !*
1100 !*    DIFFERENCES WITH EISPACK 3 -
1101 !*       THE TWO INNER LOOPS ARE REPLACED BY DDOT AND DAXPY.
1102 !*       MULTIPLICATION USED INSTEAD OF DIVISION TO FIND S.
1103 !*       OUTER LOOP RANGE CHANGED FROM 2,N TO 3,N.
1104 !*       ADDRESS POINTERS FOR  A  SIMPLIFIED.
1105 !*
1106 !*    NOTE -
1107 !*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
1108 !*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
1109 !
1110       INTEGER :: I,II,IM1,IZ,J,M,N,NM,NV
1111 !
1112       real(kind=8) :: A(NV),Z(NM,M)
1113       real(kind=8) :: H,S !,DDOT
1114 !
1115       real(kind=8),PARAMETER :: ZERO = 0.0D+00
1116 !
1117 !-----------------------------------------------------------------------
1118 !
1119       IF (M .EQ. 0) RETURN
1120       IF (N .LE. 2) RETURN
1121 !
1122       II=3
1123       DO 140 I = 3, N
1124          IZ=II+1
1125          II=II+I
1126          H = A(II)
1127          IF (H .EQ. ZERO) GO TO 140
1128             IM1 = I - 1
1129             DO 130 J = 1, M
1130                S = -( DDOT(IM1,A(IZ),1,Z(1,J),1) * H) * H
1131                CALL DAXPY(IM1,S,A(IZ),1,Z(1,J),1)
1132   130       CONTINUE
1133   140 CONTINUE
1134       return
1135       end subroutine ETRBK3
1136 !-----------------------------------------------------------------------------
1137 !*MODULE EIGEN   *DECK ETRED3
1138       subroutine ETRED3(N,NV,A,D,E,E2)
1139 !*
1140 !*    AUTHORS -
1141 !*       THIS IS A MODIFICATION OF ROUTINE TRED3 FROM EISPACK EDITION 3
1142 !*       DATED AUGUST 1983.
1143 !*       EISPACK TRED3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
1144 !*       NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
1145 !*       HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
1146 !*       THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE JUN 1986
1147 !*
1148 !*    PURPOSE -
1149 !*       THIS ROUTINE REDUCES A REAL SYMMETRIC (PACKED) MATRIX, STORED
1150 !*       AS A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
1151 !*       USING ORTHOGONAL SIMILARITY TRANSFORMATIONS, PRESERVING THE
1152 !*       INFORMATION ABOUT THE TRANSFORMATIONS IN  A.
1153 !*
1154 !*    METHOD -
1155 !*       THE TRIDIAGONAL REDUCTION IS PERFORMED IN THE FOLLOWING WAY.
1156 !*       STARTING WITH J=N, THE ELEMENTS IN THE J-TH ROW TO THE
1157 !*       LEFT OF THE DIAGONAL ARE FIRST SCALED, TO AVOID POSSIBLE
1158 !*       UNDERFLOW IN THE TRANSFORMATION THAT MIGHT RESULT IN SEVERE
1159 !*       DEPARTURE FROM ORTHOGONALITY.  THE SUM OF SQUARES  SIGMA  OF
1160 !*       THESE SCALED ELEMENTS IS NEXT FORMED.  THEN, A VECTOR  U  AND
1161 !*       A SCALAR
1162 !*                      H = U(TRANSPOSE) * U / 2
1163 !*       DEFINE A REFLECTION OPERATOR
1164 !*                      P = I - U * U(TRANSPOSE) / H
1165 !*       WHICH IS ORTHOGONAL AND SYMMETRIC AND FOR WHICH THE
1166 !*       SIMILIARITY TRANSFORMATION  PAP  ELIMINATES THE ELEMENTS IN
1167 !*       THE J-TH ROW OF  A  TO THE LEFT OF THE SUBDIAGONAL AND THE
1168 !*       SYMMETRICAL ELEMENTS IN THE J-TH COLUMN.
1169 !*
1170 !*       THE NON-ZERO COMPONENTS OF  U  ARE THE ELEMENTS OF THE J-TH
1171 !*       ROW TO THE LEFT OF THE DIAGONAL WITH THE LAST OF THEM
1172 !*       AUGMENTED BY THE SQUARE ROOT OF  SIGMA  PREFIXED BY THE SIGN
1173 !*       OF THE SUBDIAGONAL ELEMENT.  BY STORING THE TRANSFORMED SUB-
1174 !*       DIAGONAL ELEMENT IN  E(J)  AND NOT OVERWRITING THE ROW
1175 !*       ELEMENTS ELIMINATED IN THE TRANSFORMATION, FULL INFORMATION
1176 !*       ABOUT  P  IS SAVE FOR LATER USE IN  ETRBK3.
1177 !*
1178 !*       THE TRANSFORMATION SETS  E2(J)  EQUAL TO  SIGMA  AND  E(J)
1179 !*       EQUAL TO THE SQUARE ROOT OF  SIGMA  PREFIXED BY THE SIGN
1180 !*       OF THE REPLACED SUBDIAGONAL ELEMENT.
1181 !*
1182 !*       THE ABOVE STEPS ARE REPEATED ON FURTHER ROWS OF THE
1183 !*       TRANSFORMED  A  IN REVERSE ORDER UNTIL  A  IS REDUCED TO TRI-
1184 !*       DIAGONAL FORM, THAT IS, REPEATED FOR  J = N-1,N-2,...,3.
1185 !*
1186 !*    COMPLEXITY -
1187 !*       2/3 N**3
1188 !*
1189 !*    ON ENTRY-
1190 !*       N      - INTEGER
1191 !*                THE ORDER OF THE MATRIX.
1192 !*       NV     - INTEGER
1193 !*                MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
1194 !*                AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT
1195 !*       A      - W.P. REAL (NV)
1196 !*                CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
1197 !*                INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
1198 !*                ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
1199 !*
1200 !*    ON EXIT-
1201 !*       A      - W.P. REAL (NV)
1202 !*                CONTAINS INFORMATION ABOUT THE ORTHOGONAL
1203 !*                TRANSFORMATIONS USED IN THE REDUCTION.
1204 !*       D      - W.P. REAL (N)
1205 !*                CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL
1206 !*                MATRIX.
1207 !*       E      - W.P. REAL (N)
1208 !*                CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
1209 !*                MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO
1210 !*       E2     - W.P. REAL (N)
1211 !*                CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF
1212 !*                E. MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
1213 !*
1214 !*    DIFFERENCES FROM EISPACK 3 -
1215 !*       OUTER LOOP CHANGED FROM II=1,N TO I=N,3,-1
1216 !*       PARAMETER STATEMENT AND GENERIC INTRINSIC FUNCTIONS USED
1217 !*       SCALE.NE.0 TEST NOW SPOTS TRI-DIAGONAL FORM
1218 !*       VALUES LESS THAN EPSLON CLEARED TO ZERO
1219 !*       USE BLAS(1)
1220 !*       U NOT COPIED TO D, LEFT IN A
1221 !*       E2 COMPUTED FROM E
1222 !*       INNER LOOPS SPLIT INTO ROUTINES ELAU AND FREDA
1223 !*       INVERSE OF H STORED INSTEAD OF H
1224 !*
1225 !*    NOTE -
1226 !*       QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
1227 !*       B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
1228 !
1229       INTEGER :: I,IIA,IZ0,L,N,NV
1230 !
1231       real(kind=8) :: A(NV),D(N),E2(N)
1232       real(kind=8) :: E(*)!el E(L)
1233       real(kind=8) :: AIIMAX,F,G,H,HROOT,SCALE,SCALEI
1234 !      real(kind=8) :: DASUM, DNRM2
1235 !
1236       real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00
1237 !
1238 !-----------------------------------------------------------------------
1239 !
1240       IF (N .LE. 2) GO TO 310
1241       IZ0 = (N*N+N)/2
1242       AIIMAX = ABS(A(IZ0))
1243       DO 300 I = N, 3, -1
1244          L = I - 1
1245          IIA = IZ0
1246          IZ0 = IZ0 - I
1247          AIIMAX = MAX(AIIMAX, ABS(A(IIA)))
1248          SCALE = DASUM (L, A(IZ0+1), 1)
1249          IF(SCALE .EQ. ABS(A(IIA-1)) .OR. AIIMAX+SCALE .EQ. AIIMAX) THEN
1250 !
1251 !           THIS ROW IS ALREADY IN TRI-DIAGONAL FORM
1252 !
1253             D(I) = A(IIA)
1254             IF (AIIMAX+D(I) .EQ. AIIMAX) D(I) = ZERO
1255             E(I) = A(IIA-1)
1256             IF (AIIMAX+E(I) .EQ. AIIMAX) E(I) = ZERO
1257             E2(I) = E(I)*E(I)
1258             A(IIA) = ZERO
1259             GO TO 300
1260 !
1261          END IF
1262 !
1263          SCALEI = ONE / SCALE
1264          CALL DSCAL(L,SCALEI,A(IZ0+1),1)
1265          HROOT = DNRM2(L,A(IZ0+1),1)
1266 !
1267          F = A(IZ0+L)
1268          G = -SIGN(HROOT,F)
1269          E(I) = SCALE * G
1270          E2(I) = E(I)*E(I)
1271          H = HROOT*HROOT - F * G
1272          A(IZ0+L) = F - G
1273          D(I) = A(IIA)
1274          A(IIA) = ONE / SQRT(H)
1275 !           .......... FORM P THEN Q IN E(1:L) ..........
1276          CALL ELAU(ONE/H,L,A(IZ0+1),A,E)
1277 !           .......... FORM REDUCED A ..........
1278          CALL FREDA(L,A(IZ0+1),A,E)
1279 !
1280   300 CONTINUE
1281   310 CONTINUE
1282       E(1) = ZERO
1283       E2(1)= ZERO
1284       D(1) = A(1)
1285       IF(N.EQ.1) RETURN
1286 !
1287       E(2) = A(2)
1288       E2(2)= A(2)*A(2)
1289       D(2) = A(3)
1290       return
1291       end subroutine ETRED3
1292 !-----------------------------------------------------------------------------
1293 !*MODULE EIGEN   *DECK EVVRSP
1294       subroutine EVVRSP(MSGFL,N,NVECT,LENA,NV,A,B,IND,ROOT,VECT,IORDER,IERR)
1295 !*
1296 !*    AUTHOR:  S. T. ELBERT, AMES LABORATORY-USDOE, JUNE 1985
1297 !*
1298 !*    PURPOSE -
1299 !*       FINDS   (ALL) EIGENVALUES    AND    (SOME OR ALL) EIGENVECTORS
1300 !*                     *    *                                   *
1301 !*       OF A REAL SYMMETRIC PACKED MATRIX.
1302 !*            *    *         *
1303 !*
1304 !*    METHOD -
1305 !*       THE METHOD AS PRESENTED IN THIS ROUTINE CONSISTS OF FOUR STEPS:
1306 !*       FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE
1307 !*       HOUSEHOLDER TECHNIQUE (ORTHOGONAL SIMILARITY TRANSFORMATIONS).
1308 !*       SECOND, THE ROOTS ARE LOCATED USING THE RATIONAL QL METHOD.
1309 !*       THIRD, THE VECTORS OF THE TRIDIAGONAL FORM ARE EVALUATED BY THE
1310 !*       INVERSE ITERATION TECHNIQUE.  VECTORS FOR DEGENERATE OR NEAR-
1311 !*       DEGENERATE ROOTS ARE FORCED TO BE ORTHOGONAL.
1312 !*       FOURTH, THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE
1313 !*       ORIGINAL ARRAY.
1314 !*
1315 !*       THESE ROUTINES ARE MODIFICATIONS OF THE EISPACK 3
1316 !*       ROUTINES TRED3, TQLRAT, TINVIT AND TRBAK3
1317 !*
1318 !*       FOR FURTHER DETAILS, SEE EISPACK USERS GUIDE, B. T. SMITH
1319 !*       ET AL, SPRINGER-VERLAG, LECTURE NOTES IN COMPUTER SCIENCE,
1320 !*       VOL. 6, 2-ND EDITION, 1976.  ANOTHER GOOD REFERENCE IS
1321 !*       THE SYMMETRIC EIGENVALUE PROBLEM BY B. N. PARLETT
1322 !*       PUBLISHED BY PRENTICE-HALL, INC., ENGLEWOOD CLIFFS, N.J. (1980)
1323 !*
1324 !*    ON ENTRY -
1325 !*       MSGFL  - INTEGER (LOGICAL UNIT NO.)
1326 !*                FILE WHERE ERROR MESSAGES WILL BE PRINTED.
1327 !*                IF MSGFL IS 0, ERROR MESSAGES WILL BE PRINTED ON LU 6.
1328 !*                IF MSGFL IS NEGATIVE, NO ERROR MESSAGES PRINTED.
1329 !*       N      - INTEGER
1330 !*                ORDER OF MATRIX A.
1331 !*       NVECT  - INTEGER
1332 !*                NUMBER OF VECTORS DESIRED.  0 .LE. NVECT .LE. N.
1333 !*       LENA   - INTEGER
1334 !*                DIMENSION OF  A  IN CALLING ROUTINE.  MUST NOT BE LESS
1335 !*                THAN (N*N+N)/2.
1336 !*       NV     - INTEGER
1337 !*                ROW DIMENSION OF VECT IN CALLING ROUTINE.   N .LE. NV.
1338 !*       A      - WORKING PRECISION REAL (LENA)
1339 !*                INPUT MATRIX, ROWS OF THE LOWER TRIANGLE PACKED INTO
1340 !*                LINEAR ARRAY OF DIMENSION N*(N+1)/2.  THE PACKED ORDER
1341 !*                IS A(1,1), A(2,1), A(2,2), A(3,1), A(3,2), ...
1342 !*       B      - WORKING PRECISION REAL (N,8)
1343 !*                SCRATCH ARRAY, 8*N ELEMENTS
1344 !*       IND    - INTEGER (N)
1345 !*                SCRATCH ARRAY OF LENGTH N.
1346 !*       IORDER - INTEGER
1347 !*                ROOT ORDERING FLAG.
1348 !*                = 0, ROOTS WILL BE PUT IN ASCENDING ORDER.
1349 !*                = 2, ROOTS WILL BE PUT IN DESCENDING ORDER.
1350 !*
1351 !*    ON EXIT -
1352 !*       A      - DESTORYED.  NOW HOLDS REFLECTION OPERATORS.
1353 !*       ROOT   - WORKING PRECISION REAL (N)
1354 !*                ALL EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
1355 !*                  IF IORDER = 0, ROOT(1) .LE. ... .LE. ROOT(N)
1356 !*                  IF IORDER = 2, ROOT(1) .GE. ... .GE. ROOT(N)
1357 !*       VECT   - WORKING PRECISION REAL (NV,NVECT)
1358 !*                EIGENVECTORS FOR ROOT(1), ..., ROOT(NVECT).
1359 !*       IERR   - INTEGER
1360 !*                = 0 IF NO ERROR DETECTED,
1361 !*                = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
1362 !*                = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
1363 !*                (FAILURES SHOULD BE VERY RARE.  CONTACT C. MOLER.)
1364 !*
1365 !
1366       use comm_par
1367 !el      LOGICAL :: GOPARR,DSKWRK,MASWRK
1368 !
1369       integer :: MSGFL,N,NVECT,LENA,NV,IORDER,IERR
1370       real(kind=8) :: A(LENA)
1371       real(kind=8) :: B(N,8)
1372       real(kind=8) :: ROOT(N)
1373       real(kind=8) :: T
1374       real(kind=8) :: VECT(NV,*)
1375 !
1376       INTEGER :: IND(N)
1377 !
1378 !el      integer :: ME,MASTER,NPROC,IBTYP,IPTIM
1379       real(kind=8) :: DSKW
1380
1381 !el   COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1382 !
1383   900 FORMAT(26H0*** EVVRSP PARAMETERS ***/ &
1384              14H ***      N = ,I8,4H ***/ &
1385              14H ***  NVECT = ,I8,4H ***/ &
1386              14H ***   LENA = ,I8,4H ***/ &
1387              14H ***     NV = ,I8,4H ***/ &
1388              14H *** IORDER = ,I8,4H ***/ &
1389              14H ***   IERR = ,I8,4H ***)
1390   901 FORMAT(37H VALUE OF LENA IS LESS THAN (N*N+N)/2)
1391   902 FORMAT(39H EQLRAT HAS FAILED TO CONVERGE FOR ROOT,I5)
1392   903 FORMAT(18H NV IS LESS THAN N)
1393   904 FORMAT(41H EINVIT HAS FAILED TO CONVERGE FOR VECTOR,I5)
1394   905 FORMAT(51H VALUE OF IORDER MUST BE 0 (SMALLEST ROOT FIRST) OR &
1395             ,23H 2 (LARGEST ROOT FIRST))
1396   906 FORMAT(' VALUE OF N IS LESS THAN OR EQUAL ZERO')
1397
1398       integer :: LMSGFL,I,J,L,JSV,KLIM,K
1399 !
1400 !-----------------------------------------------------------------------
1401 !
1402       LMSGFL=MSGFL
1403       IF (MSGFL .EQ. 0) LMSGFL=6
1404       IERR = N - 1
1405       IF (N .LE. 0) GO TO 800
1406       IERR = N + 1
1407       IF ( (N*N+N)/2 .GT. LENA) GO TO 810
1408 !
1409 !        REDUCE REAL SYMMETRIC MATRIX A TO TRIDIAGONAL FORM
1410 !
1411       CALL ETRED3(N,LENA,A,B(1,1),B(1,2),B(1,3))
1412 !
1413 !        FIND ALL EIGENVALUES OF TRIDIAGONAL MATRIX
1414 !
1415       CALL EQLRAT(N,B(1,1),B(1,2),B(1,3),ROOT,IND,IERR,B(1,4))
1416       IF (IERR .NE. 0) GO TO 820
1417 !
1418 !         CHECK THE DESIRED ORDER OF THE EIGENVALUES
1419 !
1420       B(1,3) = IORDER
1421       IF (IORDER .EQ. 0) GO TO 300
1422          IF (IORDER .NE. 2) GO TO 850
1423 !
1424 !         ORDER ROOTS IN DESCENDING ORDER (LARGEST FIRST)...
1425 !        TURN ROOT AND IND ARRAYS END FOR END
1426 !
1427          DO 210 I = 1, N/2
1428             J = N+1-I
1429             T = ROOT(I)
1430             ROOT(I) = ROOT(J)
1431             ROOT(J) = T
1432             L = IND(I)
1433             IND(I) = IND(J)
1434             IND(J) = L
1435   210    CONTINUE
1436 !
1437 !           FIND I AND J MARKING THE START AND END OF A SEQUENCE
1438 !           OF DEGENERATE ROOTS
1439 !
1440          I=0
1441   220    CONTINUE
1442             I = I+1
1443             IF (I .GT. N) GO TO 300
1444             DO 230 J=I,N
1445                IF (ROOT(J) .NE. ROOT(I)) GO TO 240
1446   230       CONTINUE
1447             J = N+1
1448   240       CONTINUE
1449             J = J-1
1450             IF (J .EQ. I) GO TO 220
1451 !
1452 !                    TURN AROUND IND BETWEEN I AND J
1453 !
1454             JSV = J
1455             KLIM = (J-I+1)/2
1456             DO 250 K=1,KLIM
1457                L = IND(J)
1458                IND(J) = IND(I)
1459                IND(I) = L
1460                I = I+1
1461                J = J-1
1462   250       CONTINUE
1463             I = JSV
1464          GO TO 220
1465 !
1466   300 CONTINUE
1467 !
1468       IF (NVECT .LE. 0) RETURN
1469       IF (NV .LT. N) GO TO 830
1470 !
1471 !        FIND EIGENVECTORS OF TRI-DIAGONAL MATRIX VIA INVERSE ITERATION
1472 !
1473       IERR = LMSGFL
1474       CALL EINVIT(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,IND,&
1475                   VECT,IERR,B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
1476       IF (IERR .NE. 0) GO TO 840
1477 !
1478 !        FIND EIGENVECTORS OF SYMMETRIC MATRIX VIA BACK TRANSFORMATION
1479 !
1480   400 CONTINUE
1481       CALL ETRBK3(NV,N,LENA,A,NVECT,VECT)
1482       RETURN
1483 !
1484 !        ERROR MESSAGE SECTION
1485 !
1486   800 IF (LMSGFL .LT. 0) RETURN
1487       IF (MASWRK) WRITE(LMSGFL,906)
1488       GO TO 890
1489 !
1490   810 IF (LMSGFL .LT. 0) RETURN
1491       IF (MASWRK) WRITE(LMSGFL,901)
1492       GO TO 890
1493 !
1494   820 IF (LMSGFL .LT. 0) RETURN
1495       IF (MASWRK) WRITE(LMSGFL,902) IERR
1496       GO TO 890
1497 !
1498   830 IF (LMSGFL .LT. 0) RETURN
1499       IF (MASWRK) WRITE(LMSGFL,903)
1500       GO TO 890
1501 !
1502   840 CONTINUE
1503       IF ((LMSGFL .GT. 0).AND.MASWRK) WRITE(LMSGFL,904) -IERR
1504       GO TO 400
1505 !
1506   850 IERR=-1
1507       IF (LMSGFL .LT. 0) RETURN
1508       IF (MASWRK) WRITE(LMSGFL,905)
1509       GO TO 890
1510 !
1511   890 CONTINUE
1512       IF (MASWRK) WRITE(LMSGFL,900) N,NVECT,LENA,NV,IORDER,IERR
1513       return
1514       end subroutine EVVRSP
1515 !-----------------------------------------------------------------------------
1516 !*MODULE EIGEN   *DECK FREDA
1517       subroutine FREDA(L,D,A,E)
1518 !
1519       integer :: l,jk,j,k
1520       real(kind=8) :: A(*)
1521       real(kind=8) :: D(L)
1522       real(kind=8) :: E(*)!el E(L)
1523       real(kind=8) :: F
1524       real(kind=8) :: G
1525 !
1526       JK = 1
1527 !
1528 !     .......... FORM REDUCED A ..........
1529 !
1530       DO 280 J = 1, L
1531          F = D(J)
1532          G = E(J)
1533 !
1534          DO 260 K = 1, J
1535             A(JK) = A(JK) - F * E(K) - G * D(K)
1536             JK = JK + 1
1537   260    CONTINUE
1538 !
1539   280 CONTINUE
1540       return
1541       end subroutine FREDA
1542 !-----------------------------------------------------------------------------
1543 !*MODULE EIGEN   *DECK GIVEIS
1544       subroutine GIVEIS(N,NVECT,NV,A,B,INDB,ROOT,VECT,IERR)
1545
1546 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1547       integer :: N,NVECT,NV,IERR
1548       real(kind=8),DIMENSION(*) :: A
1549       real(kind=8),DIMENSION(N,8) :: B
1550       integer,DIMENSION(N) :: INDB
1551       real(kind=8),DIMENSION(N) :: ROOT
1552       real(kind=8),DIMENSION(NV,NVECT) :: VECT
1553 !
1554 !     EISPACK-BASED SUBSTITUTE FOR QCPE ROUTINE GIVENS.
1555 !     FINDS ALL EIGENVALUES AND SOME EIGENVECTORS OF A REAL SYMMETRIC
1556 !     MATRIX.   AUTHOR.. C. MOLER AND D. SPANGLER, N.R.C.C., 4/1/79.
1557 !
1558 !     INPUT..
1559 !     N     = ORDER OF MATRIX .
1560 !     NVECT = NUMBER OF VECTORS DESIRED.  0 .LE. NVECT .LE. N .
1561 !     NV    = LEADING DIMENSION OF VECT .
1562 !     A     = INPUT MATRIX, COLUMNS OF THE UPPER TRIANGLE PACKED INTO
1563 !             LINEAR ARRAY OF DIMENSION N*(N+1)/2 .
1564 !     B     = SCRATCH ARRAY, 8*N ELEMENTS (NOTE THIS IS MORE THAN
1565 !             PREVIOUS VERSIONS OF GIVENS.)
1566 !    IND    = INDEX ARRAY OF N ELEMENTS
1567 !
1568 !     OUTPUT..
1569 !     A       DESTROYED .
1570 !     ROOT  = ALL EIGENVALUES, ROOT(1) .LE. ... .LE. ROOT(N) .
1571 !             (FOR OTHER ORDERINGS, SEE BELOW.)
1572 !     VECT  = EIGENVECTORS FOR ROOT(1),..., ROOT(NVECT) .
1573 !     IERR  = 0 IF NO ERROR DETECTED,
1574 !           = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
1575 !           = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
1576 !             (FAILURES SHOULD BE VERY RARE.  CONTACT MOLER.)
1577 !
1578 !     CALLS MODIFIED EISPACK ROUTINES TRED3B, IMTQLV, TINVTB, AND
1579 !     TRBK3B.  THE ROUTINES TRED3B, TINVTB, AND TRBK3B.
1580 !     THE ORIGINAL EISPACK ROUTINES TRED3, TINVIT, AND TRBAK3
1581 !     WERE MODIFIED BY THE INTRODUCTION OF TWO ROUTINES FROM THE
1582 !     BLAS LIBRARY - DDOT AND DAXPY.
1583 !
1584 !         IF TINVIT FAILS TO CONVERGE, TQL2 IS CALLED
1585 !
1586 !         SEE EISPACK USERS GUIDE, B. T. SMITH ET AL, SPRINGER-VERLAG
1587 !     LECTURE NOTES IN COMPUTER SCIENCE, VOL. 6, 2-ND EDITION, 1976 .
1588 !     NOTE THAT IMTQLV AND TINVTB HAVE INTERNAL MACHINE
1589 !     DEPENDENT CONSTANTS.
1590 !
1591 !el      DATA ONE, ZERO /1.0D+00, 0.0D+00/
1592       real(kind=8) :: ZERO = 0.0D+00, ONE = 1.0D+00
1593
1594       integer :: i,j
1595
1596       CALL TRED3B(N,(N*N+N)/2,A,B(1,1),B(1,2),B(1,3))
1597       CALL IMTQLV(N,B(1,1),B(1,2),B(1,3),ROOT,INDB,IERR,B(1,4))
1598       IF (IERR .NE. 0) RETURN
1599 !
1600 !     TO REORDER ROOTS...
1601 !     K = N/2
1602 !     B(1,3) = 2.0D+00
1603 !     DO 50 I = 1, K
1604 !        J = N+1-I
1605 !        T = ROOT(I)
1606 !        ROOT(I) = ROOT(J)
1607 !        ROOT(J) = T
1608 ! 50  CONTINUE
1609 !
1610       IF (NVECT .LE. 0) RETURN
1611       CALL TINVTB(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,INDB,VECT,IERR,&
1612            B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
1613       IF (IERR .EQ. 0) GO TO 160
1614 !
1615 !      IF INVERSE ITERATION GIVES AN ERROR IN DETERMINING THE
1616 !      EIGENVECTORS, TRY THE QL ALGORITHM IF ALL THE EIGENVECTORS
1617 !      ARE DESIRED.
1618 !
1619       IF (NVECT .NE. N) RETURN
1620       DO 120 I = 1, NVECT
1621       DO 100 J = 1, N
1622       VECT(I,J) = ZERO
1623   100 CONTINUE
1624       VECT(I,I) = ONE
1625   120 CONTINUE
1626       CALL TQL2 (NV,N,B(1,1),B(1,2),VECT,IERR)
1627       DO 140 I = 1, NVECT
1628       ROOT(I) = B(I,1)
1629   140 CONTINUE
1630       IF (IERR .NE. 0) RETURN
1631   160 CALL TRBK3B(NV,N,(N*N+N)/2,A,NVECT,VECT)
1632       return
1633       end subroutine GIVEIS
1634 !-----------------------------------------------------------------------------
1635 !*MODULE EIGEN   *DECK GLDIAG
1636       subroutine GLDIAG(LDVECT,NVECT,N,H,WRK,EIG,VECTOR,IERR,IWRK)
1637 !
1638 !      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1639 !
1640       use comm_iofile
1641       use comm_machsw
1642       use comm_par
1643 !el      LOGICAL :: GOPARR,DSKWRK,MASWRK
1644 !
1645       integer :: LDVECT,NVECT,N,IERR
1646       real(kind=8),DIMENSION(*) :: H
1647       real(kind=8),DIMENSION(N,8) :: WRK
1648       real(kind=8),DIMENSION(N) :: EIG
1649       integer,DIMENSION(N) :: IWRK
1650       real(kind=8),DIMENSION(LDVECT,NVECT) :: VECTOR
1651 !
1652 !el      integer :: IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
1653 !el      integer :: KDIAG,ICORFL,IXDR
1654 !el      COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA
1655 !el      COMMON /MACHSW/ KDIAG,ICORFL,IXDR
1656 !el      integer :: ME,MASTER,NPROC,IBTYP,IPTIM
1657 !el      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1658 !
1659       integer :: LENH,KORDER
1660
1661 !     ----- GENERAL ROUTINE TO DIAGONALIZE A SYMMETRIC MATRIX -----
1662 !     IF KDIAG = 0, USE A ROUTINE FROM THE VECTOR LIBRARY,
1663 !                   IF AVAILABLE (SEE THE SUBROUTINE 'GLDIAG'
1664 !                   IN VECTOR.SRC), OR EVVRSP OTHERWISE
1665 !              = 1, USE EVVRSP
1666 !              = 2, USE GIVEIS
1667 !              = 3, USE JACOBI
1668 !
1669 !           N      = DIMENSION (ORDER) OF MATRIX TO BE SOLVED
1670 !           LDVECT = LEADING DIMENSION OF VECTOR
1671 !           NVECT  = NUMBER OF VECTORS DESIRED
1672 !           H      = MATRIX TO BE DIAGONALIZED
1673 !           WRK    = N*8 W.P. REAL WORDS OF SCRATCH SPACE
1674 !           EIG    = EIGENVECTORS (OUTPUT)
1675 !           VECTOR = EIGENVECTORS (OUTPUT)
1676 !           IERR   = ERROR FLAG (OUTPUT)
1677 !           IWRK   = N INTEGER WORDS OF SCRATCH SPACE
1678 !
1679       IERR = 0
1680 !
1681 !         ----- USE STEVE ELBERT'S ROUTINE -----
1682 !
1683       IF(KDIAG.LE.1  .OR.  KDIAG.GT.3) THEN
1684          LENH = (N*N+N)/2
1685          KORDER =0
1686          CALL EVVRSP(IW,N,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR, &
1687                      KORDER,IERR)
1688       END IF
1689 !
1690 !         ----- USE MODIFIED EISPAK ROUTINE -----
1691 !
1692       IF(KDIAG.EQ.2) &
1693          CALL GIVEIS(N,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
1694 !
1695 !         ----- USE JACOBI ROTATION ROUTINE -----
1696 !
1697       IF(KDIAG.EQ.3) THEN
1698          IF(NVECT.EQ.N) THEN
1699             CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,N)
1700          ELSE
1701             IF (MASWRK) WRITE(IW,9000) N,NVECT,LDVECT
1702             CALL ABRT
1703          END IF
1704       END IF
1705       RETURN
1706 !
1707  9000 FORMAT(1X,'IN -GLDIAG-, N,NVECT,LDVECT=',3I8/ &
1708              1X,'THE JACOBI CODE CANNOT COPE WITH N.NE.NVECT!'/ &
1709              1X,'SO THIS RUN DOES NOT PERMIT KDIAG=3.')
1710       end subroutine GLDIAG
1711 !-----------------------------------------------------------------------------
1712 !*MODULE EIGEN   *DECK IMTQLV
1713       subroutine IMTQLV(N,D,E,E2,W,IND,IERR,RV1)
1714
1715 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1716       INTEGER :: N,TAG,IERR
1717       real(kind=8) :: MACHEP
1718       real(kind=8),DIMENSION(N) :: D,E2,W,RV1
1719       real(kind=8) :: E(*)!el E(L)
1720       integer,DIMENSION(N) :: IND
1721       integer :: k,i,l,j,m,mml,ii
1722       real(kind=8) :: c,p,s,f,b,g,r
1723 !
1724 !     THIS ROUTINE IS A VARIANT OF  IMTQL1  WHICH IS A TRANSLATION OF
1725 !     ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND
1726 !     WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
1727 !     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
1728 !
1729 !     THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
1730 !     MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM
1731 !     THEIR CORRESPONDING SUBMATRIX INDICES.
1732 !
1733 !     ON INPUT-
1734 !
1735 !        N IS THE ORDER OF THE MATRIX,
1736 !
1737 !        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1738 !
1739 !        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1740 !          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
1741 !
1742 !        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
1743 !          E2(1) IS ARBITRARY.
1744 !
1745 !     ON OUTPUT-
1746 !
1747 !        D AND E ARE UNALTERED,
1748 !
1749 !        ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
1750 !          AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
1751 !          MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
1752 !          E2(1) IS ALSO SET TO ZERO,
1753 !
1754 !        W CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
1755 !          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
1756 !          ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
1757 !          THE SMALLEST EIGENVALUES,
1758 !
1759 !        IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE
1760 !          CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES
1761 !          BELONGING TO THE FIRST SUBMATRIX FROM THE TOP,
1762 !          2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.,
1763 !
1764 !        IERR IS SET TO
1765 !          ZERO       FOR NORMAL RETURN,
1766 !          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
1767 !                     DETERMINED AFTER 30 ITERATIONS,
1768 !
1769 !        RV1 IS A TEMPORARY STORAGE ARRAY.
1770 !
1771 !     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1772 !     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1773 !
1774 !     ------------------------------------------------------------------
1775 !
1776 !     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1777 !                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1778 !
1779 !                **********
1780       MACHEP = 2.0D+00**(-50)
1781 !
1782       IERR = 0
1783       K = 0
1784       TAG = 0
1785 !
1786       DO 100 I = 1, N
1787       W(I) = D(I)
1788       IF (I .NE. 1) RV1(I-1) = E(I)
1789   100 CONTINUE
1790 !
1791       E2(1) = 0.0D+00
1792       RV1(N) = 0.0D+00
1793 !
1794       DO 360 L = 1, N
1795       J = 0
1796 !     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
1797   120 DO 140 M = L, N
1798       IF (M .EQ. N) GO TO 160
1799       IF (ABS(RV1(M)) .LE. MACHEP * (ABS(W(M)) + ABS(W(M+1)))) GO TO &
1800            160
1801 !     ********** GUARD AGAINST UNDERFLOWED ELEMENT OF E2 **********
1802       IF (E2(M+1) .EQ. 0.0D+00) GO TO 180
1803   140 CONTINUE
1804 !
1805   160 IF (M .LE. K) GO TO 200
1806       IF (M .NE. N) E2(M+1) = 0.0D+00
1807   180 K = M
1808       TAG = TAG + 1
1809   200 P = W(L)
1810       IF (M .EQ. L) GO TO 280
1811       IF (J .EQ. 30) GO TO 380
1812       J = J + 1
1813 !     ********** FORM SHIFT **********
1814       G = (W(L+1) - P) / (2.0D+00 * RV1(L))
1815       R = SQRT(G*G+1.0D+00)
1816       G = W(M) - P + RV1(L) / (G + SIGN(R,G))
1817       S = 1.0D+00
1818       C = 1.0D+00
1819       P = 0.0D+00
1820       MML = M - L
1821 !     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
1822       DO 260 II = 1, MML
1823       I = M - II
1824       F = S * RV1(I)
1825       B = C * RV1(I)
1826       IF (ABS(F) .LT. ABS(G)) GO TO 220
1827       C = G / F
1828       R = SQRT(C*C+1.0D+00)
1829       RV1(I+1) = F * R
1830       S = 1.0D+00 / R
1831       C = C * S
1832       GO TO 240
1833   220 S = F / G
1834       R = SQRT(S*S+1.0D+00)
1835       RV1(I+1) = G * R
1836       C = 1.0D+00 / R
1837       S = S * C
1838   240 G = W(I+1) - P
1839       R = (W(I) - G) * S + 2.0D+00 * C * B
1840       P = S * R
1841       W(I+1) = G + P
1842       G = C * R - B
1843   260 CONTINUE
1844 !
1845       W(L) = W(L) - P
1846       RV1(L) = G
1847       RV1(M) = 0.0D+00
1848       GO TO 120
1849 !     ********** ORDER EIGENVALUES **********
1850   280 IF (L .EQ. 1) GO TO 320
1851 !     ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
1852       DO 300 II = 2, L
1853       I = L + 2 - II
1854       IF (P .GE. W(I-1)) GO TO 340
1855       W(I) = W(I-1)
1856       IND(I) = IND(I-1)
1857   300 CONTINUE
1858 !
1859   320 I = 1
1860   340 W(I) = P
1861       IND(I) = TAG
1862   360 CONTINUE
1863 !
1864       GO TO 400
1865 !     ********** SET ERROR -- NO CONVERGENCE TO AN
1866 !                EIGENVALUE AFTER 30 ITERATIONS **********
1867   380 IERR = L
1868   400 return
1869 !     ********** LAST CARD OF IMTQLV **********
1870       end subroutine IMTQLV
1871 !-----------------------------------------------------------------------------
1872 !*MODULE EIGEN   *DECK JACDG
1873       subroutine JACDG(A,VEC,EIG,JBIG,BIG,LDVEC,N)
1874 !
1875 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1876 !
1877       integer :: LDVEC,N
1878       real(kind=8),DIMENSION(*) :: A
1879       real(kind=8),DIMENSION(LDVEC,N) :: VEC
1880       real(kind=8),DIMENSION(N) :: EIG,BIG
1881       integer,DIMENSION(N) :: JBIG
1882 !
1883       real(kind=8),PARAMETER :: ONE=1.0D+00
1884       integer :: i,NB1,NB2,NMIN,NMAX
1885 !
1886 !     ----- JACOBI DIAGONALIZATION OF SYMMETRIC MATRIX -----
1887 !     SYMMETRIC MATRIX -A- OF DIMENSION -N- IS DESTROYED ON EXIT.
1888 !     ALL EIGENVECTORS ARE FOUND, SO -VEC- MUST BE SQUARE,
1889 !     UNLESS SOMEONE TAKES THE TROUBLE TO LOOK AT -NMAX- BELOW.
1890 !     -BIG- AND -JBIG- ARE SCRATCH WORK ARRAYS.
1891 !
1892       CALL VCLR(VEC,1,LDVEC*N)
1893       DO 20 I = 1,N
1894         VEC(I,I) = ONE
1895    20 CONTINUE
1896 !
1897       NB1 = N
1898       NB2 = (NB1*NB1+NB1)/2
1899       NMIN = 1
1900       NMAX = NB1
1901 !
1902       CALL JACDIA(A,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1903 !
1904       DO 30 I=1,N
1905         EIG(I) = A((I*I+I)/2)
1906    30 CONTINUE
1907 !
1908       CALL JACORD(VEC,EIG,NB1,LDVEC)
1909       return
1910       end subroutine JACDG
1911 !-----------------------------------------------------------------------------
1912 !*MODULE EIGEN   *DECK JACDIA
1913       subroutine JACDIA(F,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1914
1915 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1916       use comm_par
1917       integer :: NB1,NB2,LDVEC,NMIN,NMAX
1918 !el      LOGICAL :: GOPARR,DSKWRK,MASWRK
1919       real(kind=8),DIMENSION(NB2) :: F
1920       real(kind=8),DIMENSION(LDVEC,NB1) :: VEC
1921       real(kind=8),DIMENSION(NB1) :: BIG
1922       integer,DIMENSION(NB1) :: JBIG
1923 !
1924 !el      integer :: ME,MASTER,NPROC,IBTYP,IPTIM
1925 !el      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1926 !
1927       real(kind=8),PARAMETER :: ROOT2=0.707106781186548D+00
1928       real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00, D1050=1.05D+00,&
1929                  D1500=1.5D+00, D3875=3.875D+00,&
1930                  D0500=0.5D+00, D1375=1.375D+00, D0250=0.25D+00
1931       real(kind=8),PARAMETER :: C2=1.0D-12, C3=4.0D-16,&
1932                  C4=2.0D-16, C5=8.0D-09, C6=3.0D-06
1933       integer :: i,ii,j,k,jj,IEAA,IEAB,MAXIT,ITER,I1,IA,IB,IAA,IBB,IEAR,&
1934             IEBR,IR,IT,KQ,IR1
1935       real(kind=8) :: T,TT,EPS,SD,TEST,DIF,CX,SX,T2X2,T2X25,T1,T2
1936 !
1937 !      F IS THE MATRIX TO BE DIAGONALIZED, F IS STORED TRIANGULAR
1938 !      VEC IS THE ARRAY OF EIGENVECTORS, DIMENSION NB1*NB1
1939 !      BIG AND JBIG ARE TEMPORARY SCRATCH AREAS OF DIMENSION NB1
1940 !      THE ROTATIONS AMONG THE FIRST NMIN BASIS FUNCTIONS ARE NOT
1941 !      ACCOUNTED FOR.
1942 !      THE ROTATIONS AMONG THE LAST NB1-NMAX BASIS FUNCTIONS ARE NOT
1943 !      ACCOUNTED FOR.
1944 !
1945       IEAA=0
1946       IEAB=0
1947       TT=ZERO
1948       EPS = 64.0D+00*EPSLON(ONE)
1949 !
1950 !      LOOP OVER COLUMNS (K) OF TRIANGULAR MATRIX TO DETERMINE
1951 !      LARGEST OFF-DIAGONAL ELEMENTS IN ROW(I).
1952 !
1953       DO 20 I=1,NB1
1954          BIG(I)=ZERO
1955          JBIG(I)=0
1956          IF(I.LT.NMIN  .OR.  I.EQ.1) GO TO 20
1957          II = (I*I-I)/2
1958          J=MIN(I-1,NMAX)
1959          DO 10 K=1,J
1960             IF(ABS(BIG(I)).GE.ABS(F(II+K))) GO TO 10
1961             BIG(I)=F(II+K)
1962             JBIG(I)=K
1963    10    CONTINUE
1964    20 CONTINUE
1965 !
1966 !     ----- 2X2 JACOBI ITERATIONS BEGIN HERE -----
1967 !
1968       MAXIT=MAX(NB2*20,500)
1969       ITER=0
1970    30 CONTINUE
1971       ITER=ITER+1
1972 !
1973 !      FIND SMALLEST DIAGONAL ELEMENT
1974 !
1975       SD=D1050
1976       JJ=0
1977       DO 40 J=1,NB1
1978          JJ=JJ+J
1979          SD= MIN(SD,ABS(F(JJ)))
1980    40 CONTINUE
1981       TEST = MAX(EPS, C2*MAX(SD,C6))
1982 !
1983 !      FIND LARGEST OFF-DIAGONAL ELEMENT
1984 !
1985       T=ZERO
1986       I1=MAX(2,NMIN)
1987       IB = I1
1988       DO 50 I=I1,NB1
1989          IF(T.GE.ABS(BIG(I))) GO TO 50
1990          T= ABS(BIG(I))
1991          IB=I
1992    50 CONTINUE
1993 !
1994 !      TEST FOR CONVERGENCE, THEN DETERMINE ROTATION.
1995 !
1996       IF(T.LT.TEST) RETURN
1997 !                   ******
1998 !
1999       IF(ITER.GT.MAXIT) THEN
2000          IF (MASWRK) THEN
2001             WRITE(6,*) 'JACOBI DIAGONALIZATION FAILS, DIMENSION=',NB1
2002             WRITE(6,9020) ITER,T,TEST,SD
2003          ENDIF
2004          CALL ABRT
2005          STOP
2006       END IF
2007 !
2008       IA=JBIG(IB)
2009       IAA=IA*(IA-1)/2
2010       IBB=IB*(IB-1)/2
2011       DIF=F(IAA+IA)-F(IBB+IB)
2012       IF(ABS(DIF).GT.C3*T) GO TO 70
2013       SX=ROOT2
2014       CX=ROOT2
2015       GO TO 110
2016    70 T2X2=BIG(IB)/DIF
2017       T2X25=T2X2*T2X2
2018       IF(T2X25 .GT. C4) GO TO 80
2019       CX=ONE
2020       SX=T2X2
2021       GO TO 110
2022    80 IF(T2X25 .GT. C5) GO TO 90
2023       SX=T2X2*(ONE-D1500*T2X25)
2024       CX=ONE-D0500*T2X25
2025       GO TO 110
2026    90 IF(T2X25 .GT. C6) GO TO 100
2027       CX=ONE+T2X25*(T2X25*D1375 - D0500)
2028       SX= T2X2*(ONE + T2X25*(T2X25*D3875 - D1500))
2029       GO TO 110
2030   100 T=D0250  / SQRT(D0250   + T2X25)
2031       CX= SQRT(D0500   + T)
2032       SX= SIGN( SQRT(D0500   - T),T2X2)
2033   110 IEAR=IAA+1
2034       IEBR=IBB+1
2035 !
2036       DO 230 IR=1,NB1
2037          T=F(IEAR)*SX
2038          F(IEAR)=F(IEAR)*CX+F(IEBR)*SX
2039          F(IEBR)=T-F(IEBR)*CX
2040          IF(IR-IA) 220,120,130
2041   120    TT=F(IEBR)
2042          IEAA=IEAR
2043          IEAB=IEBR
2044          F(IEBR)=BIG(IB)
2045          IEAR=IEAR+IR-1
2046          IF(JBIG(IR)) 200,220,200
2047   130    T=F(IEAR)
2048          IT=IA
2049          IEAR=IEAR+IR-1
2050          IF(IR-IB) 180,150,160
2051   150    F(IEAA)=F(IEAA)*CX+F(IEAB)*SX
2052          F(IEAB)=TT*CX+F(IEBR)*SX
2053          F(IEBR)=TT*SX-F(IEBR)*CX
2054          IEBR=IEBR+IR-1
2055          GO TO 200
2056   160    IF(  ABS(T) .GE.  ABS(F(IEBR))) GO TO 170
2057          IF(IB.GT.NMAX) GO TO 170
2058          T=F(IEBR)
2059          IT=IB
2060   170    IEBR=IEBR+IR-1
2061   180    IF(  ABS(T) .LT.  ABS(BIG(IR))) GO TO 190
2062          BIG(IR) = T
2063          JBIG(IR) = IT
2064          GO TO 220
2065   190    IF(IA .NE. JBIG(IR) .AND. IB .NE. JBIG(IR))  GO TO 220
2066   200    KQ=IEAR-IR-IA+1
2067          BIG(IR)=ZERO
2068          IR1=MIN(IR-1,NMAX)
2069          DO 210 I=1,IR1
2070             K=KQ+I
2071             IF(ABS(BIG(IR)) .GE. ABS(F(K)))  GO TO 210
2072             BIG(IR) = F(K)
2073             JBIG(IR)=I
2074   210    CONTINUE
2075   220    IEAR=IEAR+1
2076   230    IEBR=IEBR+1
2077 !
2078       DO 240 I=1,NB1
2079          T1=VEC(I,IA)*CX + VEC(I,IB)*SX
2080          T2=VEC(I,IA)*SX - VEC(I,IB)*CX
2081          VEC(I,IA)=T1
2082          VEC(I,IB)=T2
2083   240 CONTINUE
2084       GO TO 30
2085 !
2086  9020 FORMAT(1X,'ITER=',I6,' T,TEST,SD=',1P,3E20.10)
2087       end subroutine JACDIA
2088 !-----------------------------------------------------------------------------
2089 !*MODULE EIGEN   *DECK JACORD
2090       subroutine JACORD(VEC,EIG,N,LDVEC)
2091
2092 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2093       integer :: N,LDVEC
2094       real(kind=8),DIMENSION(LDVEC,N) :: VEC
2095       real(kind=8),DIMENSION(N) :: EIG
2096       integer :: i,jj,j
2097       real(kind=8) :: T
2098 !
2099 !     ---- SORT EIGENDATA INTO ASCENDING ORDER -----
2100 !
2101       DO 290 I = 1, N
2102          JJ = I
2103          DO 270 J = I, N
2104             IF (EIG(J) .LT. EIG(JJ)) JJ = J
2105   270    CONTINUE
2106          IF (JJ .EQ. I) GO TO 290
2107          T = EIG(JJ)
2108          EIG(JJ) = EIG(I)
2109          EIG(I) = T
2110          DO 280 J = 1, N
2111             T = VEC(J,JJ)
2112             VEC(J,JJ) = VEC(J,I)
2113             VEC(J,I) = T
2114   280    CONTINUE
2115   290 CONTINUE
2116       return
2117       end subroutine JACORD
2118 !-----------------------------------------------------------------------------
2119 !*MODULE EIGEN   *DECK TINVTB
2120       subroutine TINVTB(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
2121
2122 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2123       integer :: NM,N,M,IERR
2124       real(kind=8),DIMENSION(N) :: D,E,E2
2125       real(kind=8),DIMENSION(M) :: W
2126       real(kind=8),DIMENSION(NM,M) :: Z
2127       real(kind=8),DIMENSION(N) :: RV1,RV2,RV3,RV4,RV6
2128       integer,DIMENSION(M) :: IND
2129       real(kind=8) :: MACHEP,NORM
2130       INTEGER :: P,Q,R,S,TAG,GROUP
2131       integer :: ii,j,jj,i,iqmp,its
2132       real(kind=8) :: ORDER,XU,UK,X0,U,EPS2,EPS3,EPS4,x1,v
2133
2134 !     ------------------------------------------------------------------
2135 !
2136 !     THIS ROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH-
2137 !     NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
2138 !     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
2139 !
2140 !     THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
2141 !     SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
2142 !     USING INVERSE ITERATION.
2143 !
2144 !     ON INPUT-
2145 !
2146 !        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2147 !          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2148 !          DIMENSION STATEMENT,
2149 !
2150 !        N IS THE ORDER OF THE MATRIX,
2151 !
2152 !        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2153 !
2154 !        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2155 !          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
2156 !
2157 !        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E,
2158 !          WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
2159 !          E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
2160 !          THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM
2161 !          OF THE MAGNITUDES OF D(I) AND D(I-1).  E2(1) MUST CONTAIN
2162 !          0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0
2163 !          IF THE EIGENVALUES ARE IN DESCENDING ORDER.  IF  BISECT,
2164 !          TRIDIB, OR  IMTQLV  HAS BEEN USED TO FIND THE EIGENVALUES,
2165 !          THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE,
2166 !
2167 !        M IS THE NUMBER OF SPECIFIED EIGENVALUES,
2168 !
2169 !        W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER,
2170 !
2171 !        IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
2172 !          ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
2173 !          1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
2174 !          THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.
2175 !
2176 !     ON OUTPUT-
2177 !
2178 !        ALL INPUT ARRAYS ARE UNALTERED,
2179 !
2180 !        Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
2181 !          ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO,
2182 !
2183 !        IERR IS SET TO
2184 !          ZERO       FOR NORMAL RETURN,
2185 !          -R         IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
2186 !                     EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS,
2187 !
2188 !        RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
2189 !
2190 !     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2191 !     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2192 !
2193 !     ------------------------------------------------------------------
2194 !
2195 !     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2196 !                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2197 !
2198 !                **********
2199       MACHEP = 2.0D+00**(-50)
2200 !
2201       IERR = 0
2202       IF (M .EQ. 0) GO TO 680
2203       TAG = 0
2204       ORDER = 1.0D+00 - E2(1)
2205       XU = 0.0D+00
2206       UK = 0.0D+00
2207       X0 = 0.0D+00
2208       U  = 0.0D+00
2209       EPS2 = 0.0D+00
2210       EPS3 = 0.0D+00
2211       EPS4 = 0.0D+00
2212       GROUP = 0
2213       Q = 0
2214 !     ********** ESTABLISH AND PROCESS NEXT SUBMATRIX **********
2215   100 P = Q + 1
2216       IP = P + 1
2217 !
2218       DO 120 Q = P, N
2219       IF (Q .EQ. N) GO TO 140
2220       IF (E2(Q+1) .EQ. 0.0D+00) GO TO 140
2221   120 CONTINUE
2222 !     ********** FIND VECTORS BY INVERSE ITERATION **********
2223   140 TAG = TAG + 1
2224       IQMP = Q - P + 1
2225       S = 0
2226 !
2227       DO 660 R = 1, M
2228       IF (IND(R) .NE. TAG) GO TO 660
2229       ITS = 1
2230       X1 = W(R)
2231       IF (S .NE. 0) GO TO 220
2232 !     ********** CHECK FOR ISOLATED ROOT **********
2233       XU = 1.0D+00
2234       IF (P .NE. Q) GO TO 160
2235       RV6(P) = 1.0D+00
2236       GO TO 600
2237   160 NORM = ABS(D(P))
2238 !
2239       DO 180 I = IP, Q
2240   180 NORM = NORM + ABS(D(I)) + ABS(E(I))
2241 !     ********** EPS2 IS THE CRITERION FOR GROUPING,
2242 !                EPS3 REPLACES ZERO PIVOTS AND EQUAL
2243 !                ROOTS ARE MODIFIED BY EPS3,
2244 !                EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW **********
2245       EPS2 = 1.0D-03 * NORM
2246       EPS3 = MACHEP * NORM
2247       UK = IQMP
2248       EPS4 = UK * EPS3
2249       UK = EPS4 / SQRT(UK)
2250       S = P
2251   200 GROUP = 0
2252       GO TO 240
2253 !     ********** LOOK FOR CLOSE OR COINCIDENT ROOTS **********
2254   220 IF (ABS(X1-X0) .GE. EPS2) GO TO 200
2255       GROUP = GROUP + 1
2256       IF (ORDER * (X1 - X0) .LE. 0.0D+00) X1 = X0 + ORDER * EPS3
2257 !     ********** ELIMINATION WITH INTERCHANGES AND
2258 !                INITIALIZATION OF VECTOR **********
2259   240 V = 0.0D+00
2260 !
2261       DO 300 I = P, Q
2262       RV6(I) = UK
2263       IF (I .EQ. P) GO TO 280
2264       IF (ABS(E(I)) .LT. ABS(U)) GO TO 260
2265 !     ********** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
2266 !                E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY **********
2267       XU = U / E(I)
2268       RV4(I) = XU
2269       RV1(I-1) = E(I)
2270       RV2(I-1) = D(I) - X1
2271       RV3(I-1) = 0.0D+00
2272       IF (I .NE. Q) RV3(I-1) = E(I+1)
2273       U = V - XU * RV2(I-1)
2274       V = -XU * RV3(I-1)
2275       GO TO 300
2276   260 XU = E(I) / U
2277       RV4(I) = XU
2278       RV1(I-1) = U
2279       RV2(I-1) = V
2280       RV3(I-1) = 0.0D+00
2281   280 U = D(I) - X1 - XU * V
2282       IF (I .NE. Q) V = E(I+1)
2283   300 CONTINUE
2284 !
2285       IF (U .EQ. 0.0D+00) U = EPS3
2286       RV1(Q) = U
2287       RV2(Q) = 0.0D+00
2288       RV3(Q) = 0.0D+00
2289 !     ********** BACK SUBSTITUTION
2290 !                FOR I=Q STEP -1 UNTIL P DO -- **********
2291   320 DO 340 II = P, Q
2292       I = P + Q - II
2293       RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
2294       V = U
2295       U = RV6(I)
2296   340 CONTINUE
2297 !     ********** ORTHOGONALIZE WITH RESPECT TO PREVIOUS
2298 !                MEMBERS OF GROUP **********
2299       IF (GROUP .EQ. 0) GO TO 400
2300       J = R
2301 !
2302       DO 380 JJ = 1, GROUP
2303   360 J = J - 1
2304       IF (IND(J) .NE. TAG) GO TO 360
2305       XU = DDOT(IQMP,RV6(P),1,Z(P,J),1)
2306 !
2307       CALL DAXPY(IQMP,-XU,Z(P,J),1,RV6(P),1)
2308 !
2309   380 CONTINUE
2310 !
2311   400 NORM = 0.0D+00
2312 !
2313       DO 420 I = P, Q
2314   420 NORM = NORM + ABS(RV6(I))
2315 !
2316       IF (NORM .GE. 1.0D+00) GO TO 560
2317 !     ********** FORWARD SUBSTITUTION **********
2318       IF (ITS .EQ. 5) GO TO 540
2319       IF (NORM .NE. 0.0D+00) GO TO 440
2320       RV6(S) = EPS4
2321       S = S + 1
2322       IF (S .GT. Q) S = P
2323       GO TO 480
2324   440 XU = EPS4 / NORM
2325 !
2326       DO 460 I = P, Q
2327   460 RV6(I) = RV6(I) * XU
2328 !     ********** ELIMINATION OPERATIONS ON NEXT VECTOR
2329 !                ITERATE **********
2330   480 DO 520 I = IP, Q
2331       U = RV6(I)
2332 !     ********** IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
2333 !                WAS PERFORMED EARLIER IN THE
2334 !                TRIANGULARIZATION PROCESS **********
2335       IF (RV1(I-1) .NE. E(I)) GO TO 500
2336       U = RV6(I-1)
2337       RV6(I-1) = RV6(I)
2338   500 RV6(I) = U - RV4(I) * RV6(I-1)
2339   520 CONTINUE
2340 !
2341       ITS = ITS + 1
2342       GO TO 320
2343 !     ********** SET ERROR -- NON-CONVERGED EIGENVECTOR **********
2344   540 IERR = -R
2345       XU = 0.0D+00
2346       GO TO 600
2347 !     ********** NORMALIZE SO THAT SUM OF SQUARES IS
2348 !                1 AND EXPAND TO FULL ORDER **********
2349   560 U = 0.0D+00
2350 !
2351       DO 580 I = P, Q
2352       RV6(I) = RV6(I) / NORM
2353   580 U = U + RV6(I)**2
2354 !
2355       XU = 1.0D+00 / SQRT(U)
2356 !
2357   600 DO 620 I = 1, N
2358   620 Z(I,R) = 0.0D+00
2359 !
2360       DO 640 I = P, Q
2361   640 Z(I,R) = RV6(I) * XU
2362 !
2363       X0 = X1
2364   660 CONTINUE
2365 !
2366       IF (Q .LT. N) GO TO 100
2367   680 return
2368 !     ********** LAST CARD OF TINVIT **********
2369       end subroutine TINVTB
2370 !-----------------------------------------------------------------------------
2371 !*MODULE EIGEN   *DECK TQL2
2372       subroutine TQL2(NM,N,D,E,Z,IERR)
2373
2374 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2375       integer :: NM,N,IERR
2376       real(kind=8) :: MACHEP
2377       real(kind=8),DIMENSION(N) :: D!,E
2378       real(kind=8) :: E(*)!el E(L)
2379       real(kind=8),DIMENSION(NM,N) :: Z
2380       integer :: ii,i,j,mml,m,l1,k,l
2381       real(kind=8) :: c,f,b,h,g,p,r,s
2382 !
2383 !     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
2384 !     NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
2385 !     WILKINSON.
2386 !     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
2387 !
2388 !     THIS ROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
2389 !     OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
2390 !     THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
2391 !     BE FOUND IF  TRED2  HAS BEEN USED TO REDUCE THIS
2392 !     FULL MATRIX TO TRIDIAGONAL FORM.
2393 !
2394 !     ON INPUT-
2395 !
2396 !        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2397 !          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2398 !          DIMENSION STATEMENT,
2399 !
2400 !        N IS THE ORDER OF THE MATRIX,
2401 !
2402 !        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2403 !
2404 !        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2405 !          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY,
2406 !
2407 !        Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
2408 !          REDUCTION BY  TRED2, IF PERFORMED.  IF THE EIGENVECTORS
2409 !          OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
2410 !          THE IDENTITY MATRIX.
2411 !
2412 !      ON OUTPUT-
2413 !
2414 !        D CONTAINS THE EIGENVALUES IN ASCENDING ORDER.  IF AN
2415 !          ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
2416 !          UNORDERED FOR INDICES 1,2,...,IERR-1,
2417 !
2418 !        E HAS BEEN DESTROYED,
2419 !
2420 !        Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
2421 !          TRIDIAGONAL (OR FULL) MATRIX.  IF AN ERROR EXIT IS MADE,
2422 !          Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
2423 !          EIGENVALUES,
2424 !
2425 !        IERR IS SET TO
2426 !          ZERO       FOR NORMAL RETURN,
2427 !          J          IF THE J-TH EIGENVALUE HAS NOT BEEN
2428 !                     DETERMINED AFTER 30 ITERATIONS.
2429 !
2430 !     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2431 !     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2432 !
2433 !     ------------------------------------------------------------------
2434 !
2435 !     ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2436 !                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2437 !
2438 !                **********
2439       MACHEP = 2.0D+00**(-50)
2440 !
2441       IERR = 0
2442       IF (N .EQ. 1) GO TO 400
2443 !
2444       DO 100 I = 2, N
2445   100 E(I-1) = E(I)
2446 !
2447       F = 0.0D+00
2448       B = 0.0D+00
2449       E(N) = 0.0D+00
2450 !
2451       DO 300 L = 1, N
2452       J = 0
2453       H = MACHEP * (ABS(D(L)) + ABS(E(L)))
2454       IF (B .LT. H) B = H
2455 !     ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
2456       DO 120 M = L, N
2457       IF (ABS(E(M)) .LE. B) GO TO 140
2458 !     ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
2459 !                THROUGH THE BOTTOM OF THE LOOP **********
2460   120 CONTINUE
2461 !
2462   140 IF (M .EQ. L) GO TO 280
2463   160 IF (J .EQ. 30) GO TO 380
2464       J = J + 1
2465 !     ********** FORM SHIFT **********
2466       L1 = L + 1
2467       G = D(L)
2468       P = (D(L1) - G) / (2.0D+00 * E(L))
2469       R = SQRT(P*P+1.0D+00)
2470       D(L) = E(L) / (P + SIGN(R,P))
2471       H = G - D(L)
2472 !
2473       DO 180 I = L1, N
2474   180 D(I) = D(I) - H
2475 !
2476       F = F + H
2477 !     ********** QL TRANSFORMATION **********
2478       P = D(M)
2479       C = 1.0D+00
2480       S = 0.0D+00
2481       MML = M - L
2482 !     ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
2483       DO 260 II = 1, MML
2484       I = M - II
2485       G = C * E(I)
2486       H = C * P
2487       IF (ABS(P) .LT. ABS(E(I))) GO TO 200
2488       C = E(I) / P
2489       R = SQRT(C*C+1.0D+00)
2490       E(I+1) = S * P * R
2491       S = C / R
2492       C = 1.0D+00 / R
2493       GO TO 220
2494   200 C = P / E(I)
2495       R = SQRT(C*C+1.0D+00)
2496       E(I+1) = S * E(I) * R
2497       S = 1.0D+00 / R
2498       C = C * S
2499   220 P = C * D(I) - S * G
2500       D(I+1) = H + S * (C * G + S * D(I))
2501 !     ********** FORM VECTOR **********
2502       CALL DROT(N,Z(1,I+1),1,Z(1,I),1,C,S)
2503 !
2504   260 CONTINUE
2505 !
2506       E(L) = S * P
2507       D(L) = C * P
2508       IF (ABS(E(L)) .GT. B) GO TO 160
2509   280 D(L) = D(L) + F
2510   300 CONTINUE
2511 !     ********** ORDER EIGENVALUES AND EIGENVECTORS **********
2512       DO 360 II = 2, N
2513       I = II - 1
2514       K = I
2515       P = D(I)
2516 !
2517       DO 320 J = II, N
2518       IF (D(J) .GE. P) GO TO 320
2519       K = J
2520       P = D(J)
2521   320 CONTINUE
2522 !
2523       IF (K .EQ. I) GO TO 360
2524       D(K) = D(I)
2525       D(I) = P
2526 !
2527       CALL DSWAP(N,Z(1,I),1,Z(1,K),1)
2528 !
2529   360 CONTINUE
2530 !
2531       GO TO 400
2532 !     ********** SET ERROR -- NO CONVERGENCE TO AN
2533 !                EIGENVALUE AFTER 30 ITERATIONS **********
2534   380 IERR = L
2535   400 return
2536 !     ********** LAST CARD OF TQL2 **********
2537       end subroutine TQL2
2538 !-----------------------------------------------------------------------------
2539 !*MODULE EIGEN   *DECK TRBK3B
2540       subroutine TRBK3B(NM,N,NV,A,M,Z)
2541
2542 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2543       integer :: NM,N,NV,M
2544       real(kind=8),DIMENSION(NV) :: A
2545       real(kind=8),DIMENSION(NM,M) :: Z
2546       integer :: i,l,iz,ik,j
2547       real(kind=8) :: h,s
2548 !
2549 !     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
2550 !     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
2551 !     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
2552 !
2553 !     THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
2554 !     MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
2555 !     SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY  TRED3B.
2556 !
2557 !     ON INPUT-
2558 !
2559 !        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2560 !          ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2561 !          DIMENSION STATEMENT,
2562 !
2563 !        N IS THE ORDER OF THE MATRIX,
2564 !
2565 !        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2566 !          AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2567 !
2568 !        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS
2569 !          USED IN THE REDUCTION BY  TRED3B IN ITS FIRST
2570 !          N*(N+1)/2 POSITIONS,
2571 !
2572 !        M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
2573 !
2574 !        Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
2575 !          IN ITS FIRST M COLUMNS.
2576 !
2577 !     ON OUTPUT-
2578 !
2579 !        Z CONTAINS THE TRANSFORMED EIGENVECTORS
2580 !          IN ITS FIRST M COLUMNS.
2581 !
2582 !     NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
2583 !
2584 !     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2585 !     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2586 !
2587 !     ------------------------------------------------------------------
2588 !
2589       IF (M .EQ. 0) GO TO 140
2590       IF (N .EQ. 1) GO TO 140
2591 !
2592       DO 120 I = 2, N
2593       L = I - 1
2594       IZ = (I * L) / 2
2595       IK = IZ + I
2596       H = A(IK)
2597       IF (H .EQ. 0.0D+00) GO TO 120
2598 !
2599       DO 100 J = 1, M
2600       S = -DDOT(L,A(IZ+1),1,Z(1,J),1)
2601 !
2602 !     ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
2603       S = (S / H) / H
2604 !
2605       CALL DAXPY(L,S,A(IZ+1),1,Z(1,J),1)
2606 !
2607   100 CONTINUE
2608 !
2609   120 CONTINUE
2610 !
2611   140 return
2612 !     ********** LAST CARD OF TRBAK3 **********
2613       end subroutine TRBK3B
2614 !-----------------------------------------------------------------------------
2615 !*MODULE EIGEN   *DECK TRED3B
2616       subroutine TRED3B(N,NV,A,D,E,E2)
2617
2618 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2619       integer :: N,NV
2620       real(kind=8),DIMENSION(NV) :: A
2621       real(kind=8),DIMENSION(N) :: D,E2
2622       real(kind=8) :: E(*)!el E(L)
2623       integer :: ii,i,l,iz,k,jk,j,jm1
2624       real(kind=8) :: h,f,g,scale,dt,hh
2625 !
2626 !     THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
2627 !     NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
2628 !     HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
2629 !
2630 !     THIS ROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS
2631 !     A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
2632 !     USING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
2633 !
2634 !     ON INPUT-
2635 !
2636 !        N IS THE ORDER OF THE MATRIX,
2637 !
2638 !        NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2639 !          AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2640 !
2641 !        A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
2642 !          INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
2643 !          ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
2644 !
2645 !     ON OUTPUT-
2646 !
2647 !        A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
2648 !          TRANSFORMATIONS USED IN THE REDUCTION,
2649 !
2650 !        D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
2651 !
2652 !        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
2653 !          MATRIX IN ITS LAST N-1 POSITIONS.  E(1) IS SET TO ZERO,
2654 !
2655 !        E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
2656 !          E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
2657 !
2658 !     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2659 !     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2660 !
2661 !     ------------------------------------------------------------------
2662 !
2663 !     ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
2664       DO 300 II = 1, N
2665       I = N + 1 - II
2666       L = I - 1
2667       IZ = (I * L) / 2
2668       H = 0.0D+00
2669       SCALE = 0.0D+00
2670       IF (L .LT. 1) GO TO 120
2671 !     ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
2672       DO 100 K = 1, L
2673       IZ = IZ + 1
2674       D(K) = A(IZ)
2675       SCALE = SCALE + ABS(D(K))
2676   100 CONTINUE
2677 !
2678       IF (SCALE .NE. 0.0D+00) GO TO 140
2679   120 E(I) = 0.0D+00
2680       E2(I) = 0.0D+00
2681       GO TO 280
2682 !
2683   140 DO 160 K = 1, L
2684       D(K) = D(K) / SCALE
2685       H = H + D(K) * D(K)
2686   160 CONTINUE
2687 !
2688       E2(I) = SCALE * SCALE * H
2689       F = D(L)
2690       G = -SIGN(SQRT(H),F)
2691       E(I) = SCALE * G
2692       H = H - F * G
2693       D(L) = F - G
2694       A(IZ) = SCALE * D(L)
2695       IF (L .EQ. 1) GO TO 280
2696       F = 0.0D+00
2697 !
2698       JK = 1
2699       DO 220 J = 1, L
2700       JM1 = J - 1
2701       DT = D(J)
2702       G = 0.0D+00
2703 !     ********** FORM ELEMENT OF A*U **********
2704       IF (JM1 .EQ. 0) GO TO 200
2705       DO 180 K = 1, JM1
2706       E(K) = E(K) + DT * A(JK)
2707       G = G + D(K) * A(JK)
2708       JK = JK + 1
2709   180 CONTINUE
2710   200 E(J) = G + A(JK) * DT
2711       JK = JK + 1
2712 !     ********** FORM ELEMENT OF P **********
2713   220 CONTINUE
2714       F = 0.0D+00
2715       DO 240 J = 1, L
2716       E(J) = E(J) / H
2717       F = F + E(J) * D(J)
2718   240 CONTINUE
2719 !
2720       HH = F / (H + H)
2721       JK = 0
2722 !     ********** FORM REDUCED A **********
2723       DO 260 J = 1, L
2724       F = D(J)
2725       G = E(J) - HH * F
2726       E(J) = G
2727 !
2728       DO 260 K = 1, J
2729       JK = JK + 1
2730       A(JK) = A(JK) - F * E(K) - G * D(K)
2731   260 CONTINUE
2732 !
2733   280 D(I) = A(IZ+1)
2734       A(IZ+1) = SCALE * SQRT(H)
2735   300 CONTINUE
2736 !
2737       return
2738 !     ********** LAST CARD OF TRED3 **********
2739       end subroutine TRED3B
2740 !-----------------------------------------------------------------------------
2741 ! blas.f
2742 !-----------------------------------------------------------------------------
2743 ! 10 NOV 94 - MWS - DNRM2: REMOVE FTNCHECK WARNINGS
2744 ! 11 JUN 94 - MWS - INCLUDE A COPY OF DGEMV (LEVEL TWO ROUTINE)
2745 ! 11 AUG 87 - MWS - SANITIZE FLOATING POINT CONSTANTS IN DNRM2
2746 ! 26 MAR 87 - MWS - USE GENERIC SIGN IN DROTG
2747 ! 28 NOV 86 - STE - SUPPLY ALL LEVEL ONE BLAS
2748 !  7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
2749 !
2750 ! BASIC LINEAR ALGEBRA SUBPROGRAMS (BLAS) FROM LINPACK  (LEVEL 1)
2751 !
2752 !   THIS MODULE SHOULD BE COMPILED ONLY IF SPECIALLY CODED
2753 !   VERSIONS OF THESE ROUTINES ARE NOT AVAILABLE ON THE TARGET MACHINE
2754 !
2755 !*MODULE BLAS1   *DECK DASUM
2756       real(kind=8) function DASUM(N,DX,INCX)
2757 !
2758 !     TAKES THE SUM OF THE ABSOLUTE VALUES.
2759 !     JACK DONGARRA, LINPACK, 3/11/78.
2760 !
2761       real(kind=8) :: DX(1),DTEMP
2762       INTEGER :: I,INCX,M,MP1,N,NINCX
2763 !
2764       DASUM = 0.0D+00
2765       DTEMP = 0.0D+00
2766       IF(N.LE.0)RETURN
2767       IF(INCX.EQ.1)GO TO 20
2768 !
2769 !        CODE FOR INCREMENT NOT EQUAL TO 1
2770 !
2771       NINCX = N*INCX
2772       DO 10 I = 1,NINCX,INCX
2773         DTEMP = DTEMP + ABS(DX(I))
2774    10 CONTINUE
2775       DASUM = DTEMP
2776       RETURN
2777 !
2778 !        CODE FOR INCREMENT EQUAL TO 1
2779 !
2780 !
2781 !        CLEAN-UP LOOP
2782 !
2783    20 M = MOD(N,6)
2784       IF( M .EQ. 0 ) GO TO 40
2785       DO 30 I = 1,M
2786         DTEMP = DTEMP + ABS(DX(I))
2787    30 CONTINUE
2788       IF( N .LT. 6 ) GO TO 60
2789    40 MP1 = M + 1
2790       DO 50 I = MP1,N,6
2791         DTEMP = DTEMP + ABS(DX(I)) + ABS(DX(I + 1)) + ABS(DX(I + 2)) &
2792         + ABS(DX(I + 3)) + ABS(DX(I + 4)) + ABS(DX(I + 5))
2793    50 CONTINUE
2794    60 DASUM = DTEMP
2795       return
2796       end function DASUM
2797 !-----------------------------------------------------------------------------
2798 !*MODULE BLAS1   *DECK DAXPY
2799       subroutine DAXPY(N,DA,DX,INCX,DY,INCY)
2800
2801 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2802       integer :: N,INCX,INCY
2803       real(kind=8),DIMENSION(1) :: DX,DY
2804       real(kind=8) :: DA
2805 !
2806 !     CONSTANT TIMES A VECTOR PLUS A VECTOR.
2807 !           DY(I) = DY(I) + DA * DX(I)
2808 !     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
2809 !     JACK DONGARRA, LINPACK, 3/11/78.
2810 !
2811       integer :: ix,iy,i,m,mp1
2812       IF(N.LE.0)RETURN
2813       IF (DA .EQ. 0.0D+00) RETURN
2814       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2815 !
2816 !        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2817 !          NOT EQUAL TO 1
2818 !
2819       IX = 1
2820       IY = 1
2821       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2822       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2823       DO 10 I = 1,N
2824         DY(IY) = DY(IY) + DA*DX(IX)
2825         IX = IX + INCX
2826         IY = IY + INCY
2827    10 CONTINUE
2828       RETURN
2829 !
2830 !        CODE FOR BOTH INCREMENTS EQUAL TO 1
2831 !
2832 !
2833 !        CLEAN-UP LOOP
2834 !
2835    20 M = MOD(N,4)
2836       IF( M .EQ. 0 ) GO TO 40
2837       DO 30 I = 1,M
2838         DY(I) = DY(I) + DA*DX(I)
2839    30 CONTINUE
2840       IF( N .LT. 4 ) RETURN
2841    40 MP1 = M + 1
2842       DO 50 I = MP1,N,4
2843         DY(I) = DY(I) + DA*DX(I)
2844         DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
2845         DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
2846         DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
2847    50 CONTINUE
2848       return
2849       end subroutine DAXPY
2850 !-----------------------------------------------------------------------------
2851 !*MODULE BLAS1   *DECK DCOPY
2852       subroutine DCOPY(N,DX,INCX,DY,INCY)
2853
2854 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2855       integer :: N,INCX,INCY
2856       real(kind=8),DIMENSION(*) :: DX,DY
2857 !
2858 !     COPIES A VECTOR.
2859 !           DY(I) <== DX(I)
2860 !     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
2861 !     JACK DONGARRA, LINPACK, 3/11/78.
2862 !
2863       integer :: ix,iy,m,i,mp1
2864       IF(N.LE.0)RETURN
2865       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2866 !
2867 !        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2868 !          NOT EQUAL TO 1
2869 !
2870       IX = 1
2871       IY = 1
2872       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2873       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2874       DO 10 I = 1,N
2875         DY(IY) = DX(IX)
2876         IX = IX + INCX
2877         IY = IY + INCY
2878    10 CONTINUE
2879       RETURN
2880 !
2881 !        CODE FOR BOTH INCREMENTS EQUAL TO 1
2882 !
2883 !
2884 !        CLEAN-UP LOOP
2885 !
2886    20 M = MOD(N,7)
2887       IF( M .EQ. 0 ) GO TO 40
2888       DO 30 I = 1,M
2889         DY(I) = DX(I)
2890    30 CONTINUE
2891       IF( N .LT. 7 ) RETURN
2892    40 MP1 = M + 1
2893       DO 50 I = MP1,N,7
2894         DY(I) = DX(I)
2895         DY(I + 1) = DX(I + 1)
2896         DY(I + 2) = DX(I + 2)
2897         DY(I + 3) = DX(I + 3)
2898         DY(I + 4) = DX(I + 4)
2899         DY(I + 5) = DX(I + 5)
2900         DY(I + 6) = DX(I + 6)
2901    50 CONTINUE
2902       return
2903       end subroutine DCOPY
2904 !-----------------------------------------------------------------------------
2905 !*MODULE BLAS1   *DECK DDOT
2906       real(kind=8) function DDOT(N,DX,INCX,DY,INCY)
2907
2908 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2909       integer :: N,INCX,INCY
2910       real(kind=8),DIMENSION(1) :: DX,DY
2911 !
2912 !     FORMS THE DOT PRODUCT OF TWO VECTORS.
2913 !           DOT = DX(I) * DY(I)
2914 !     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
2915 !     JACK DONGARRA, LINPACK, 3/11/78.
2916 !
2917       integer ::ix,iy,m,mp1,i
2918       real(kind=8) :: DTEMP
2919       DDOT = 0.0D+00
2920       DTEMP = 0.0D+00
2921       IF(N.LE.0)RETURN
2922       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2923 !
2924 !        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2925 !          NOT EQUAL TO 1
2926 !
2927       IX = 1
2928       IY = 1
2929       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2930       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2931       DO 10 I = 1,N
2932         DTEMP = DTEMP + DX(IX)*DY(IY)
2933         IX = IX + INCX
2934         IY = IY + INCY
2935    10 CONTINUE
2936       DDOT = DTEMP
2937       RETURN
2938 !
2939 !        CODE FOR BOTH INCREMENTS EQUAL TO 1
2940 !
2941 !
2942 !        CLEAN-UP LOOP
2943 !
2944    20 M = MOD(N,5)
2945       IF( M .EQ. 0 ) GO TO 40
2946       DO 30 I = 1,M
2947         DTEMP = DTEMP + DX(I)*DY(I)
2948    30 CONTINUE
2949       IF( N .LT. 5 ) GO TO 60
2950    40 MP1 = M + 1
2951       DO 50 I = MP1,N,5
2952         DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + &
2953          DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
2954    50 CONTINUE
2955    60 DDOT = DTEMP
2956       return
2957       end function DDOT
2958 !-----------------------------------------------------------------------------
2959 !*MODULE BLAS1   *DECK DNRM2
2960       real(kind=8) function DNRM2(N,DX,INCX)
2961
2962       INTEGER      :: NEXT,N,INCX
2963       real(kind=8) :: DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
2964       DATA   ZERO, ONE /0.0D+00, 1.0D+00/
2965
2966       integer :: i,j,nn
2967 !
2968 !     EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
2969 !     INCREMENT INCX .
2970 !     IF    N .LE. 0 RETURN WITH RESULT = 0.
2971 !     IF N .GE. 1 THEN INCX MUST BE .GE. 1
2972 !
2973 !           C.L.LAWSON, 1978 JAN 08
2974 !
2975 !     FOUR PHASE METHOD     USING TWO BUILT-IN CONSTANTS THAT ARE
2976 !     HOPEFULLY APPLICABLE TO ALL MACHINES.
2977 !         CUTLO = MAXIMUM OF  SQRT(U/EPS)  OVER ALL KNOWN MACHINES.
2978 !         CUTHI = MINIMUM OF  SQRT(V)      OVER ALL KNOWN MACHINES.
2979 !     WHERE
2980 !         EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
2981 !         U   = SMALLEST POSITIVE NO.   (UNDERFLOW LIMIT)
2982 !         V   = LARGEST  NO.            (OVERFLOW  LIMIT)
2983 !
2984 !     BRIEF OUTLINE OF ALGORITHM..
2985 !
2986 !     PHASE 1    SCANS ZERO COMPONENTS.
2987 !     MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
2988 !     MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
2989 !     MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
2990 !     WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
2991 !
2992 !     VALUES FOR CUTLO AND CUTHI..
2993 !     FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
2994 !     DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
2995 !     CUTLO, S.P.   U/EPS = 2**(-102) FOR  HONEYWELL.  CLOSE SECONDS ARE
2996 !                   UNIVAC AND DEC AT 2**(-103)
2997 !                   THUS CUTLO = 2**(-51) = 4.44089E-16
2998 !     CUTHI, S.P.   V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
2999 !                   THUS CUTHI = 2**(63.5) = 1.30438E19
3000 !     CUTLO, D.P.   U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
3001 !                   THUS CUTLO = 2**(-33.5) = 8.23181D-11
3002 !     CUTHI, D.P.   SAME AS S.P.  CUTHI = 1.30438D+19
3003 !     DATA CUTLO, CUTHI / 8.232D-11,  1.304D+19 /
3004 !     DATA CUTLO, CUTHI / 4.441E-16,  1.304E19 /
3005       DATA CUTLO, CUTHI / 8.232D-11,  1.304D+19 /
3006 !
3007       J=0
3008       IF(N .GT. 0) GO TO 10
3009          DNRM2  = ZERO
3010          GO TO 300
3011 !
3012    10 ASSIGN 30 TO NEXT
3013       SUM = ZERO
3014       NN = N * INCX
3015 !                                                 BEGIN MAIN LOOP
3016       I = 1
3017    20    GO TO NEXT,(30, 50, 70, 110)
3018    30 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
3019       ASSIGN 50 TO NEXT
3020       XMAX = ZERO
3021 !
3022 !                        PHASE 1.  SUM IS ZERO
3023 !
3024    50 IF( DX(I) .EQ. ZERO) GO TO 200
3025       IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
3026 !
3027 !                                PREPARE FOR PHASE 2.
3028       ASSIGN 70 TO NEXT
3029       GO TO 105
3030 !
3031 !                                PREPARE FOR PHASE 4.
3032 !
3033   100 I = J
3034       ASSIGN 110 TO NEXT
3035       SUM = (SUM / DX(I)) / DX(I)
3036   105 XMAX = ABS(DX(I))
3037       GO TO 115
3038 !
3039 !                   PHASE 2.  SUM IS SMALL.
3040 !                             SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
3041 !
3042    70 IF( ABS(DX(I)) .GT. CUTLO ) GO TO 75
3043 !
3044 !                     COMMON CODE FOR PHASES 2 AND 4.
3045 !                     IN PHASE 4 SUM IS LARGE.  SCALE TO AVOID OVERFLOW.
3046 !
3047   110 IF( ABS(DX(I)) .LE. XMAX ) GO TO 115
3048          SUM = ONE + SUM * (XMAX / DX(I))**2
3049          XMAX = ABS(DX(I))
3050          GO TO 200
3051 !
3052   115 SUM = SUM + (DX(I)/XMAX)**2
3053       GO TO 200
3054 !
3055 !
3056 !                  PREPARE FOR PHASE 3.
3057 !
3058    75 SUM = (SUM * XMAX) * XMAX
3059 !
3060 !
3061 !     FOR REAL OR D.P. SET HITEST = CUTHI/N
3062 !     FOR COMPLEX      SET HITEST = CUTHI/(2*N)
3063 !
3064    85 HITEST = CUTHI/N
3065 !
3066 !                   PHASE 3.  SUM IS MID-RANGE.  NO SCALING.
3067 !
3068       DO 95 J =I,NN,INCX
3069       IF(ABS(DX(J)) .GE. HITEST) GO TO 100
3070    95    SUM = SUM + DX(J)**2
3071       DNRM2 = SQRT( SUM )
3072       GO TO 300
3073 !
3074   200 CONTINUE
3075       I = I + INCX
3076       IF ( I .LE. NN ) GO TO 20
3077 !
3078 !              END OF MAIN LOOP.
3079 !
3080 !              COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
3081 !
3082       DNRM2 = XMAX * SQRT(SUM)
3083   300 CONTINUE
3084       return
3085       end function DNRM2
3086 !-----------------------------------------------------------------------------
3087 !*MODULE BLAS1   *DECK DROT
3088       subroutine DROT(N,DX,INCX,DY,INCY,C,S)
3089
3090 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3091       integer :: N,INCX,INCY
3092       real(kind=8),DIMENSION(1) :: DX,DY
3093       real(kind=8) :: C,S
3094 !
3095 !     APPLIES A PLANE ROTATION.
3096 !           DX(I) =  C*DX(I) + S*DY(I)
3097 !           DY(I) = -S*DX(I) + C*DY(I)
3098 !     JACK DONGARRA, LINPACK, 3/11/78.
3099 !
3100       integer :: ix,iy,i
3101       real(kind=8) :: DTEMP
3102       IF(N.LE.0)RETURN
3103       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
3104 !
3105 !       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
3106 !         TO 1
3107 !
3108       IX = 1
3109       IY = 1
3110       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
3111       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
3112       DO 10 I = 1,N
3113         DTEMP = C*DX(IX) + S*DY(IY)
3114         DY(IY) = C*DY(IY) - S*DX(IX)
3115         DX(IX) = DTEMP
3116         IX = IX + INCX
3117         IY = IY + INCY
3118    10 CONTINUE
3119       RETURN
3120 !
3121 !       CODE FOR BOTH INCREMENTS EQUAL TO 1
3122 !
3123    20 DO 30 I = 1,N
3124         DTEMP = C*DX(I) + S*DY(I)
3125         DY(I) = C*DY(I) - S*DX(I)
3126         DX(I) = DTEMP
3127    30 CONTINUE
3128       return
3129       end subroutine DROT
3130 !-----------------------------------------------------------------------------
3131 !*MODULE BLAS1   *DECK DROTG
3132       subroutine DROTG(DA,DB,C,S)
3133 !
3134 !     CONSTRUCT GIVENS PLANE ROTATION.
3135 !     JACK DONGARRA, LINPACK, 3/11/78.
3136 !
3137       real(kind=8) :: DA,DB,C,S,ROE,SCALE,R,Z
3138 !
3139       real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00
3140 !
3141 !-----------------------------------------------------------------------
3142 !
3143 !
3144       ROE = DB
3145       IF( ABS(DA) .GT. ABS(DB) ) ROE = DA
3146       SCALE = ABS(DA) + ABS(DB)
3147       IF( SCALE .NE. ZERO ) GO TO 10
3148          C = ONE
3149          S = ZERO
3150          R = ZERO
3151          GO TO 20
3152 !
3153    10 R = SCALE*SQRT((DA/SCALE)**2 + (DB/SCALE)**2)
3154       R = SIGN(ONE,ROE)*R
3155       C = DA/R
3156       S = DB/R
3157    20 Z = ONE
3158       IF( ABS(DA) .GT. ABS(DB) ) Z = S
3159       IF( ABS(DB) .GE. ABS(DA) .AND. C .NE. ZERO ) Z = ONE/C
3160       DA = R
3161       DB = Z
3162       return
3163       end subroutine DROTG
3164 !-----------------------------------------------------------------------------
3165 !*MODULE BLAS1   *DECK DSCAL
3166       subroutine DSCAL(N,DA,DX,INCX)
3167
3168 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3169       integer :: N,INCX
3170       real(kind=8),DIMENSION(1) :: DX
3171       real(kind=8) :: DA
3172 !
3173 !     SCALES A VECTOR BY A CONSTANT.
3174 !           DX(I) = DA * DX(I)
3175 !     USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
3176 !     JACK DONGARRA, LINPACK, 3/11/78.
3177 !
3178       integer :: NINCX,m,mp1,i
3179       IF(N.LE.0)RETURN
3180       IF(INCX.EQ.1)GO TO 20
3181 !
3182 !        CODE FOR INCREMENT NOT EQUAL TO 1
3183 !
3184       NINCX = N*INCX
3185       DO 10 I = 1,NINCX,INCX
3186         DX(I) = DA*DX(I)
3187    10 CONTINUE
3188       RETURN
3189 !
3190 !        CODE FOR INCREMENT EQUAL TO 1
3191 !
3192 !
3193 !        CLEAN-UP LOOP
3194 !
3195    20 M = MOD(N,5)
3196       IF( M .EQ. 0 ) GO TO 40
3197       DO 30 I = 1,M
3198         DX(I) = DA*DX(I)
3199    30 CONTINUE
3200       IF( N .LT. 5 ) RETURN
3201    40 MP1 = M + 1
3202       DO 50 I = MP1,N,5
3203         DX(I) = DA*DX(I)
3204         DX(I + 1) = DA*DX(I + 1)
3205         DX(I + 2) = DA*DX(I + 2)
3206         DX(I + 3) = DA*DX(I + 3)
3207         DX(I + 4) = DA*DX(I + 4)
3208    50 CONTINUE
3209       return
3210       end subroutine DSCAL
3211 !-----------------------------------------------------------------------------
3212 !*MODULE BLAS1   *DECK DSWAP
3213       subroutine DSWAP(N,DX,INCX,DY,INCY)
3214
3215 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3216       integer :: N,INCX,INCY
3217       real(kind=8),DIMENSION(1) :: DX,DY
3218 !
3219 !     INTERCHANGES TWO VECTORS.
3220 !           DX(I) <==> DY(I)
3221 !     USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
3222 !     JACK DONGARRA, LINPACK, 3/11/78.
3223 !
3224       integer :: ix,iy,i,m,mp1
3225       real(kind=8) :: DTEMP
3226       IF(N.LE.0)RETURN
3227       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
3228 !
3229 !       CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
3230 !         TO 1
3231 !
3232       IX = 1
3233       IY = 1
3234       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
3235       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
3236       DO 10 I = 1,N
3237         DTEMP = DX(IX)
3238         DX(IX) = DY(IY)
3239         DY(IY) = DTEMP
3240         IX = IX + INCX
3241         IY = IY + INCY
3242    10 CONTINUE
3243       RETURN
3244 !
3245 !       CODE FOR BOTH INCREMENTS EQUAL TO 1
3246 !
3247 !
3248 !       CLEAN-UP LOOP
3249 !
3250    20 M = MOD(N,3)
3251       IF( M .EQ. 0 ) GO TO 40
3252       DO 30 I = 1,M
3253         DTEMP = DX(I)
3254         DX(I) = DY(I)
3255         DY(I) = DTEMP
3256    30 CONTINUE
3257       IF( N .LT. 3 ) RETURN
3258    40 MP1 = M + 1
3259       DO 50 I = MP1,N,3
3260         DTEMP = DX(I)
3261         DX(I) = DY(I)
3262         DY(I) = DTEMP
3263         DTEMP = DX(I + 1)
3264         DX(I + 1) = DY(I + 1)
3265         DY(I + 1) = DTEMP
3266         DTEMP = DX(I + 2)
3267         DX(I + 2) = DY(I + 2)
3268         DY(I + 2) = DTEMP
3269    50 CONTINUE
3270       return
3271       end subroutine DSWAP
3272 !-----------------------------------------------------------------------------
3273 !*MODULE BLAS1   *DECK IDAMAX
3274       integer function IDAMAX(N,DX,INCX)
3275
3276 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3277       integer :: N,INCX
3278       real(kind=8),DIMENSION(1) :: DX
3279 !
3280 !     FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
3281 !     JACK DONGARRA, LINPACK, 3/11/78.
3282 !
3283       integer :: ix,iy,i
3284       real(kind=8) :: RMAX
3285       IDAMAX = 0
3286       IF( N .LT. 1 ) RETURN
3287       IDAMAX = 1
3288       IF(N.EQ.1)RETURN
3289       IF(INCX.EQ.1)GO TO 20
3290 !
3291 !        CODE FOR INCREMENT NOT EQUAL TO 1
3292 !
3293       IX = 1
3294       RMAX = ABS(DX(1))
3295       IX = IX + INCX
3296       DO 10 I = 2,N
3297          IF(ABS(DX(IX)).LE.RMAX) GO TO 5
3298          IDAMAX = I
3299          RMAX = ABS(DX(IX))
3300     5    IX = IX + INCX
3301    10 CONTINUE
3302       RETURN
3303 !
3304 !        CODE FOR INCREMENT EQUAL TO 1
3305 !
3306    20 RMAX = ABS(DX(1))
3307       DO 30 I = 2,N
3308          IF(ABS(DX(I)).LE.RMAX) GO TO 30
3309          IDAMAX = I
3310          RMAX = ABS(DX(I))
3311    30 CONTINUE
3312       return
3313       end function IDAMAX
3314 !-----------------------------------------------------------------------------
3315 !*MODULE BLAS    *DECK DGEMV
3316       subroutine DGEMV(FORMA,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
3317
3318 !      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3319       integer :: M,N,INCX,INCY,LDA
3320       CHARACTER(len=1) :: FORMA
3321       real(kind=8),DIMENSION(LDA,*) :: A
3322       real(kind=8),DIMENSION(*) :: X,Y
3323       real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00
3324       real(kind=8) :: ALPHA,BETA
3325       integer :: i,locy
3326 !
3327 !        CLONE OF -DGEMV- WRITTEN BY MIKE SCHMIDT
3328 !
3329       LOCY = 1
3330       IF(FORMA.EQ.'T') GO TO 200
3331 !
3332 !                  Y = ALPHA * A * X + BETA * Y
3333 !
3334       IF(ALPHA.EQ.ONE  .AND.  BETA.EQ.ZERO) THEN
3335          DO 110 I=1,M
3336             Y(LOCY) =       DDOT(N,A(I,1),LDA,X,INCX)
3337             LOCY = LOCY+INCY
3338   110    CONTINUE
3339       ELSE
3340          DO 120 I=1,M
3341             Y(LOCY) = ALPHA*DDOT(N,A(I,1),LDA,X,INCX) + BETA*Y(LOCY)
3342             LOCY = LOCY+INCY
3343   120    CONTINUE
3344       END IF
3345       RETURN
3346 !
3347 !                  Y = ALPHA * A-TRANSPOSE * X + BETA * Y
3348 !
3349   200 CONTINUE
3350       IF(ALPHA.EQ.ONE  .AND.  BETA.EQ.ZERO) THEN
3351          DO 210 I=1,N
3352             Y(LOCY) =       DDOT(M,A(1,I),1,X,INCX)
3353             LOCY = LOCY+INCY
3354   210    CONTINUE
3355       ELSE
3356          DO 220 I=1,N
3357             Y(LOCY) = ALPHA*DDOT(M,A(1,I),1,X,INCX) + BETA*Y(LOCY)
3358             LOCY = LOCY+INCY
3359   220    CONTINUE
3360       END IF
3361       return
3362       end subroutine DGEMV
3363 !-----------------------------------------------------------------------------
3364 !-----------------------------------------------------------------------------
3365       end module MD_calc