2 !-----------------------------------------------------------------------------
4 use MD_data, only:D_ban,IP
6 ! use prng ! prng.f90 or prng_32.f90
9 !-----------------------------------------------------------------------------
11 !-----------------------------------------------------------------------------
13 !-----------------------------------------------------------------------------
17 !-----------------------------------------------------------------------------
18 !*MODULE MTHLIB *DECK VCLR
19 subroutine VCLR(A,INCA,N)
21 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
23 real(kind=8),DIMENSION(*) :: A
25 real(kind=8),PARAMETER :: ZERO=0.0D+00
29 ! ----- ZERO OUT VECTOR -A-, USING INCREMENT -INCA- -----
31 IF (INCA .NE. 1) GO TO 200
45 !-----------------------------------------------------------------------------
47 !-----------------------------------------------------------------------------
48 subroutine BANACH(N,NMAX,A,X,osob)
49 !**********************
51 ! implicit real*8 (a-h,o-z)
52 ! include 'DIMENSIONS'
54 real(kind=8),DIMENSION(NMAX,NMAX) :: A
55 real(kind=8),DIMENSION(NMAX) :: X
56 !el real(kind=8),DIMENSION(6*nres) :: D !(MAXRES6) maxres6=6*maxres
59 real(kind=8) :: xx,aij,aijd
62 !el allocate(D_ban(6*nres))
65 if (dabs(a(1,1)).lt.1.0d-15) then
87 if (dabs(xx).lt.1.0d-15) then
94 CALL BANAII(N,NMAX,A,X)
97 !-----------------------------------------------------------------------------
98 subroutine BANAII(N,NMAX,A,X)
99 !************************
100 ! implicit real*8 (a-h,o-z)
101 ! include 'DIMENSIONS'
103 real(kind=8),DIMENSION(NMAX,NMAX) :: A
104 real(kind=8),DIMENSION(NMAX) :: X
105 !el real(kind=8),DIMENSION(6*nres) :: D !(MAXRES6) maxres6=6*maxres
106 !el COMMON /BANII/ D ---> D_ban
126 end subroutine BANAII
127 !-----------------------------------------------------------------------------
128 subroutine MATINVERT(N,NMAX,A,A1,osob)
130 ! implicit real*8 (a-h,o-z)
131 ! include 'DIMENSIONS'
133 real(kind=8),DIMENSION(NMAX,NMAX) :: A,A1
134 !el real(kind=8),DIMENSION(6*nres) :: D !(MAXRES6) maxres6=6*maxres
136 real(kind=8),DIMENSION(NMAX) :: X
143 CALL BANACH(N,NMAX,A,X,osob)
153 CALL BANAII(N,NMAX,A,X)
159 end subroutine MATINVERT
160 !-----------------------------------------------------------------------------
162 !-----------------------------------------------------------------------------
163 subroutine bond_move(nbond,nstart,psi,lprint,error)
165 use mcm_data, only:print_mc
166 use geometry, only:alpha,beta,refsys,matmult
167 ! Move NBOND fragment starting from the CA(nstart) by angle PSI.
168 ! implicit real*8 (a-h,o-z)
169 ! include 'DIMENSIONS'
170 integer :: nbond,nstart
172 logical :: fail,error,lprint
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.CHAIN'
175 ! include 'COMMON.VAR'
176 ! include 'COMMON.IOUNITS'
177 ! include 'COMMON.MCM'
178 real(kind=8),dimension(3) :: x,e1,e2,e3
179 real(kind=8),dimension(3,3) :: e,rot,trans
180 real(kind=8) :: cospsi,sinpsi,rij
181 integer :: i,j,nend,i2,i3,i4,k
184 if (print_mc.gt.2) then
185 write (iout,*) 'nstart=',nstart,' nend=',nend,' nbond=',nbond
186 write (iout,*) 'psi=',psi
187 write (iout,'(a)') 'Original coordinates of the fragment'
189 write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
192 if (nstart.lt.1 .or. nend .gt.nres .or. nbond.lt.2 .or. &
193 nbond.ge.nres-1) then
194 write (iout,'(a)') 'Bad data in BOND_MOVE.'
198 ! Generate the reference system.
202 call refsys(i2,i3,i4,e1,e2,e3,error)
203 ! Return, if couldn't define the reference system.
205 ! Compute the transformation matrix.
223 if (print_mc.gt.2) then
224 write (iout,'(a)') 'Reference system and matrix r:'
226 write(iout,'(i5,2(3f10.5,5x))')i,(e(i,j),j=1,3),(rot(i,j),j=1,3)
230 call matmult(rot,e,trans)
238 call matmult(e,trans,trans)
241 write (iout,'(a)') 'The trans matrix:'
243 write (iout,'(i5,3f10.5)') i,(trans(i,j),j=1,3)
251 rij=rij+trans(j,k)*(c(k,i)-c(k,nstart))
261 write (iout,'(a)') 'Rotated coordinates of the fragment'
263 write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
267 ! call int_from_cart(.false.,lprint)
268 if (nstart.gt.1) then
269 theta(nstart+1)=alpha(nstart-1,nstart,nstart+1)
270 phi(nstart+2)=beta(nstart-1,nstart,nstart+1,nstart+2)
271 if (nstart.gt.2) phi(nstart+1)= &
272 beta(nstart-2,nstart-1,nstart,nstart+1)
274 if (nend.lt.nres) then
275 theta(nend+1)=alpha(nend-1,nend,nend+1)
276 phi(nend+1)=beta(nend-2,nend-1,nend,nend+1)
277 if (nend.lt.nres-1) phi(nend+2)= &
278 beta(nend-1,nend,nend+1,nend+2)
280 if (print_mc.gt.2) then
281 write (iout,'(/a,i3,a,i3,a/)') &
282 'Moved internal coordinates of the ',nstart,'-',nend,&
284 do i=nstart+1,nstart+2
285 write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
288 write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
292 end subroutine bond_move
293 !-----------------------------------------------------------------------------
295 !-----------------------------------------------------------------------------
296 ! 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
297 ! 31 MAR 94 - MWS - ADD A VARIABLE TO END OF MACHSW COMMON
298 ! 26 JUN 93 - MWS - ETRED3: ADD RETURN FOR SPECIAL CASE N=1
299 ! 4 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
300 ! 30 AUG 91 - MWS - JACDIA: LIMIT ITERATIONS, USE EPSLON IN TEST.
301 ! 14 JUL 91 - MWS - JACOBI DIAGONALIZATION ALLOWS FOR LDVEC.NE.N
302 ! 29 JAN 91 - TLW - GLDIAG: CHANGED COMMON DIAGSW TO MACHSW
303 ! 29 OCT 90 - STE - FIX JACDIA UNDEFINED VARIABLE BUG
304 ! 14 SEP 90 - MK - NEW JACOBI DIAGONALIZATION (KDIAG=3)
305 ! 27 MAR 88 - MWS - ALLOW FOR VECTOR ROUTINE IN GLDIAG
306 ! 11 AUG 87 - MWS - SANITIZE CONSTANTS IN EQLRAT
307 ! 15 FEB 87 - STE - FIX EINVIT SUB-MATRIX LOOP LIMIT
308 ! SCRATCH ARRAYS ARE N*8 REAL AND N INTEGER
309 ! 8 DEC 86 - STE - USE PERF INDEX FROM ESTPI1 TO JUDGE EINVIT FAILURE
310 ! 30 NOV 86 - STE - DELETE LIGENB, MAKE EVVRSP DEFAULT
311 ! (GIVEIS FAILS ON CRAY FOR BENCHMC AND BENCHCI)
312 ! 7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
313 ! 11 OCT 85 - STE - LIGENB,TQL2: USE DROT,DSWAP; TINVTB: SCALE VECTOR
314 ! BEFORE NORMALIZING; GENERIC FUNCTIONS
315 ! 24 FEB 84 - STE - INITIALIZE INDEX ARRAY FOR LIGENB IN GLDIAG
316 ! 1 DEC 83 - STE - CHANGE MACHEP FROM 2**-54 TO 2**-50
317 ! 28 SEP 82 - MWS - CONVERT TO IBM
319 !*MODULE EIGEN *DECK EINVIT
320 subroutine EINVIT(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
323 !* THIS IS A MODIFICATION OF TINVIT FROM EISPACK EDITION 3
324 !* DATED AUGUST 1983.
325 !* TINVIT IS A TRANSLATION OF THE INVERSE ITERATION TECHNIQUE
326 !* IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
327 !* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
328 !* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
331 !* THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
332 !* SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES.
335 !* INVERSE ITERATION.
339 !* MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
340 !* ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
341 !* DIMENSION STATEMENT.
344 !* CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
346 !* CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
347 !* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY.
348 !* E2 - W.P. REAL (N)
349 !* CONTAINS THE SQUARES OF CORRESPONDING ELEMENTS OF E,
350 !* WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
351 !* E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
352 !* THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE
353 !* SUM OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST
354 !* CONTAIN 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER,
355 !* OR 2.0 IF THE EIGENVALUES ARE IN DESCENDING ORDER.
356 !* IF TQLRAT, BISECT, TRIDIB, OR IMTQLV
357 !* HAS BEEN USED TO FIND THE EIGENVALUES, THEIR
358 !* OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE.
360 !* THE NUMBER OF SPECIFIED EIGENVECTORS.
362 !* CONTAINS THE M EIGENVALUES IN ASCENDING
363 !* OR DESCENDING ORDER.
365 !* CONTAINS IN FIRST M POSITIONS THE SUBMATRIX INDICES
366 !* ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
367 !* 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX
368 !* FROM THE TOP, 2 FOR THOSE BELONGING TO THE SECOND
370 !* IERR - INTEGER (LOGICAL UNIT NUMBER)
371 !* LOGICAL UNIT FOR ERROR MESSAGES
374 !* ALL INPUT ARRAYS ARE UNALTERED.
375 !* Z - W.P. REAL (NM,M)
376 !* CONTAINS THE ASSOCIATED SET OF ORTHONORMAL
377 !* EIGENVECTORS. ANY VECTOR WHICH WHICH FAILS TO CONVERGE
378 !* IS LEFT AS IS (BUT NORMALIZED) WHEN ITERATING STOPPED.
381 !* ZERO FOR NORMAL RETURN,
382 !* -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
383 !* EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS.
384 !* (ONLY LAST FAILURE TO CONVERGE IS REPORTED)
386 !* RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
388 !* RV1 - W.P. REAL (N)
389 !* DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
390 !* RV2 - W.P. REAL (N)
391 !* SUPER(1)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
392 !* RV3 - W.P. REAL (N)
393 !* SUPER(2)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
394 !* RV4 - W.P. REAL (N)
395 !* ELEMENTS DEFINING L IN LU DECOMPOSITION
396 !* RV6 - W.P. REAL (N)
397 !* APPROXIMATE EIGENVECTOR
399 !* DIFFERENCES FROM EISPACK 3 -
400 !* EPS3 IS SCALED BY EPSCAL (ENHANCES CONVERGENCE, BUT
402 !* ONE MORE ITERATION (MINIMUM 2) IS PERFORMED AFTER CONVERGENCE
403 !* (ENHANCES ACCURACY)!
404 !* REPLACE LOOP WITH PYTHAG WITH SINGLE CALL TO DNRM2!
405 !* IF NOT CONVERGED, USE PERFORMANCE INDEX TO DECIDE ON ERROR
406 !* VALUE SETTING, BUT DO NOT STOP!
407 !* L.U. FOR ERROR MESSAGES PASSED THROUGH IERR
408 !* USE PARAMETER STATEMENTS AND GENERIC INTRINSIC FUNCTIONS
410 !* USE IF-THEN-ELSE TO CLARIFY LOGIC
411 !* LOOP OVER SUBSPACES MADE INTO DO LOOP.
412 !* LOOP OVER INVERSE ITERATIONS MADE INTO DO LOOP
413 !* ZERO ONLY REQUIRED PORTIONS OF OUTPUT VECTOR
416 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
417 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
421 LOGICAL :: CONVGD !el,GOPARR,DSKWRK,MASWRK
423 INTEGER :: GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG
426 real(kind=8),dimension(N) :: D,E2
427 real(kind=8) :: E(*)!el E(L)
428 real(kind=8) :: W(M),Z(NM,M)
429 real(kind=8),dimension(N) :: RV1,RV2,RV3,RV4,RV6
430 real(kind=8) :: ANORM,EPS2,EPS3,EPS4,NORM,ORDER,RHO,U,UK,V
431 real(kind=8) :: X0,X1,XU
432 ! real(kind=8) :: ESTPI1 !, DASUM, DDOT, DNRM2 EPSLON,
434 !el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
435 !el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
437 real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00, GRPTOL = 0.001D+00
438 real(kind=8),PARAMETER :: EPSCAL = 0.5D+00, HUNDRD = 100.0D+00, TEN = 10.0D+00
440 001 FORMAT(' EIGENVECTOR ROUTINE EINVIT DID NOT CONVERGE FOR VECTOR' &
441 ,I5,'. NORM =',1P,E10.2,' PERFORMANCE INDEX =',E10.2/ &
442 ' (AN ERROR HALT WILL OCCUR IF THE PI IS GREATER THAN 100)')
445 !-----------------------------------------------------------------------
462 ! .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
465 IF (E2(Q+1) .EQ. ZERO) GO TO 140
469 ! .......... FIND VECTORS BY INVERSE ITERATION ..........
477 IF (IND(R) .NE. TAG) GO TO 920
480 IF (S .NE. 0) GO TO 510
482 ! .......... CHECK FOR ISOLATED ROOT ..........
493 NORM = MAX( NORM, ABS(D(I)) + ABS(E(I)) )
496 ! .......... EPS2 IS THE CRITERION FOR GROUPING,
497 ! EPS3 REPLACES ZERO PIVOTS AND EQUAL
498 ! ROOTS ARE MODIFIED BY EPS3,
499 ! EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
502 EPS3 = EPSCAL * EPSLON(NORM)
510 ! .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
512 510 IF (ABS(X1-X0) .GE. EPS2) THEN
522 IF (ORDER * (X1 - X0) .LE. EPS3) X1 = X0 + ORDER * EPS3
525 ! .......... ELIMINATION WITH INTERCHANGES AND
526 ! INITIALIZATION OF VECTOR ..........
535 IF (ABS(E(I)) .GT. ABS(U)) THEN
537 ! EXCHANGE ROWS BEFORE ELIMINATION
539 ! *** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
540 ! E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......
547 U = V - XU * RV2(I-1)
552 ! STRAIGHT ELIMINATION
559 U = D(I) - X1 - XU * V
564 IF (ABS(U) .LE. EPS3) U = EPS3
569 ! DO INVERSE ITERATIONS
573 IF (ITS .EQ. 1) GO TO 600
575 ! .......... FORWARD SUBSTITUTION ..........
577 IF (NORM .EQ. ZERO) THEN
583 CALL DSCAL (Q-P+1, XU, RV6(P), 1)
586 ! ... ELIMINATION OPERATIONS ON NEXT VECTOR
591 ! IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
592 ! WAS PERFORMED EARLIER IN THE
593 ! TRIANGULARIZATION PROCESS ..........
595 IF (RV1(I-1) .EQ. E(I)) THEN
601 RV6(I) = U - RV4(I) * RV6(I-1)
605 ! .......... BACK SUBSTITUTION
607 RV6(Q) = RV6(Q) / RV1(Q)
611 DO 620 I = Q-1, P, -1
612 RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
617 IF (GROUP .EQ. 0) GO TO 700
619 ! ....... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
620 ! MEMBERS OF GROUP ..........
625 IF (IND(J) .NE. TAG) GO TO 630
626 CALL DAXPY(Q-P+1, -DDOT(Q-P+1,RV6(P),1,Z(P,J),1),&
629 NORM = DASUM(Q-P+1, RV6(P), 1)
632 IF (CONVGD) GO TO 840
633 IF (NORM .GE. ONE) CONVGD = .TRUE.
636 ! .......... NORMALIZE SO THAT SUM OF SQUARES IS
637 ! 1 AND EXPAND TO FULL ORDER ..........
641 XU = ONE / DNRM2(Q-P+1,RV6(P),1)
654 IF (.NOT.CONVGD) THEN
655 RHO = ESTPI1(Q-P+1,X1,D(P),E(P),Z(P,R),ANORM)
656 IF (RHO .GE. TEN .AND. LUEMSG .GT. 0 .AND. MASWRK) &
657 WRITE(LUEMSG,001) R,NORM,RHO
659 ! *** SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
661 IF (RHO .GT. HUNDRD) IERR = -R
667 IF (Q .EQ. N) GO TO 940
671 end subroutine EINVIT
672 !-----------------------------------------------------------------------------
673 !*MODULE EIGEN *DECK ELAUM
674 subroutine ELAU(HINV,L,D,A,E)
676 integer :: L,JL,JK,J,JM1,K
679 !el real(kind=8) :: E(L)
680 real(kind=8) :: E(*)!el E(L)
686 real(kind=8),PARAMETER :: ZERO = 0.0D+00, HALF = 0.5D+00
698 E(K) = E(K) + A(JK) * F
706 ! .......... FORM P ..........
714 ! .......... FORM Q ..........
718 250 E(J) = E(J) - HH * D(J)
722 !-----------------------------------------------------------------------------
723 !*MODULE EIGEN *DECK EPSLON
724 real(kind=8) function EPSLON(X)
727 !* THIS ROUTINE WAS TAKEN FROM EISPACK EDITION 3 DATED 4/6/83
728 !* THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE NOV 1986
731 !* ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
734 !* X - WORKING PRECISION REAL
735 !* VALUES TO FIND EPSLON FOR
738 !* EPSLON - WORKING PRECISION REAL
739 !* SMALLEST POSITIVE VALUE SUCH THAT X+EPSLON .NE. ZERO
742 !* THIS ROUTINE SHOULD PERFORM PROPERLY ON ALL SYSTEMS
743 !* SATISFYING THE FOLLOWING TWO ASSUMPTIONS,
744 !* 1. THE BASE USED IN REPRESENTING FLOATING POINT
745 !* NUMBERS IS NOT A POWER OF THREE.
746 !* 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO
747 !* THE ACCURACY USED IN FLOATING POINT VARIABLES
748 !* THAT ARE STORED IN MEMORY.
749 !* THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO
750 !* FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING
752 !* UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT,
753 !* A IS NOT EXACTLY EQUAL TO FOUR-THIRDS,
754 !* B HAS A ZERO FOR ITS LAST BIT OR DIGIT,
755 !* C IS NOT EXACTLY EQUAL TO ONE,
756 !* EPS MEASURES THE SEPARATION OF 1.0 FROM
757 !* THE NEXT LARGER FLOATING POINT NUMBER.
758 !* THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED
759 !* ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD.
761 !* DIFFERENCES FROM EISPACK 3 -
762 !* USE IS MADE OF PARAMETER STATEMENTS AND INTRINSIC FUNCTIONS
763 !* --NO EXECUTEABLE CODE CHANGES--
766 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
767 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
769 real(kind=8) :: A,B,C,EPS,X
771 real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, FOUR=4.0D+00
773 !-----------------------------------------------------------------------
779 IF (EPS .EQ. ZERO) GO TO 10
783 !-----------------------------------------------------------------------------
784 !*MODULE EIGEN *DECK EQLRAT
785 subroutine EQLRAT(N,DIAG,E,E2IN,D,IND,IERR,E2)
788 !* THIS IS A MODIFICATION OF ROUTINE EQLRAT FROM EISPACK EDITION 3
789 !* DATED AUGUST 1983.
790 !* TQLRAT IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
791 !* ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
792 !* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
795 !* THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
796 !* TRIDIAGONAL MATRIX
803 !* THE ORDER OF THE MATRIX.
805 !* CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
806 !* E2 - W.P. REAL (N)
807 !* CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF
808 !* THE INPUT MATRIX IN ITS LAST N-1 POSITIONS.
809 !* E2(1) IS ARBITRARY.
813 !* CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
814 !* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
815 !* ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
816 !* THE SMALLEST EIGENVALUES.
817 !* E2 - W.P. REAL (N)
821 !* ZERO FOR NORMAL RETURN,
822 !* J IF THE J-TH EIGENVALUE HAS NOT BEEN
823 !* DETERMINED AFTER 30 ITERATIONS.
825 !* DIFFERENCES FROM EISPACK 3 -
826 !* G=G+B INSTEAD OF IF(G.EQ.0) G=B ; B=B/4
827 !* F77 BACKWARD LOOPS INSTEAD OF F66 CONSTRUCT
828 !* GENERIC INTRINSIC FUNCTIONS
829 !* ARRARY IND ADDED FOR USE BY EINVIT
832 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
833 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
835 INTEGER :: I,J,L,M,N,II,L1,IERR
836 INTEGER,dimension(N) :: IND
838 real(kind=8),dimension(N) :: D,E,E2,DIAG,E2IN
839 real(kind=8) :: B,C,F,G,H,P,R,S,T !,EPSLON
841 real(kind=8),PARAMETER :: ZERO = 0.0D+00, SCALE= 1.0D+00/64.0D+00, ONE = 1.0D+00
844 !-----------------------------------------------------------------------
850 IF (N .EQ. 1) GO TO 1001
854 100 E2(I-1) = E2IN(I)
864 H = ABS(D(L)) + ABS(E(L))
865 IF (T .GE. H) GO TO 105
871 ! .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
874 IF (E2(M) .GT. C) GO TO 110
875 ! .......... E2(N) IS ALWAYS ZERO, SO THERE IS AN EXIT
876 ! FROM THE LOOP ..........
878 IF (M .LE. K) GO TO 125
879 IF (M .NE. N) E2IN(M+1) = ZERO
883 IF (M .EQ. L) GO TO 210
888 ! .......... FORM SHIFT ..........
892 P = (D(L1) - G) / (2.0D+00 * S)
893 R = SQRT(P*P+1.0D+00)
894 D(L) = S / (P + SIGN(R,P))
901 ! .......... RATIONAL QL TRANSFORMATION ..........
910 D(I+1) = H + S * (H + D(I))
911 G = D(I) - E2(I) / G + B
917 ! .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST
918 IF (H .EQ. ZERO) GO TO 210
919 IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
921 IF (E2(L) .EQ. ZERO) GO TO 210
923 ! .......... SET ERROR -- NO CONVERGENCE TO AN
924 ! EIGENVALUE AFTER 30 ITERATIONS ..........
931 ! .......... ORDER EIGENVALUES ..........
933 IF (L .EQ. 1) GO TO 250
934 IF (P .LT. D(1)) GO TO 230
936 ! .......... LOOP TO FIND ORDERED POSITION
938 IF (P .LT. D(I)) GO TO 220
941 IF (I .EQ. L) GO TO 250
943 DO 240 II = L, I+1, -1
954 end subroutine EQLRAT
955 !-----------------------------------------------------------------------------
956 !*MODULE EIGEN *DECK ESTPI1
957 real(kind=8) function ESTPI1(N,EVAL,D,E,X,ANORM)
960 !* STEPHEN T. ELBERT (AMES LABORATORY-USDOE) DATE: 5 DEC 1986
963 !* EVALUATE SYMMETRIC TRIDIAGONAL MATRIX PERFORMANCE INDEX
969 !* THIS ROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX A*X-X*EVAL
970 !* WHERE A IS A SYMMETRIC TRIDIAGONAL MATRIX STORED
971 !* IN THE DIAGONAL (D) AND SUB-DIAGONAL (E) VECTORS, EVAL IS THE
972 !* EIGENVALUE OF AN EIGENVECTOR OF A, NAMELY X.
973 !* THIS NORM IS SCALED BY MACHINE ACCURACY FOR THE PROBLEM SIZE.
974 !* ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
978 !* THE ORDER OF THE MATRIX A.
980 !* THE EIGENVALUE CORRESPONDING TO VECTOR X.
982 !* THE DIAGONAL VECTOR OF A.
984 !* THE SUB-DIAGONAL VECTOR OF A.
986 !* AN EIGENVECTOR OF A.
988 !* THE NORM OF A IF IT HAS BEEN PREVIOUSLY COMPUTED.
992 !* THE NORM OF A, COMPUTED IF INITIALLY ZERO.
993 !* ESTPI1 - W.P. REAL
994 !* !!A*X-X*EVAL!! / (EPSLON(10*N)*!!A!!*!!X!!);
995 !* WHERE EPSLON(X) IS THE SMALLEST NUMBER SUCH THAT
996 !* X + EPSLON(X) .NE. X
998 !* ESTPI1 .LT. 1 == SATISFACTORY PERFORMANCE
999 !* .GE. 1 AND .LE. 100 == MARGINAL PERFORMANCE
1000 !* .GT. 100 == POOR PERFORMANCE
1001 !* (SEE LECT. NOTES IN COMP. SCI. VOL.6 PP 124-125)
1004 real(kind=8) :: ANORM,EVAL,RNORM,SIZE,XNORM
1005 real(kind=8),dimension(N) :: D,X
1006 real(kind=8) :: E(*)!el E(L)
1007 ! real(kind=8) :: EPSLON
1009 real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00
1011 !-----------------------------------------------------------------------
1014 IF( N .LE. 1 ) RETURN
1016 IF (ANORM .EQ. ZERO) THEN
1020 ANORM = MAX( ABS(D(1)) + ABS(E(2)), &
1021 ABS(D(N)) + ABS(E(N)))
1023 ANORM = MAX( ANORM, ABS(E(I))+ABS(D(I))+ABS(E(I+1)))
1025 IF(ANORM .EQ. ZERO) ANORM = ONE
1028 ! COMPUTE NORMS OF RESIDUAL AND EIGENVECTOR
1030 XNORM = ABS(X(1)) + ABS(X(N))
1031 RNORM = ABS( (D(1)-EVAL)*X(1) + E(2)*X(2)) &
1032 +ABS( (D(N)-EVAL)*X(N) + E(N)*X(N-1))
1034 XNORM = XNORM + ABS(X(I))
1035 RNORM = RNORM + ABS(E(I)*X(I-1) + (D(I)-EVAL)*X(I) &
1039 ESTPI1 = RNORM / (EPSLON(SIZE)*ANORM*XNORM)
1042 !-----------------------------------------------------------------------------
1043 !*MODULE EIGEN *DECK ETRBK3
1044 subroutine ETRBK3(NM,N,NV,A,M,Z)
1047 !* THIS IS A MODIFICATION OF ROUTINE TRBAK3 FROM EISPACK EDITION 3
1048 !* DATED AUGUST 1983.
1049 !* EISPACK TRBAK3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
1050 !* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
1051 !* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
1052 !* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
1055 !* THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
1056 !* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
1057 !* SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY ETRED3.
1060 !* THE CALCULATION IS CARRIED OUT BY FORMING THE MATRIX PRODUCT
1062 !* WHERE Q IS A PRODUCT OF THE ORTHOGONAL SYMMETRIC MATRICES
1063 !* Q = PROD(I)[1 - U(I)*.TRANSPOSE.U(I)*H(I)]
1064 !* U IS THE AUGMENTED SUB-DIAGONAL ROWS OF A AND
1065 !* Z IS THE SET OF EIGENVECTORS OF THE TRIDIAGONAL
1066 !* MATRIX F WHICH WAS FORMED FROM THE ORIGINAL SYMMETRIC
1067 !* MATRIX C BY THE SIMILARITY TRANSFORMATION
1068 !* F = Q(TRANSPOSE) C Q
1069 !* NOTE THAT ETRBK3 PRESERVES VECTOR EUCLIDEAN NORMS.
1077 !* MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1078 !* ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
1079 !* DIMENSION STATEMENT.
1081 !* THE ORDER OF THE MATRIX A.
1083 !* MUST BE SET TO THE DIMENSION OF THE ARRAY A AS
1084 !* DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT.
1085 !* A - W.P. REAL (NV)
1086 !* CONTAINS INFORMATION ABOUT THE ORTHOGONAL
1087 !* TRANSFORMATIONS USED IN THE REDUCTION BY ETRED3 IN
1088 !* ITS FIRST NV = N*(N+1)/2 POSITIONS.
1090 !* THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
1091 !* Z - W.P REAL (NM,M)
1092 !* CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
1093 !* IN ITS FIRST M COLUMNS.
1096 !* Z - W.P. REAL (NM,M)
1097 !* CONTAINS THE TRANSFORMED EIGENVECTORS
1098 !* IN ITS FIRST M COLUMNS.
1100 !* DIFFERENCES WITH EISPACK 3 -
1101 !* THE TWO INNER LOOPS ARE REPLACED BY DDOT AND DAXPY.
1102 !* MULTIPLICATION USED INSTEAD OF DIVISION TO FIND S.
1103 !* OUTER LOOP RANGE CHANGED FROM 2,N TO 3,N.
1104 !* ADDRESS POINTERS FOR A SIMPLIFIED.
1107 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
1108 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
1110 INTEGER :: I,II,IM1,IZ,J,M,N,NM,NV
1112 real(kind=8) :: A(NV),Z(NM,M)
1113 real(kind=8) :: H,S !,DDOT
1115 real(kind=8),PARAMETER :: ZERO = 0.0D+00
1117 !-----------------------------------------------------------------------
1119 IF (M .EQ. 0) RETURN
1120 IF (N .LE. 2) RETURN
1127 IF (H .EQ. ZERO) GO TO 140
1130 S = -( DDOT(IM1,A(IZ),1,Z(1,J),1) * H) * H
1131 CALL DAXPY(IM1,S,A(IZ),1,Z(1,J),1)
1135 end subroutine ETRBK3
1136 !-----------------------------------------------------------------------------
1137 !*MODULE EIGEN *DECK ETRED3
1138 subroutine ETRED3(N,NV,A,D,E,E2)
1141 !* THIS IS A MODIFICATION OF ROUTINE TRED3 FROM EISPACK EDITION 3
1142 !* DATED AUGUST 1983.
1143 !* EISPACK TRED3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
1144 !* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
1145 !* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
1146 !* THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE JUN 1986
1149 !* THIS ROUTINE REDUCES A REAL SYMMETRIC (PACKED) MATRIX, STORED
1150 !* AS A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
1151 !* USING ORTHOGONAL SIMILARITY TRANSFORMATIONS, PRESERVING THE
1152 !* INFORMATION ABOUT THE TRANSFORMATIONS IN A.
1155 !* THE TRIDIAGONAL REDUCTION IS PERFORMED IN THE FOLLOWING WAY.
1156 !* STARTING WITH J=N, THE ELEMENTS IN THE J-TH ROW TO THE
1157 !* LEFT OF THE DIAGONAL ARE FIRST SCALED, TO AVOID POSSIBLE
1158 !* UNDERFLOW IN THE TRANSFORMATION THAT MIGHT RESULT IN SEVERE
1159 !* DEPARTURE FROM ORTHOGONALITY. THE SUM OF SQUARES SIGMA OF
1160 !* THESE SCALED ELEMENTS IS NEXT FORMED. THEN, A VECTOR U AND
1162 !* H = U(TRANSPOSE) * U / 2
1163 !* DEFINE A REFLECTION OPERATOR
1164 !* P = I - U * U(TRANSPOSE) / H
1165 !* WHICH IS ORTHOGONAL AND SYMMETRIC AND FOR WHICH THE
1166 !* SIMILIARITY TRANSFORMATION PAP ELIMINATES THE ELEMENTS IN
1167 !* THE J-TH ROW OF A TO THE LEFT OF THE SUBDIAGONAL AND THE
1168 !* SYMMETRICAL ELEMENTS IN THE J-TH COLUMN.
1170 !* THE NON-ZERO COMPONENTS OF U ARE THE ELEMENTS OF THE J-TH
1171 !* ROW TO THE LEFT OF THE DIAGONAL WITH THE LAST OF THEM
1172 !* AUGMENTED BY THE SQUARE ROOT OF SIGMA PREFIXED BY THE SIGN
1173 !* OF THE SUBDIAGONAL ELEMENT. BY STORING THE TRANSFORMED SUB-
1174 !* DIAGONAL ELEMENT IN E(J) AND NOT OVERWRITING THE ROW
1175 !* ELEMENTS ELIMINATED IN THE TRANSFORMATION, FULL INFORMATION
1176 !* ABOUT P IS SAVE FOR LATER USE IN ETRBK3.
1178 !* THE TRANSFORMATION SETS E2(J) EQUAL TO SIGMA AND E(J)
1179 !* EQUAL TO THE SQUARE ROOT OF SIGMA PREFIXED BY THE SIGN
1180 !* OF THE REPLACED SUBDIAGONAL ELEMENT.
1182 !* THE ABOVE STEPS ARE REPEATED ON FURTHER ROWS OF THE
1183 !* TRANSFORMED A IN REVERSE ORDER UNTIL A IS REDUCED TO TRI-
1184 !* DIAGONAL FORM, THAT IS, REPEATED FOR J = N-1,N-2,...,3.
1191 !* THE ORDER OF THE MATRIX.
1193 !* MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
1194 !* AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT
1195 !* A - W.P. REAL (NV)
1196 !* CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
1197 !* INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
1198 !* ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
1201 !* A - W.P. REAL (NV)
1202 !* CONTAINS INFORMATION ABOUT THE ORTHOGONAL
1203 !* TRANSFORMATIONS USED IN THE REDUCTION.
1204 !* D - W.P. REAL (N)
1205 !* CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL
1207 !* E - W.P. REAL (N)
1208 !* CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
1209 !* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO
1210 !* E2 - W.P. REAL (N)
1211 !* CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF
1212 !* E. MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
1214 !* DIFFERENCES FROM EISPACK 3 -
1215 !* OUTER LOOP CHANGED FROM II=1,N TO I=N,3,-1
1216 !* PARAMETER STATEMENT AND GENERIC INTRINSIC FUNCTIONS USED
1217 !* SCALE.NE.0 TEST NOW SPOTS TRI-DIAGONAL FORM
1218 !* VALUES LESS THAN EPSLON CLEARED TO ZERO
1220 !* U NOT COPIED TO D, LEFT IN A
1221 !* E2 COMPUTED FROM E
1222 !* INNER LOOPS SPLIT INTO ROUTINES ELAU AND FREDA
1223 !* INVERSE OF H STORED INSTEAD OF H
1226 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
1227 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
1229 INTEGER :: I,IIA,IZ0,L,N,NV
1231 real(kind=8) :: A(NV),D(N),E2(N)
1232 real(kind=8) :: E(*)!el E(L)
1233 real(kind=8) :: AIIMAX,F,G,H,HROOT,SCALE,SCALEI
1234 ! real(kind=8) :: DASUM, DNRM2
1236 real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00
1238 !-----------------------------------------------------------------------
1240 IF (N .LE. 2) GO TO 310
1242 AIIMAX = ABS(A(IZ0))
1247 AIIMAX = MAX(AIIMAX, ABS(A(IIA)))
1248 SCALE = DASUM (L, A(IZ0+1), 1)
1249 IF(SCALE .EQ. ABS(A(IIA-1)) .OR. AIIMAX+SCALE .EQ. AIIMAX) THEN
1251 ! THIS ROW IS ALREADY IN TRI-DIAGONAL FORM
1254 IF (AIIMAX+D(I) .EQ. AIIMAX) D(I) = ZERO
1256 IF (AIIMAX+E(I) .EQ. AIIMAX) E(I) = ZERO
1263 SCALEI = ONE / SCALE
1264 CALL DSCAL(L,SCALEI,A(IZ0+1),1)
1265 HROOT = DNRM2(L,A(IZ0+1),1)
1271 H = HROOT*HROOT - F * G
1274 A(IIA) = ONE / SQRT(H)
1275 ! .......... FORM P THEN Q IN E(1:L) ..........
1276 CALL ELAU(ONE/H,L,A(IZ0+1),A,E)
1277 ! .......... FORM REDUCED A ..........
1278 CALL FREDA(L,A(IZ0+1),A,E)
1291 end subroutine ETRED3
1292 !-----------------------------------------------------------------------------
1293 !*MODULE EIGEN *DECK EVVRSP
1294 subroutine EVVRSP(MSGFL,N,NVECT,LENA,NV,A,B,IND,ROOT,VECT,IORDER,IERR)
1296 !* AUTHOR: S. T. ELBERT, AMES LABORATORY-USDOE, JUNE 1985
1299 !* FINDS (ALL) EIGENVALUES AND (SOME OR ALL) EIGENVECTORS
1301 !* OF A REAL SYMMETRIC PACKED MATRIX.
1305 !* THE METHOD AS PRESENTED IN THIS ROUTINE CONSISTS OF FOUR STEPS:
1306 !* FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE
1307 !* HOUSEHOLDER TECHNIQUE (ORTHOGONAL SIMILARITY TRANSFORMATIONS).
1308 !* SECOND, THE ROOTS ARE LOCATED USING THE RATIONAL QL METHOD.
1309 !* THIRD, THE VECTORS OF THE TRIDIAGONAL FORM ARE EVALUATED BY THE
1310 !* INVERSE ITERATION TECHNIQUE. VECTORS FOR DEGENERATE OR NEAR-
1311 !* DEGENERATE ROOTS ARE FORCED TO BE ORTHOGONAL.
1312 !* FOURTH, THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE
1315 !* THESE ROUTINES ARE MODIFICATIONS OF THE EISPACK 3
1316 !* ROUTINES TRED3, TQLRAT, TINVIT AND TRBAK3
1318 !* FOR FURTHER DETAILS, SEE EISPACK USERS GUIDE, B. T. SMITH
1319 !* ET AL, SPRINGER-VERLAG, LECTURE NOTES IN COMPUTER SCIENCE,
1320 !* VOL. 6, 2-ND EDITION, 1976. ANOTHER GOOD REFERENCE IS
1321 !* THE SYMMETRIC EIGENVALUE PROBLEM BY B. N. PARLETT
1322 !* PUBLISHED BY PRENTICE-HALL, INC., ENGLEWOOD CLIFFS, N.J. (1980)
1325 !* MSGFL - INTEGER (LOGICAL UNIT NO.)
1326 !* FILE WHERE ERROR MESSAGES WILL BE PRINTED.
1327 !* IF MSGFL IS 0, ERROR MESSAGES WILL BE PRINTED ON LU 6.
1328 !* IF MSGFL IS NEGATIVE, NO ERROR MESSAGES PRINTED.
1330 !* ORDER OF MATRIX A.
1332 !* NUMBER OF VECTORS DESIRED. 0 .LE. NVECT .LE. N.
1334 !* DIMENSION OF A IN CALLING ROUTINE. MUST NOT BE LESS
1337 !* ROW DIMENSION OF VECT IN CALLING ROUTINE. N .LE. NV.
1338 !* A - WORKING PRECISION REAL (LENA)
1339 !* INPUT MATRIX, ROWS OF THE LOWER TRIANGLE PACKED INTO
1340 !* LINEAR ARRAY OF DIMENSION N*(N+1)/2. THE PACKED ORDER
1341 !* IS A(1,1), A(2,1), A(2,2), A(3,1), A(3,2), ...
1342 !* B - WORKING PRECISION REAL (N,8)
1343 !* SCRATCH ARRAY, 8*N ELEMENTS
1344 !* IND - INTEGER (N)
1345 !* SCRATCH ARRAY OF LENGTH N.
1347 !* ROOT ORDERING FLAG.
1348 !* = 0, ROOTS WILL BE PUT IN ASCENDING ORDER.
1349 !* = 2, ROOTS WILL BE PUT IN DESCENDING ORDER.
1352 !* A - DESTORYED. NOW HOLDS REFLECTION OPERATORS.
1353 !* ROOT - WORKING PRECISION REAL (N)
1354 !* ALL EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
1355 !* IF IORDER = 0, ROOT(1) .LE. ... .LE. ROOT(N)
1356 !* IF IORDER = 2, ROOT(1) .GE. ... .GE. ROOT(N)
1357 !* VECT - WORKING PRECISION REAL (NV,NVECT)
1358 !* EIGENVECTORS FOR ROOT(1), ..., ROOT(NVECT).
1360 !* = 0 IF NO ERROR DETECTED,
1361 !* = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
1362 !* = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
1363 !* (FAILURES SHOULD BE VERY RARE. CONTACT C. MOLER.)
1367 !el LOGICAL :: GOPARR,DSKWRK,MASWRK
1369 integer :: MSGFL,N,NVECT,LENA,NV,IORDER,IERR
1370 real(kind=8) :: A(LENA)
1371 real(kind=8) :: B(N,8)
1372 real(kind=8) :: ROOT(N)
1374 real(kind=8) :: VECT(NV,*)
1378 !el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
1379 real(kind=8) :: DSKW
1381 !el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1383 900 FORMAT(26H0*** EVVRSP PARAMETERS ***/ &
1384 14H *** N = ,I8,4H ***/ &
1385 14H *** NVECT = ,I8,4H ***/ &
1386 14H *** LENA = ,I8,4H ***/ &
1387 14H *** NV = ,I8,4H ***/ &
1388 14H *** IORDER = ,I8,4H ***/ &
1389 14H *** IERR = ,I8,4H ***)
1390 901 FORMAT(37H VALUE OF LENA IS LESS THAN (N*N+N)/2)
1391 902 FORMAT(39H EQLRAT HAS FAILED TO CONVERGE FOR ROOT,I5)
1392 903 FORMAT(18H NV IS LESS THAN N)
1393 904 FORMAT(41H EINVIT HAS FAILED TO CONVERGE FOR VECTOR,I5)
1394 905 FORMAT(51H VALUE OF IORDER MUST BE 0 (SMALLEST ROOT FIRST) OR &
1395 ,23H 2 (LARGEST ROOT FIRST))
1396 906 FORMAT(' VALUE OF N IS LESS THAN OR EQUAL ZERO')
1398 integer :: LMSGFL,I,J,L,JSV,KLIM,K
1400 !-----------------------------------------------------------------------
1403 IF (MSGFL .EQ. 0) LMSGFL=6
1405 IF (N .LE. 0) GO TO 800
1407 IF ( (N*N+N)/2 .GT. LENA) GO TO 810
1409 ! REDUCE REAL SYMMETRIC MATRIX A TO TRIDIAGONAL FORM
1411 CALL ETRED3(N,LENA,A,B(1,1),B(1,2),B(1,3))
1413 ! FIND ALL EIGENVALUES OF TRIDIAGONAL MATRIX
1415 CALL EQLRAT(N,B(1,1),B(1,2),B(1,3),ROOT,IND,IERR,B(1,4))
1416 IF (IERR .NE. 0) GO TO 820
1418 ! CHECK THE DESIRED ORDER OF THE EIGENVALUES
1421 IF (IORDER .EQ. 0) GO TO 300
1422 IF (IORDER .NE. 2) GO TO 850
1424 ! ORDER ROOTS IN DESCENDING ORDER (LARGEST FIRST)...
1425 ! TURN ROOT AND IND ARRAYS END FOR END
1437 ! FIND I AND J MARKING THE START AND END OF A SEQUENCE
1438 ! OF DEGENERATE ROOTS
1443 IF (I .GT. N) GO TO 300
1445 IF (ROOT(J) .NE. ROOT(I)) GO TO 240
1450 IF (J .EQ. I) GO TO 220
1452 ! TURN AROUND IND BETWEEN I AND J
1468 IF (NVECT .LE. 0) RETURN
1469 IF (NV .LT. N) GO TO 830
1471 ! FIND EIGENVECTORS OF TRI-DIAGONAL MATRIX VIA INVERSE ITERATION
1474 CALL EINVIT(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,IND,&
1475 VECT,IERR,B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
1476 IF (IERR .NE. 0) GO TO 840
1478 ! FIND EIGENVECTORS OF SYMMETRIC MATRIX VIA BACK TRANSFORMATION
1481 CALL ETRBK3(NV,N,LENA,A,NVECT,VECT)
1484 ! ERROR MESSAGE SECTION
1486 800 IF (LMSGFL .LT. 0) RETURN
1487 IF (MASWRK) WRITE(LMSGFL,906)
1490 810 IF (LMSGFL .LT. 0) RETURN
1491 IF (MASWRK) WRITE(LMSGFL,901)
1494 820 IF (LMSGFL .LT. 0) RETURN
1495 IF (MASWRK) WRITE(LMSGFL,902) IERR
1498 830 IF (LMSGFL .LT. 0) RETURN
1499 IF (MASWRK) WRITE(LMSGFL,903)
1503 IF ((LMSGFL .GT. 0).AND.MASWRK) WRITE(LMSGFL,904) -IERR
1507 IF (LMSGFL .LT. 0) RETURN
1508 IF (MASWRK) WRITE(LMSGFL,905)
1512 IF (MASWRK) WRITE(LMSGFL,900) N,NVECT,LENA,NV,IORDER,IERR
1514 end subroutine EVVRSP
1515 !-----------------------------------------------------------------------------
1516 !*MODULE EIGEN *DECK FREDA
1517 subroutine FREDA(L,D,A,E)
1520 real(kind=8) :: A(*)
1521 real(kind=8) :: D(L)
1522 real(kind=8) :: E(*)!el E(L)
1528 ! .......... FORM REDUCED A ..........
1535 A(JK) = A(JK) - F * E(K) - G * D(K)
1541 end subroutine FREDA
1542 !-----------------------------------------------------------------------------
1543 !*MODULE EIGEN *DECK GIVEIS
1544 subroutine GIVEIS(N,NVECT,NV,A,B,INDB,ROOT,VECT,IERR)
1546 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1547 integer :: N,NVECT,NV,IERR
1548 real(kind=8),DIMENSION(*) :: A
1549 real(kind=8),DIMENSION(N,8) :: B
1550 integer,DIMENSION(N) :: INDB
1551 real(kind=8),DIMENSION(N) :: ROOT
1552 real(kind=8),DIMENSION(NV,NVECT) :: VECT
1554 ! EISPACK-BASED SUBSTITUTE FOR QCPE ROUTINE GIVENS.
1555 ! FINDS ALL EIGENVALUES AND SOME EIGENVECTORS OF A REAL SYMMETRIC
1556 ! MATRIX. AUTHOR.. C. MOLER AND D. SPANGLER, N.R.C.C., 4/1/79.
1559 ! N = ORDER OF MATRIX .
1560 ! NVECT = NUMBER OF VECTORS DESIRED. 0 .LE. NVECT .LE. N .
1561 ! NV = LEADING DIMENSION OF VECT .
1562 ! A = INPUT MATRIX, COLUMNS OF THE UPPER TRIANGLE PACKED INTO
1563 ! LINEAR ARRAY OF DIMENSION N*(N+1)/2 .
1564 ! B = SCRATCH ARRAY, 8*N ELEMENTS (NOTE THIS IS MORE THAN
1565 ! PREVIOUS VERSIONS OF GIVENS.)
1566 ! IND = INDEX ARRAY OF N ELEMENTS
1570 ! ROOT = ALL EIGENVALUES, ROOT(1) .LE. ... .LE. ROOT(N) .
1571 ! (FOR OTHER ORDERINGS, SEE BELOW.)
1572 ! VECT = EIGENVECTORS FOR ROOT(1),..., ROOT(NVECT) .
1573 ! IERR = 0 IF NO ERROR DETECTED,
1574 ! = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
1575 ! = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
1576 ! (FAILURES SHOULD BE VERY RARE. CONTACT MOLER.)
1578 ! CALLS MODIFIED EISPACK ROUTINES TRED3B, IMTQLV, TINVTB, AND
1579 ! TRBK3B. THE ROUTINES TRED3B, TINVTB, AND TRBK3B.
1580 ! THE ORIGINAL EISPACK ROUTINES TRED3, TINVIT, AND TRBAK3
1581 ! WERE MODIFIED BY THE INTRODUCTION OF TWO ROUTINES FROM THE
1582 ! BLAS LIBRARY - DDOT AND DAXPY.
1584 ! IF TINVIT FAILS TO CONVERGE, TQL2 IS CALLED
1586 ! SEE EISPACK USERS GUIDE, B. T. SMITH ET AL, SPRINGER-VERLAG
1587 ! LECTURE NOTES IN COMPUTER SCIENCE, VOL. 6, 2-ND EDITION, 1976 .
1588 ! NOTE THAT IMTQLV AND TINVTB HAVE INTERNAL MACHINE
1589 ! DEPENDENT CONSTANTS.
1591 !el DATA ONE, ZERO /1.0D+00, 0.0D+00/
1592 real(kind=8) :: ZERO = 0.0D+00, ONE = 1.0D+00
1596 CALL TRED3B(N,(N*N+N)/2,A,B(1,1),B(1,2),B(1,3))
1597 CALL IMTQLV(N,B(1,1),B(1,2),B(1,3),ROOT,INDB,IERR,B(1,4))
1598 IF (IERR .NE. 0) RETURN
1600 ! TO REORDER ROOTS...
1610 IF (NVECT .LE. 0) RETURN
1611 CALL TINVTB(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,INDB,VECT,IERR,&
1612 B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
1613 IF (IERR .EQ. 0) GO TO 160
1615 ! IF INVERSE ITERATION GIVES AN ERROR IN DETERMINING THE
1616 ! EIGENVECTORS, TRY THE QL ALGORITHM IF ALL THE EIGENVECTORS
1619 IF (NVECT .NE. N) RETURN
1626 CALL TQL2 (NV,N,B(1,1),B(1,2),VECT,IERR)
1630 IF (IERR .NE. 0) RETURN
1631 160 CALL TRBK3B(NV,N,(N*N+N)/2,A,NVECT,VECT)
1633 end subroutine GIVEIS
1634 !-----------------------------------------------------------------------------
1635 !*MODULE EIGEN *DECK GLDIAG
1636 subroutine GLDIAG(LDVECT,NVECT,N,H,WRK,EIG,VECTOR,IERR,IWRK)
1638 ! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1643 !el LOGICAL :: GOPARR,DSKWRK,MASWRK
1645 integer :: LDVECT,NVECT,N,IERR
1646 real(kind=8),DIMENSION(*) :: H
1647 real(kind=8),DIMENSION(N,8) :: WRK
1648 real(kind=8),DIMENSION(N) :: EIG
1649 integer,DIMENSION(N) :: IWRK
1650 real(kind=8),DIMENSION(LDVECT,NVECT) :: VECTOR
1652 !el integer :: IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
1653 !el integer :: KDIAG,ICORFL,IXDR
1654 !el COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA
1655 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
1656 !el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
1657 !el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1659 integer :: LENH,KORDER
1661 ! ----- GENERAL ROUTINE TO DIAGONALIZE A SYMMETRIC MATRIX -----
1662 ! IF KDIAG = 0, USE A ROUTINE FROM THE VECTOR LIBRARY,
1663 ! IF AVAILABLE (SEE THE SUBROUTINE 'GLDIAG'
1664 ! IN VECTOR.SRC), OR EVVRSP OTHERWISE
1669 ! N = DIMENSION (ORDER) OF MATRIX TO BE SOLVED
1670 ! LDVECT = LEADING DIMENSION OF VECTOR
1671 ! NVECT = NUMBER OF VECTORS DESIRED
1672 ! H = MATRIX TO BE DIAGONALIZED
1673 ! WRK = N*8 W.P. REAL WORDS OF SCRATCH SPACE
1674 ! EIG = EIGENVECTORS (OUTPUT)
1675 ! VECTOR = EIGENVECTORS (OUTPUT)
1676 ! IERR = ERROR FLAG (OUTPUT)
1677 ! IWRK = N INTEGER WORDS OF SCRATCH SPACE
1681 ! ----- USE STEVE ELBERT'S ROUTINE -----
1683 IF(KDIAG.LE.1 .OR. KDIAG.GT.3) THEN
1686 CALL EVVRSP(IW,N,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR, &
1690 ! ----- USE MODIFIED EISPAK ROUTINE -----
1693 CALL GIVEIS(N,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
1695 ! ----- USE JACOBI ROTATION ROUTINE -----
1699 CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,N)
1701 IF (MASWRK) WRITE(IW,9000) N,NVECT,LDVECT
1707 9000 FORMAT(1X,'IN -GLDIAG-, N,NVECT,LDVECT=',3I8/ &
1708 1X,'THE JACOBI CODE CANNOT COPE WITH N.NE.NVECT!'/ &
1709 1X,'SO THIS RUN DOES NOT PERMIT KDIAG=3.')
1710 end subroutine GLDIAG
1711 !-----------------------------------------------------------------------------
1712 !*MODULE EIGEN *DECK IMTQLV
1713 subroutine IMTQLV(N,D,E,E2,W,IND,IERR,RV1)
1715 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1716 INTEGER :: N,TAG,IERR
1717 real(kind=8) :: MACHEP
1718 real(kind=8),DIMENSION(N) :: D,E2,W,RV1
1719 real(kind=8) :: E(*)!el E(L)
1720 integer,DIMENSION(N) :: IND
1721 integer :: k,i,l,j,m,mml,ii
1722 real(kind=8) :: c,p,s,f,b,g,r
1724 ! THIS ROUTINE IS A VARIANT OF IMTQL1 WHICH IS A TRANSLATION OF
1725 ! ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND
1726 ! WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
1727 ! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
1729 ! THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
1730 ! MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM
1731 ! THEIR CORRESPONDING SUBMATRIX INDICES.
1735 ! N IS THE ORDER OF THE MATRIX,
1737 ! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1739 ! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1740 ! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
1742 ! E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
1743 ! E2(1) IS ARBITRARY.
1747 ! D AND E ARE UNALTERED,
1749 ! ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
1750 ! AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
1751 ! MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
1752 ! E2(1) IS ALSO SET TO ZERO,
1754 ! W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
1755 ! ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
1756 ! ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
1757 ! THE SMALLEST EIGENVALUES,
1759 ! IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE
1760 ! CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES
1761 ! BELONGING TO THE FIRST SUBMATRIX FROM THE TOP,
1762 ! 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.,
1765 ! ZERO FOR NORMAL RETURN,
1766 ! J IF THE J-TH EIGENVALUE HAS NOT BEEN
1767 ! DETERMINED AFTER 30 ITERATIONS,
1769 ! RV1 IS A TEMPORARY STORAGE ARRAY.
1771 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1772 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1774 ! ------------------------------------------------------------------
1776 ! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1777 ! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1780 MACHEP = 2.0D+00**(-50)
1788 IF (I .NE. 1) RV1(I-1) = E(I)
1796 ! ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
1798 IF (M .EQ. N) GO TO 160
1799 IF (ABS(RV1(M)) .LE. MACHEP * (ABS(W(M)) + ABS(W(M+1)))) GO TO &
1801 ! ********** GUARD AGAINST UNDERFLOWED ELEMENT OF E2 **********
1802 IF (E2(M+1) .EQ. 0.0D+00) GO TO 180
1805 160 IF (M .LE. K) GO TO 200
1806 IF (M .NE. N) E2(M+1) = 0.0D+00
1810 IF (M .EQ. L) GO TO 280
1811 IF (J .EQ. 30) GO TO 380
1813 ! ********** FORM SHIFT **********
1814 G = (W(L+1) - P) / (2.0D+00 * RV1(L))
1815 R = SQRT(G*G+1.0D+00)
1816 G = W(M) - P + RV1(L) / (G + SIGN(R,G))
1821 ! ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
1826 IF (ABS(F) .LT. ABS(G)) GO TO 220
1828 R = SQRT(C*C+1.0D+00)
1834 R = SQRT(S*S+1.0D+00)
1839 R = (W(I) - G) * S + 2.0D+00 * C * B
1849 ! ********** ORDER EIGENVALUES **********
1850 280 IF (L .EQ. 1) GO TO 320
1851 ! ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
1854 IF (P .GE. W(I-1)) GO TO 340
1865 ! ********** SET ERROR -- NO CONVERGENCE TO AN
1866 ! EIGENVALUE AFTER 30 ITERATIONS **********
1869 ! ********** LAST CARD OF IMTQLV **********
1870 end subroutine IMTQLV
1871 !-----------------------------------------------------------------------------
1872 !*MODULE EIGEN *DECK JACDG
1873 subroutine JACDG(A,VEC,EIG,JBIG,BIG,LDVEC,N)
1875 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1878 real(kind=8),DIMENSION(*) :: A
1879 real(kind=8),DIMENSION(LDVEC,N) :: VEC
1880 real(kind=8),DIMENSION(N) :: EIG,BIG
1881 integer,DIMENSION(N) :: JBIG
1883 real(kind=8),PARAMETER :: ONE=1.0D+00
1884 integer :: i,NB1,NB2,NMIN,NMAX
1886 ! ----- JACOBI DIAGONALIZATION OF SYMMETRIC MATRIX -----
1887 ! SYMMETRIC MATRIX -A- OF DIMENSION -N- IS DESTROYED ON EXIT.
1888 ! ALL EIGENVECTORS ARE FOUND, SO -VEC- MUST BE SQUARE,
1889 ! UNLESS SOMEONE TAKES THE TROUBLE TO LOOK AT -NMAX- BELOW.
1890 ! -BIG- AND -JBIG- ARE SCRATCH WORK ARRAYS.
1892 CALL VCLR(VEC,1,LDVEC*N)
1898 NB2 = (NB1*NB1+NB1)/2
1902 CALL JACDIA(A,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1905 EIG(I) = A((I*I+I)/2)
1908 CALL JACORD(VEC,EIG,NB1,LDVEC)
1910 end subroutine JACDG
1911 !-----------------------------------------------------------------------------
1912 !*MODULE EIGEN *DECK JACDIA
1913 subroutine JACDIA(F,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1915 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
1917 integer :: NB1,NB2,LDVEC,NMIN,NMAX
1918 !el LOGICAL :: GOPARR,DSKWRK,MASWRK
1919 real(kind=8),DIMENSION(NB2) :: F
1920 real(kind=8),DIMENSION(LDVEC,NB1) :: VEC
1921 real(kind=8),DIMENSION(NB1) :: BIG
1922 integer,DIMENSION(NB1) :: JBIG
1924 !el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
1925 !el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
1927 real(kind=8),PARAMETER :: ROOT2=0.707106781186548D+00
1928 real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00, D1050=1.05D+00,&
1929 D1500=1.5D+00, D3875=3.875D+00,&
1930 D0500=0.5D+00, D1375=1.375D+00, D0250=0.25D+00
1931 real(kind=8),PARAMETER :: C2=1.0D-12, C3=4.0D-16,&
1932 C4=2.0D-16, C5=8.0D-09, C6=3.0D-06
1933 integer :: i,ii,j,k,jj,IEAA,IEAB,MAXIT,ITER,I1,IA,IB,IAA,IBB,IEAR,&
1935 real(kind=8) :: T,TT,EPS,SD,TEST,DIF,CX,SX,T2X2,T2X25,T1,T2
1937 ! F IS THE MATRIX TO BE DIAGONALIZED, F IS STORED TRIANGULAR
1938 ! VEC IS THE ARRAY OF EIGENVECTORS, DIMENSION NB1*NB1
1939 ! BIG AND JBIG ARE TEMPORARY SCRATCH AREAS OF DIMENSION NB1
1940 ! THE ROTATIONS AMONG THE FIRST NMIN BASIS FUNCTIONS ARE NOT
1942 ! THE ROTATIONS AMONG THE LAST NB1-NMAX BASIS FUNCTIONS ARE NOT
1948 EPS = 64.0D+00*EPSLON(ONE)
1950 ! LOOP OVER COLUMNS (K) OF TRIANGULAR MATRIX TO DETERMINE
1951 ! LARGEST OFF-DIAGONAL ELEMENTS IN ROW(I).
1956 IF(I.LT.NMIN .OR. I.EQ.1) GO TO 20
1960 IF(ABS(BIG(I)).GE.ABS(F(II+K))) GO TO 10
1966 ! ----- 2X2 JACOBI ITERATIONS BEGIN HERE -----
1968 MAXIT=MAX(NB2*20,500)
1973 ! FIND SMALLEST DIAGONAL ELEMENT
1979 SD= MIN(SD,ABS(F(JJ)))
1981 TEST = MAX(EPS, C2*MAX(SD,C6))
1983 ! FIND LARGEST OFF-DIAGONAL ELEMENT
1989 IF(T.GE.ABS(BIG(I))) GO TO 50
1994 ! TEST FOR CONVERGENCE, THEN DETERMINE ROTATION.
1996 IF(T.LT.TEST) RETURN
1999 IF(ITER.GT.MAXIT) THEN
2001 WRITE(6,*) 'JACOBI DIAGONALIZATION FAILS, DIMENSION=',NB1
2002 WRITE(6,9020) ITER,T,TEST,SD
2011 DIF=F(IAA+IA)-F(IBB+IB)
2012 IF(ABS(DIF).GT.C3*T) GO TO 70
2018 IF(T2X25 .GT. C4) GO TO 80
2022 80 IF(T2X25 .GT. C5) GO TO 90
2023 SX=T2X2*(ONE-D1500*T2X25)
2026 90 IF(T2X25 .GT. C6) GO TO 100
2027 CX=ONE+T2X25*(T2X25*D1375 - D0500)
2028 SX= T2X2*(ONE + T2X25*(T2X25*D3875 - D1500))
2030 100 T=D0250 / SQRT(D0250 + T2X25)
2032 SX= SIGN( SQRT(D0500 - T),T2X2)
2038 F(IEAR)=F(IEAR)*CX+F(IEBR)*SX
2039 F(IEBR)=T-F(IEBR)*CX
2040 IF(IR-IA) 220,120,130
2046 IF(JBIG(IR)) 200,220,200
2050 IF(IR-IB) 180,150,160
2051 150 F(IEAA)=F(IEAA)*CX+F(IEAB)*SX
2052 F(IEAB)=TT*CX+F(IEBR)*SX
2053 F(IEBR)=TT*SX-F(IEBR)*CX
2056 160 IF( ABS(T) .GE. ABS(F(IEBR))) GO TO 170
2057 IF(IB.GT.NMAX) GO TO 170
2061 180 IF( ABS(T) .LT. ABS(BIG(IR))) GO TO 190
2065 190 IF(IA .NE. JBIG(IR) .AND. IB .NE. JBIG(IR)) GO TO 220
2071 IF(ABS(BIG(IR)) .GE. ABS(F(K))) GO TO 210
2079 T1=VEC(I,IA)*CX + VEC(I,IB)*SX
2080 T2=VEC(I,IA)*SX - VEC(I,IB)*CX
2086 9020 FORMAT(1X,'ITER=',I6,' T,TEST,SD=',1P,3E20.10)
2087 end subroutine JACDIA
2088 !-----------------------------------------------------------------------------
2089 !*MODULE EIGEN *DECK JACORD
2090 subroutine JACORD(VEC,EIG,N,LDVEC)
2092 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2094 real(kind=8),DIMENSION(LDVEC,N) :: VEC
2095 real(kind=8),DIMENSION(N) :: EIG
2099 ! ---- SORT EIGENDATA INTO ASCENDING ORDER -----
2104 IF (EIG(J) .LT. EIG(JJ)) JJ = J
2106 IF (JJ .EQ. I) GO TO 290
2112 VEC(J,JJ) = VEC(J,I)
2117 end subroutine JACORD
2118 !-----------------------------------------------------------------------------
2119 !*MODULE EIGEN *DECK TINVTB
2120 subroutine TINVTB(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
2122 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2123 integer :: NM,N,M,IERR
2124 real(kind=8),DIMENSION(N) :: D,E,E2
2125 real(kind=8),DIMENSION(M) :: W
2126 real(kind=8),DIMENSION(NM,M) :: Z
2127 real(kind=8),DIMENSION(N) :: RV1,RV2,RV3,RV4,RV6
2128 integer,DIMENSION(M) :: IND
2129 real(kind=8) :: MACHEP,NORM
2130 INTEGER :: P,Q,R,S,TAG,GROUP
2131 integer :: ii,j,jj,i,iqmp,its
2132 real(kind=8) :: ORDER,XU,UK,X0,U,EPS2,EPS3,EPS4,x1,v
2134 ! ------------------------------------------------------------------
2136 ! THIS ROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH-
2137 ! NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
2138 ! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
2140 ! THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
2141 ! SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
2142 ! USING INVERSE ITERATION.
2146 ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2147 ! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2148 ! DIMENSION STATEMENT,
2150 ! N IS THE ORDER OF THE MATRIX,
2152 ! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2154 ! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2155 ! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
2157 ! E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E,
2158 ! WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
2159 ! E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
2160 ! THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM
2161 ! OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST CONTAIN
2162 ! 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0
2163 ! IF THE EIGENVALUES ARE IN DESCENDING ORDER. IF BISECT,
2164 ! TRIDIB, OR IMTQLV HAS BEEN USED TO FIND THE EIGENVALUES,
2165 ! THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE,
2167 ! M IS THE NUMBER OF SPECIFIED EIGENVALUES,
2169 ! W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER,
2171 ! IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
2172 ! ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
2173 ! 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
2174 ! THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.
2178 ! ALL INPUT ARRAYS ARE UNALTERED,
2180 ! Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
2181 ! ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO,
2184 ! ZERO FOR NORMAL RETURN,
2185 ! -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
2186 ! EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS,
2188 ! RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
2190 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2191 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2193 ! ------------------------------------------------------------------
2195 ! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2196 ! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2199 MACHEP = 2.0D+00**(-50)
2202 IF (M .EQ. 0) GO TO 680
2204 ORDER = 1.0D+00 - E2(1)
2214 ! ********** ESTABLISH AND PROCESS NEXT SUBMATRIX **********
2219 IF (Q .EQ. N) GO TO 140
2220 IF (E2(Q+1) .EQ. 0.0D+00) GO TO 140
2222 ! ********** FIND VECTORS BY INVERSE ITERATION **********
2228 IF (IND(R) .NE. TAG) GO TO 660
2231 IF (S .NE. 0) GO TO 220
2232 ! ********** CHECK FOR ISOLATED ROOT **********
2234 IF (P .NE. Q) GO TO 160
2237 160 NORM = ABS(D(P))
2240 180 NORM = NORM + ABS(D(I)) + ABS(E(I))
2241 ! ********** EPS2 IS THE CRITERION FOR GROUPING,
2242 ! EPS3 REPLACES ZERO PIVOTS AND EQUAL
2243 ! ROOTS ARE MODIFIED BY EPS3,
2244 ! EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW **********
2245 EPS2 = 1.0D-03 * NORM
2246 EPS3 = MACHEP * NORM
2249 UK = EPS4 / SQRT(UK)
2253 ! ********** LOOK FOR CLOSE OR COINCIDENT ROOTS **********
2254 220 IF (ABS(X1-X0) .GE. EPS2) GO TO 200
2256 IF (ORDER * (X1 - X0) .LE. 0.0D+00) X1 = X0 + ORDER * EPS3
2257 ! ********** ELIMINATION WITH INTERCHANGES AND
2258 ! INITIALIZATION OF VECTOR **********
2263 IF (I .EQ. P) GO TO 280
2264 IF (ABS(E(I)) .LT. ABS(U)) GO TO 260
2265 ! ********** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
2266 ! E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY **********
2270 RV2(I-1) = D(I) - X1
2272 IF (I .NE. Q) RV3(I-1) = E(I+1)
2273 U = V - XU * RV2(I-1)
2281 280 U = D(I) - X1 - XU * V
2282 IF (I .NE. Q) V = E(I+1)
2285 IF (U .EQ. 0.0D+00) U = EPS3
2289 ! ********** BACK SUBSTITUTION
2290 ! FOR I=Q STEP -1 UNTIL P DO -- **********
2291 320 DO 340 II = P, Q
2293 RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
2297 ! ********** ORTHOGONALIZE WITH RESPECT TO PREVIOUS
2298 ! MEMBERS OF GROUP **********
2299 IF (GROUP .EQ. 0) GO TO 400
2302 DO 380 JJ = 1, GROUP
2304 IF (IND(J) .NE. TAG) GO TO 360
2305 XU = DDOT(IQMP,RV6(P),1,Z(P,J),1)
2307 CALL DAXPY(IQMP,-XU,Z(P,J),1,RV6(P),1)
2314 420 NORM = NORM + ABS(RV6(I))
2316 IF (NORM .GE. 1.0D+00) GO TO 560
2317 ! ********** FORWARD SUBSTITUTION **********
2318 IF (ITS .EQ. 5) GO TO 540
2319 IF (NORM .NE. 0.0D+00) GO TO 440
2324 440 XU = EPS4 / NORM
2327 460 RV6(I) = RV6(I) * XU
2328 ! ********** ELIMINATION OPERATIONS ON NEXT VECTOR
2329 ! ITERATE **********
2330 480 DO 520 I = IP(1), Q
2332 ! ********** IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
2333 ! WAS PERFORMED EARLIER IN THE
2334 ! TRIANGULARIZATION PROCESS **********
2335 IF (RV1(I-1) .NE. E(I)) GO TO 500
2338 500 RV6(I) = U - RV4(I) * RV6(I-1)
2343 ! ********** SET ERROR -- NON-CONVERGED EIGENVECTOR **********
2347 ! ********** NORMALIZE SO THAT SUM OF SQUARES IS
2348 ! 1 AND EXPAND TO FULL ORDER **********
2352 RV6(I) = RV6(I) / NORM
2353 580 U = U + RV6(I)**2
2355 XU = 1.0D+00 / SQRT(U)
2358 620 Z(I,R) = 0.0D+00
2361 640 Z(I,R) = RV6(I) * XU
2366 IF (Q .LT. N) GO TO 100
2368 ! ********** LAST CARD OF TINVIT **********
2369 end subroutine TINVTB
2370 !-----------------------------------------------------------------------------
2371 !*MODULE EIGEN *DECK TQL2
2372 subroutine TQL2(NM,N,D,E,Z,IERR)
2374 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2375 integer :: NM,N,IERR
2376 real(kind=8) :: MACHEP
2377 real(kind=8),DIMENSION(N) :: D!,E
2378 real(kind=8) :: E(*)!el E(L)
2379 real(kind=8),DIMENSION(NM,N) :: Z
2380 integer :: ii,i,j,mml,m,l1,k,l
2381 real(kind=8) :: c,f,b,h,g,p,r,s
2383 ! THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
2384 ! NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
2386 ! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
2388 ! THIS ROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
2389 ! OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
2390 ! THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
2391 ! BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS
2392 ! FULL MATRIX TO TRIDIAGONAL FORM.
2396 ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2397 ! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2398 ! DIMENSION STATEMENT,
2400 ! N IS THE ORDER OF THE MATRIX,
2402 ! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2404 ! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2405 ! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
2407 ! Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
2408 ! REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS
2409 ! OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
2410 ! THE IDENTITY MATRIX.
2414 ! D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
2415 ! ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
2416 ! UNORDERED FOR INDICES 1,2,...,IERR-1,
2418 ! E HAS BEEN DESTROYED,
2420 ! Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
2421 ! TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE,
2422 ! Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
2426 ! ZERO FOR NORMAL RETURN,
2427 ! J IF THE J-TH EIGENVALUE HAS NOT BEEN
2428 ! DETERMINED AFTER 30 ITERATIONS.
2430 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2431 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2433 ! ------------------------------------------------------------------
2435 ! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2436 ! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2439 MACHEP = 2.0D+00**(-50)
2442 IF (N .EQ. 1) GO TO 400
2453 H = MACHEP * (ABS(D(L)) + ABS(E(L)))
2455 ! ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
2457 IF (ABS(E(M)) .LE. B) GO TO 140
2458 ! ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
2459 ! THROUGH THE BOTTOM OF THE LOOP **********
2462 140 IF (M .EQ. L) GO TO 280
2463 160 IF (J .EQ. 30) GO TO 380
2465 ! ********** FORM SHIFT **********
2468 P = (D(L1) - G) / (2.0D+00 * E(L))
2469 R = SQRT(P*P+1.0D+00)
2470 D(L) = E(L) / (P + SIGN(R,P))
2477 ! ********** QL TRANSFORMATION **********
2482 ! ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
2487 IF (ABS(P) .LT. ABS(E(I))) GO TO 200
2489 R = SQRT(C*C+1.0D+00)
2495 R = SQRT(C*C+1.0D+00)
2496 E(I+1) = S * E(I) * R
2499 220 P = C * D(I) - S * G
2500 D(I+1) = H + S * (C * G + S * D(I))
2501 ! ********** FORM VECTOR **********
2502 CALL DROT(N,Z(1,I+1),1,Z(1,I),1,C,S)
2508 IF (ABS(E(L)) .GT. B) GO TO 160
2511 ! ********** ORDER EIGENVALUES AND EIGENVECTORS **********
2518 IF (D(J) .GE. P) GO TO 320
2523 IF (K .EQ. I) GO TO 360
2527 CALL DSWAP(N,Z(1,I),1,Z(1,K),1)
2532 ! ********** SET ERROR -- NO CONVERGENCE TO AN
2533 ! EIGENVALUE AFTER 30 ITERATIONS **********
2536 ! ********** LAST CARD OF TQL2 **********
2538 !-----------------------------------------------------------------------------
2539 !*MODULE EIGEN *DECK TRBK3B
2540 subroutine TRBK3B(NM,N,NV,A,M,Z)
2542 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2543 integer :: NM,N,NV,M
2544 real(kind=8),DIMENSION(NV) :: A
2545 real(kind=8),DIMENSION(NM,M) :: Z
2546 integer :: i,l,iz,ik,j
2549 ! THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
2550 ! NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
2551 ! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
2553 ! THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
2554 ! MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
2555 ! SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRED3B.
2559 ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2560 ! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2561 ! DIMENSION STATEMENT,
2563 ! N IS THE ORDER OF THE MATRIX,
2565 ! NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2566 ! AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2568 ! A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS
2569 ! USED IN THE REDUCTION BY TRED3B IN ITS FIRST
2570 ! N*(N+1)/2 POSITIONS,
2572 ! M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
2574 ! Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
2575 ! IN ITS FIRST M COLUMNS.
2579 ! Z CONTAINS THE TRANSFORMED EIGENVECTORS
2580 ! IN ITS FIRST M COLUMNS.
2582 ! NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
2584 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2585 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2587 ! ------------------------------------------------------------------
2589 IF (M .EQ. 0) GO TO 140
2590 IF (N .EQ. 1) GO TO 140
2597 IF (H .EQ. 0.0D+00) GO TO 120
2600 S = -DDOT(L,A(IZ+1),1,Z(1,J),1)
2602 ! ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
2605 CALL DAXPY(L,S,A(IZ+1),1,Z(1,J),1)
2612 ! ********** LAST CARD OF TRBAK3 **********
2613 end subroutine TRBK3B
2614 !-----------------------------------------------------------------------------
2615 !*MODULE EIGEN *DECK TRED3B
2616 subroutine TRED3B(N,NV,A,D,E,E2)
2618 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2620 real(kind=8),DIMENSION(NV) :: A
2621 real(kind=8),DIMENSION(N) :: D,E2
2622 real(kind=8) :: E(*)!el E(L)
2623 integer :: ii,i,l,iz,k,jk,j,jm1
2624 real(kind=8) :: h,f,g,scale,dt,hh
2626 ! THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
2627 ! NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
2628 ! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
2630 ! THIS ROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS
2631 ! A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
2632 ! USING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
2636 ! N IS THE ORDER OF THE MATRIX,
2638 ! NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2639 ! AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2641 ! A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
2642 ! INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
2643 ! ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
2647 ! A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
2648 ! TRANSFORMATIONS USED IN THE REDUCTION,
2650 ! D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
2652 ! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
2653 ! MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO,
2655 ! E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
2656 ! E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
2658 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2659 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2661 ! ------------------------------------------------------------------
2663 ! ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
2670 IF (L .LT. 1) GO TO 120
2671 ! ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
2675 SCALE = SCALE + ABS(D(K))
2678 IF (SCALE .NE. 0.0D+00) GO TO 140
2688 E2(I) = SCALE * SCALE * H
2690 G = -SIGN(SQRT(H),F)
2694 A(IZ) = SCALE * D(L)
2695 IF (L .EQ. 1) GO TO 280
2703 ! ********** FORM ELEMENT OF A*U **********
2704 IF (JM1 .EQ. 0) GO TO 200
2706 E(K) = E(K) + DT * A(JK)
2707 G = G + D(K) * A(JK)
2710 200 E(J) = G + A(JK) * DT
2712 ! ********** FORM ELEMENT OF P **********
2722 ! ********** FORM REDUCED A **********
2730 A(JK) = A(JK) - F * E(K) - G * D(K)
2734 A(IZ+1) = SCALE * SQRT(H)
2738 ! ********** LAST CARD OF TRED3 **********
2739 end subroutine TRED3B
2740 !-----------------------------------------------------------------------------
2742 !-----------------------------------------------------------------------------
2743 ! 10 NOV 94 - MWS - DNRM2: REMOVE FTNCHECK WARNINGS
2744 ! 11 JUN 94 - MWS - INCLUDE A COPY OF DGEMV (LEVEL TWO ROUTINE)
2745 ! 11 AUG 87 - MWS - SANITIZE FLOATING POINT CONSTANTS IN DNRM2
2746 ! 26 MAR 87 - MWS - USE GENERIC SIGN IN DROTG
2747 ! 28 NOV 86 - STE - SUPPLY ALL LEVEL ONE BLAS
2748 ! 7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
2750 ! BASIC LINEAR ALGEBRA SUBPROGRAMS (BLAS) FROM LINPACK (LEVEL 1)
2752 ! THIS MODULE SHOULD BE COMPILED ONLY IF SPECIALLY CODED
2753 ! VERSIONS OF THESE ROUTINES ARE NOT AVAILABLE ON THE TARGET MACHINE
2755 !*MODULE BLAS1 *DECK DASUM
2756 real(kind=8) function DASUM(N,DX,INCX)
2758 ! TAKES THE SUM OF THE ABSOLUTE VALUES.
2759 ! JACK DONGARRA, LINPACK, 3/11/78.
2761 real(kind=8) :: DX(1),DTEMP
2762 INTEGER :: I,INCX,M,MP1,N,NINCX
2767 IF(INCX.EQ.1)GO TO 20
2769 ! CODE FOR INCREMENT NOT EQUAL TO 1
2772 DO 10 I = 1,NINCX,INCX
2773 DTEMP = DTEMP + ABS(DX(I))
2778 ! CODE FOR INCREMENT EQUAL TO 1
2784 IF( M .EQ. 0 ) GO TO 40
2786 DTEMP = DTEMP + ABS(DX(I))
2788 IF( N .LT. 6 ) GO TO 60
2791 DTEMP = DTEMP + ABS(DX(I)) + ABS(DX(I + 1)) + ABS(DX(I + 2)) &
2792 + ABS(DX(I + 3)) + ABS(DX(I + 4)) + ABS(DX(I + 5))
2797 !-----------------------------------------------------------------------------
2798 !*MODULE BLAS1 *DECK DAXPY
2799 subroutine DAXPY(N,DA,DX,INCX,DY,INCY)
2801 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2802 integer :: N,INCX,INCY
2803 real(kind=8),DIMENSION(1) :: DX,DY
2806 ! CONSTANT TIMES A VECTOR PLUS A VECTOR.
2807 ! DY(I) = DY(I) + DA * DX(I)
2808 ! USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
2809 ! JACK DONGARRA, LINPACK, 3/11/78.
2811 integer :: ix,iy,i,m,mp1
2813 IF (DA .EQ. 0.0D+00) RETURN
2814 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2816 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2821 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2822 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2824 DY(IY) = DY(IY) + DA*DX(IX)
2830 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
2836 IF( M .EQ. 0 ) GO TO 40
2838 DY(I) = DY(I) + DA*DX(I)
2840 IF( N .LT. 4 ) RETURN
2843 DY(I) = DY(I) + DA*DX(I)
2844 DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
2845 DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
2846 DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
2849 end subroutine DAXPY
2850 !-----------------------------------------------------------------------------
2851 !*MODULE BLAS1 *DECK DCOPY
2852 subroutine DCOPY(N,DX,INCX,DY,INCY)
2854 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2855 integer :: N,INCX,INCY
2856 real(kind=8),DIMENSION(*) :: DX,DY
2860 ! USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
2861 ! JACK DONGARRA, LINPACK, 3/11/78.
2863 integer :: ix,iy,m,i,mp1
2865 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2867 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2872 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2873 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2881 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
2887 IF( M .EQ. 0 ) GO TO 40
2891 IF( N .LT. 7 ) RETURN
2895 DY(I + 1) = DX(I + 1)
2896 DY(I + 2) = DX(I + 2)
2897 DY(I + 3) = DX(I + 3)
2898 DY(I + 4) = DX(I + 4)
2899 DY(I + 5) = DX(I + 5)
2900 DY(I + 6) = DX(I + 6)
2903 end subroutine DCOPY
2904 !-----------------------------------------------------------------------------
2905 !*MODULE BLAS1 *DECK DDOT
2906 real(kind=8) function DDOT(N,DX,INCX,DY,INCY)
2908 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2909 integer :: N,INCX,INCY
2910 real(kind=8),DIMENSION(1) :: DX,DY
2912 ! FORMS THE DOT PRODUCT OF TWO VECTORS.
2913 ! DOT = DX(I) * DY(I)
2914 ! USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
2915 ! JACK DONGARRA, LINPACK, 3/11/78.
2917 integer ::ix,iy,m,mp1,i
2918 real(kind=8) :: DTEMP
2922 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2924 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2929 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2930 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2932 DTEMP = DTEMP + DX(IX)*DY(IY)
2939 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
2945 IF( M .EQ. 0 ) GO TO 40
2947 DTEMP = DTEMP + DX(I)*DY(I)
2949 IF( N .LT. 5 ) GO TO 60
2952 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + &
2953 DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
2958 !-----------------------------------------------------------------------------
2959 !*MODULE BLAS1 *DECK DNRM2
2960 real(kind=8) function DNRM2(N,DX,INCX)
2962 INTEGER :: NEXT,N,INCX
2963 real(kind=8) :: DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
2964 DATA ZERO, ONE /0.0D+00, 1.0D+00/
2968 ! EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
2970 ! IF N .LE. 0 RETURN WITH RESULT = 0.
2971 ! IF N .GE. 1 THEN INCX MUST BE .GE. 1
2973 ! C.L.LAWSON, 1978 JAN 08
2975 ! FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE
2976 ! HOPEFULLY APPLICABLE TO ALL MACHINES.
2977 ! CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES.
2978 ! CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES.
2980 ! EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
2981 ! U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT)
2982 ! V = LARGEST NO. (OVERFLOW LIMIT)
2984 ! BRIEF OUTLINE OF ALGORITHM..
2986 ! PHASE 1 SCANS ZERO COMPONENTS.
2987 ! MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
2988 ! MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
2989 ! MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
2990 ! WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
2992 ! VALUES FOR CUTLO AND CUTHI..
2993 ! FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
2994 ! DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
2995 ! CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE
2996 ! UNIVAC AND DEC AT 2**(-103)
2997 ! THUS CUTLO = 2**(-51) = 4.44089E-16
2998 ! CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
2999 ! THUS CUTHI = 2**(63.5) = 1.30438E19
3000 ! CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
3001 ! THUS CUTLO = 2**(-33.5) = 8.23181D-11
3002 ! CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D+19
3003 ! DATA CUTLO, CUTHI / 8.232D-11, 1.304D+19 /
3004 ! DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
3005 DATA CUTLO, CUTHI / 8.232D-11, 1.304D+19 /
3008 IF(N .GT. 0) GO TO 10
3012 10 ASSIGN 30 TO NEXT
3017 20 GO TO NEXT,(30, 50, 70, 110)
3018 30 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
3022 ! PHASE 1. SUM IS ZERO
3024 50 IF( DX(I) .EQ. ZERO) GO TO 200
3025 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
3027 ! PREPARE FOR PHASE 2.
3031 ! PREPARE FOR PHASE 4.
3035 SUM = (SUM / DX(I)) / DX(I)
3036 105 XMAX = ABS(DX(I))
3039 ! PHASE 2. SUM IS SMALL.
3040 ! SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
3042 70 IF( ABS(DX(I)) .GT. CUTLO ) GO TO 75
3044 ! COMMON CODE FOR PHASES 2 AND 4.
3045 ! IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
3047 110 IF( ABS(DX(I)) .LE. XMAX ) GO TO 115
3048 SUM = ONE + SUM * (XMAX / DX(I))**2
3052 115 SUM = SUM + (DX(I)/XMAX)**2
3056 ! PREPARE FOR PHASE 3.
3058 75 SUM = (SUM * XMAX) * XMAX
3061 ! FOR REAL OR D.P. SET HITEST = CUTHI/N
3062 ! FOR COMPLEX SET HITEST = CUTHI/(2*N)
3066 ! PHASE 3. SUM IS MID-RANGE. NO SCALING.
3069 IF(ABS(DX(J)) .GE. HITEST) GO TO 100
3070 95 SUM = SUM + DX(J)**2
3076 IF ( I .LE. NN ) GO TO 20
3080 ! COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
3082 DNRM2 = XMAX * SQRT(SUM)
3086 !-----------------------------------------------------------------------------
3087 !*MODULE BLAS1 *DECK DROT
3088 subroutine DROT(N,DX,INCX,DY,INCY,C,S)
3090 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3091 integer :: N,INCX,INCY
3092 real(kind=8),DIMENSION(1) :: DX,DY
3095 ! APPLIES A PLANE ROTATION.
3096 ! DX(I) = C*DX(I) + S*DY(I)
3097 ! DY(I) = -S*DX(I) + C*DY(I)
3098 ! JACK DONGARRA, LINPACK, 3/11/78.
3101 real(kind=8) :: DTEMP
3103 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
3105 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
3110 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
3111 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
3113 DTEMP = C*DX(IX) + S*DY(IY)
3114 DY(IY) = C*DY(IY) - S*DX(IX)
3121 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
3124 DTEMP = C*DX(I) + S*DY(I)
3125 DY(I) = C*DY(I) - S*DX(I)
3130 !-----------------------------------------------------------------------------
3131 !*MODULE BLAS1 *DECK DROTG
3132 subroutine DROTG(DA,DB,C,S)
3134 ! CONSTRUCT GIVENS PLANE ROTATION.
3135 ! JACK DONGARRA, LINPACK, 3/11/78.
3137 real(kind=8) :: DA,DB,C,S,ROE,SCALE,R,Z
3139 real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00
3141 !-----------------------------------------------------------------------
3145 IF( ABS(DA) .GT. ABS(DB) ) ROE = DA
3146 SCALE = ABS(DA) + ABS(DB)
3147 IF( SCALE .NE. ZERO ) GO TO 10
3153 10 R = SCALE*SQRT((DA/SCALE)**2 + (DB/SCALE)**2)
3158 IF( ABS(DA) .GT. ABS(DB) ) Z = S
3159 IF( ABS(DB) .GE. ABS(DA) .AND. C .NE. ZERO ) Z = ONE/C
3163 end subroutine DROTG
3164 !-----------------------------------------------------------------------------
3165 !*MODULE BLAS1 *DECK DSCAL
3166 subroutine DSCAL(N,DA,DX,INCX)
3168 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3170 real(kind=8),DIMENSION(1) :: DX
3173 ! SCALES A VECTOR BY A CONSTANT.
3174 ! DX(I) = DA * DX(I)
3175 ! USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
3176 ! JACK DONGARRA, LINPACK, 3/11/78.
3178 integer :: NINCX,m,mp1,i
3180 IF(INCX.EQ.1)GO TO 20
3182 ! CODE FOR INCREMENT NOT EQUAL TO 1
3185 DO 10 I = 1,NINCX,INCX
3190 ! CODE FOR INCREMENT EQUAL TO 1
3196 IF( M .EQ. 0 ) GO TO 40
3200 IF( N .LT. 5 ) RETURN
3204 DX(I + 1) = DA*DX(I + 1)
3205 DX(I + 2) = DA*DX(I + 2)
3206 DX(I + 3) = DA*DX(I + 3)
3207 DX(I + 4) = DA*DX(I + 4)
3210 end subroutine DSCAL
3211 !-----------------------------------------------------------------------------
3212 !*MODULE BLAS1 *DECK DSWAP
3213 subroutine DSWAP(N,DX,INCX,DY,INCY)
3215 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3216 integer :: N,INCX,INCY
3217 real(kind=8),DIMENSION(1) :: DX,DY
3219 ! INTERCHANGES TWO VECTORS.
3221 ! USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
3222 ! JACK DONGARRA, LINPACK, 3/11/78.
3224 integer :: ix,iy,i,m,mp1
3225 real(kind=8) :: DTEMP
3227 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
3229 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
3234 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
3235 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
3245 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
3251 IF( M .EQ. 0 ) GO TO 40
3257 IF( N .LT. 3 ) RETURN
3264 DX(I + 1) = DY(I + 1)
3267 DX(I + 2) = DY(I + 2)
3271 end subroutine DSWAP
3272 !-----------------------------------------------------------------------------
3273 !*MODULE BLAS1 *DECK IDAMAX
3274 integer function IDAMAX(N,DX,INCX)
3276 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3278 real(kind=8),DIMENSION(1) :: DX
3280 ! FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
3281 ! JACK DONGARRA, LINPACK, 3/11/78.
3284 real(kind=8) :: RMAX
3286 IF( N .LT. 1 ) RETURN
3289 IF(INCX.EQ.1)GO TO 20
3291 ! CODE FOR INCREMENT NOT EQUAL TO 1
3297 IF(ABS(DX(IX)).LE.RMAX) GO TO 5
3304 ! CODE FOR INCREMENT EQUAL TO 1
3306 20 RMAX = ABS(DX(1))
3308 IF(ABS(DX(I)).LE.RMAX) GO TO 30
3314 !-----------------------------------------------------------------------------
3315 !*MODULE BLAS *DECK DGEMV
3316 subroutine DGEMV(FORMA,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
3318 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3319 integer :: M,N,INCX,INCY,LDA
3320 CHARACTER(len=1) :: FORMA
3321 real(kind=8),DIMENSION(LDA,*) :: A
3322 real(kind=8),DIMENSION(*) :: X,Y
3323 real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00
3324 real(kind=8) :: ALPHA,BETA
3327 ! CLONE OF -DGEMV- WRITTEN BY MIKE SCHMIDT
3330 IF(FORMA.EQ.'T') GO TO 200
3332 ! Y = ALPHA * A * X + BETA * Y
3334 IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN
3336 Y(LOCY) = DDOT(N,A(I,1),LDA,X,INCX)
3341 Y(LOCY) = ALPHA*DDOT(N,A(I,1),LDA,X,INCX) + BETA*Y(LOCY)
3347 ! Y = ALPHA * A-TRANSPOSE * X + BETA * Y
3350 IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN
3352 Y(LOCY) = DDOT(M,A(1,I),1,X,INCX)
3357 Y(LOCY) = ALPHA*DDOT(M,A(1,I),1,X,INCX) + BETA*Y(LOCY)
3362 end subroutine DGEMV
3363 !-----------------------------------------------------------------------------
3364 !-----------------------------------------------------------------------------