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