2 !-----------------------------------------------------------------------------
4 use MD_data, only:D_ban,IP
8 !-----------------------------------------------------------------------------
10 !-----------------------------------------------------------------------------
12 !-----------------------------------------------------------------------------
16 !-----------------------------------------------------------------------------
17 !*MODULE MTHLIB *DECK VCLR
18 subroutine VCLR(A,INCA,N)
20 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
22 real(kind=8),DIMENSION(*) :: A
24 real(kind=8),PARAMETER :: ZERO=0.0D+00
28 ! ----- ZERO OUT VECTOR -A-, USING INCREMENT -INCA- -----
30 IF (INCA .NE. 1) GO TO 200
44 !-----------------------------------------------------------------------------
46 !-----------------------------------------------------------------------------
47 subroutine BANACH(N,NMAX,A,X,osob)
48 !**********************
50 ! implicit real*8 (a-h,o-z)
51 ! include 'DIMENSIONS'
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
58 real(kind=8) :: xx,aij,aijd
61 !el allocate(D_ban(6*nres))
64 if (dabs(a(1,1)).lt.1.0d-15) then
86 if (dabs(xx).lt.1.0d-15) then
93 CALL BANAII(N,NMAX,A,X)
96 !-----------------------------------------------------------------------------
97 subroutine BANAII(N,NMAX,A,X)
98 !************************
99 ! implicit real*8 (a-h,o-z)
100 ! include 'DIMENSIONS'
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
125 end subroutine BANAII
126 !-----------------------------------------------------------------------------
127 subroutine MATINVERT(N,NMAX,A,A1,osob)
129 ! implicit real*8 (a-h,o-z)
130 ! include 'DIMENSIONS'
132 real(kind=8),DIMENSION(NMAX,NMAX) :: A,A1
133 !el real(kind=8),DIMENSION(6*nres) :: D !(MAXRES6) maxres6=6*maxres
135 real(kind=8),DIMENSION(NMAX) :: X
142 CALL BANACH(N,NMAX,A,X,osob)
152 CALL BANAII(N,NMAX,A,X)
158 end subroutine MATINVERT
159 !-----------------------------------------------------------------------------
161 !-----------------------------------------------------------------------------
162 subroutine bond_move(nbond,nstart,psi,lprint,error)
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
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
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'
188 write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
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.'
197 ! Generate the reference system.
201 call refsys(i2,i3,i4,e1,e2,e3,error)
202 ! Return, if couldn't define the reference system.
204 ! Compute the transformation matrix.
222 if (print_mc.gt.2) then
223 write (iout,'(a)') 'Reference system and matrix r:'
225 write(iout,'(i5,2(3f10.5,5x))')i,(e(i,j),j=1,3),(rot(i,j),j=1,3)
229 call matmult(rot,e,trans)
237 call matmult(e,trans,trans)
240 write (iout,'(a)') 'The trans matrix:'
242 write (iout,'(i5,3f10.5)') i,(trans(i,j),j=1,3)
250 rij=rij+trans(j,k)*(c(k,i)-c(k,nstart))
260 write (iout,'(a)') 'Rotated coordinates of the fragment'
262 write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
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)
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)
279 if (print_mc.gt.2) then
280 write (iout,'(/a,i3,a,i3,a/)') &
281 'Moved internal coordinates of the ',nstart,'-',nend,&
283 do i=nstart+1,nstart+2
284 write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
287 write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
291 end subroutine bond_move
292 !-----------------------------------------------------------------------------
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
318 !*MODULE EIGEN *DECK EINVIT
319 subroutine EINVIT(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
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)
330 !* THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
331 !* SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES.
334 !* INVERSE ITERATION.
338 !* MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
339 !* ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
340 !* DIMENSION STATEMENT.
343 !* CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
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.
359 !* THE NUMBER OF SPECIFIED EIGENVECTORS.
361 !* CONTAINS THE M EIGENVALUES IN ASCENDING
362 !* OR DESCENDING ORDER.
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
369 !* IERR - INTEGER (LOGICAL UNIT NUMBER)
370 !* LOGICAL UNIT FOR ERROR MESSAGES
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.
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)
385 !* RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
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
398 !* DIFFERENCES FROM EISPACK 3 -
399 !* EPS3 IS SCALED BY EPSCAL (ENHANCES CONVERGENCE, BUT
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
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
415 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
416 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
420 LOGICAL :: CONVGD !el,GOPARR,DSKWRK,MASWRK
422 INTEGER :: GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG
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,
433 !el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
434 !el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
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
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)')
444 !-----------------------------------------------------------------------
461 ! .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
464 IF (E2(Q+1) .EQ. ZERO) GO TO 140
468 ! .......... FIND VECTORS BY INVERSE ITERATION ..........
476 IF (IND(R) .NE. TAG) GO TO 920
479 IF (S .NE. 0) GO TO 510
481 ! .......... CHECK FOR ISOLATED ROOT ..........
492 NORM = MAX( NORM, ABS(D(I)) + ABS(E(I)) )
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 .........
501 EPS3 = EPSCAL * EPSLON(NORM)
509 ! .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
511 510 IF (ABS(X1-X0) .GE. EPS2) THEN
521 IF (ORDER * (X1 - X0) .LE. EPS3) X1 = X0 + ORDER * EPS3
524 ! .......... ELIMINATION WITH INTERCHANGES AND
525 ! INITIALIZATION OF VECTOR ..........
534 IF (ABS(E(I)) .GT. ABS(U)) THEN
536 ! EXCHANGE ROWS BEFORE ELIMINATION
538 ! *** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
539 ! E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......
546 U = V - XU * RV2(I-1)
551 ! STRAIGHT ELIMINATION
558 U = D(I) - X1 - XU * V
563 IF (ABS(U) .LE. EPS3) U = EPS3
568 ! DO INVERSE ITERATIONS
572 IF (ITS .EQ. 1) GO TO 600
574 ! .......... FORWARD SUBSTITUTION ..........
576 IF (NORM .EQ. ZERO) THEN
582 CALL DSCAL (Q-P+1, XU, RV6(P), 1)
585 ! ... ELIMINATION OPERATIONS ON NEXT VECTOR
590 ! IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
591 ! WAS PERFORMED EARLIER IN THE
592 ! TRIANGULARIZATION PROCESS ..........
594 IF (RV1(I-1) .EQ. E(I)) THEN
600 RV6(I) = U - RV4(I) * RV6(I-1)
604 ! .......... BACK SUBSTITUTION
606 RV6(Q) = RV6(Q) / RV1(Q)
610 DO 620 I = Q-1, P, -1
611 RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
616 IF (GROUP .EQ. 0) GO TO 700
618 ! ....... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
619 ! MEMBERS OF GROUP ..........
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),&
628 NORM = DASUM(Q-P+1, RV6(P), 1)
631 IF (CONVGD) GO TO 840
632 IF (NORM .GE. ONE) CONVGD = .TRUE.
635 ! .......... NORMALIZE SO THAT SUM OF SQUARES IS
636 ! 1 AND EXPAND TO FULL ORDER ..........
640 XU = ONE / DNRM2(Q-P+1,RV6(P),1)
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
658 ! *** SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
660 IF (RHO .GT. HUNDRD) IERR = -R
666 IF (Q .EQ. N) GO TO 940
670 end subroutine EINVIT
671 !-----------------------------------------------------------------------------
672 !*MODULE EIGEN *DECK ELAUM
673 subroutine ELAU(HINV,L,D,A,E)
675 integer :: L,JL,JK,J,JM1,K
678 real(kind=8) :: E(*)!el E(L)
684 real(kind=8),PARAMETER :: ZERO = 0.0D+00, HALF = 0.5D+00
696 E(K) = E(K) + A(JK) * F
704 ! .......... FORM P ..........
712 ! .......... FORM Q ..........
716 250 E(J) = E(J) - HH * D(J)
720 !-----------------------------------------------------------------------------
721 !*MODULE EIGEN *DECK EPSLON
722 real(kind=8) function EPSLON(X)
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
729 !* ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
732 !* X - WORKING PRECISION REAL
733 !* VALUES TO FIND EPSLON FOR
736 !* EPSLON - WORKING PRECISION REAL
737 !* SMALLEST POSITIVE VALUE SUCH THAT X+EPSLON .NE. ZERO
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
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.
759 !* DIFFERENCES FROM EISPACK 3 -
760 !* USE IS MADE OF PARAMETER STATEMENTS AND INTRINSIC FUNCTIONS
761 !* --NO EXECUTEABLE CODE CHANGES--
764 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
765 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
767 real(kind=8) :: A,B,C,EPS,X
769 real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, FOUR=4.0D+00
771 !-----------------------------------------------------------------------
777 IF (EPS .EQ. ZERO) GO TO 10
781 !-----------------------------------------------------------------------------
782 !*MODULE EIGEN *DECK EQLRAT
783 subroutine EQLRAT(N,DIAG,E,E2IN,D,IND,IERR,E2)
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)
793 !* THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
794 !* TRIDIAGONAL MATRIX
801 !* THE ORDER OF THE MATRIX.
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.
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)
819 !* ZERO FOR NORMAL RETURN,
820 !* J IF THE J-TH EIGENVALUE HAS NOT BEEN
821 !* DETERMINED AFTER 30 ITERATIONS.
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
830 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
831 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
833 INTEGER :: I,J,L,M,N,II,L1,IERR
834 INTEGER,dimension(N) :: IND
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
839 real(kind=8),PARAMETER :: ZERO = 0.0D+00, SCALE= 1.0D+00/64.0D+00, ONE = 1.0D+00
842 !-----------------------------------------------------------------------
848 IF (N .EQ. 1) GO TO 1001
852 100 E2(I-1) = E2IN(I)
862 H = ABS(D(L)) + ABS(E(L))
863 IF (T .GE. H) GO TO 105
869 ! .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
872 IF (E2(M) .GT. C) GO TO 110
873 ! .......... E2(N) IS ALWAYS ZERO, SO THERE IS AN EXIT
874 ! FROM THE LOOP ..........
876 IF (M .LE. K) GO TO 125
877 IF (M .NE. N) E2IN(M+1) = ZERO
881 IF (M .EQ. L) GO TO 210
886 ! .......... FORM SHIFT ..........
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))
899 ! .......... RATIONAL QL TRANSFORMATION ..........
908 D(I+1) = H + S * (H + D(I))
909 G = D(I) - E2(I) / G + B
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
919 IF (E2(L) .EQ. ZERO) GO TO 210
921 ! .......... SET ERROR -- NO CONVERGENCE TO AN
922 ! EIGENVALUE AFTER 30 ITERATIONS ..........
929 ! .......... ORDER EIGENVALUES ..........
931 IF (L .EQ. 1) GO TO 250
932 IF (P .LT. D(1)) GO TO 230
934 ! .......... LOOP TO FIND ORDERED POSITION
936 IF (P .LT. D(I)) GO TO 220
939 IF (I .EQ. L) GO TO 250
941 DO 240 II = L, I+1, -1
952 end subroutine EQLRAT
953 !-----------------------------------------------------------------------------
954 !*MODULE EIGEN *DECK ESTPI1
955 real(kind=8) function ESTPI1(N,EVAL,D,E,X,ANORM)
958 !* STEPHEN T. ELBERT (AMES LABORATORY-USDOE) DATE: 5 DEC 1986
961 !* EVALUATE SYMMETRIC TRIDIAGONAL MATRIX PERFORMANCE INDEX
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.
976 !* THE ORDER OF THE MATRIX A.
978 !* THE EIGENVALUE CORRESPONDING TO VECTOR X.
980 !* THE DIAGONAL VECTOR OF A.
982 !* THE SUB-DIAGONAL VECTOR OF A.
984 !* AN EIGENVECTOR OF A.
986 !* THE NORM OF A IF IT HAS BEEN PREVIOUSLY COMPUTED.
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
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)
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
1007 real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00
1009 !-----------------------------------------------------------------------
1012 IF( N .LE. 1 ) RETURN
1014 IF (ANORM .EQ. ZERO) THEN
1018 ANORM = MAX( ABS(D(1)) + ABS(E(2)), &
1019 ABS(D(N)) + ABS(E(N)))
1021 ANORM = MAX( ANORM, ABS(E(I))+ABS(D(I))+ABS(E(I+1)))
1023 IF(ANORM .EQ. ZERO) ANORM = ONE
1026 ! COMPUTE NORMS OF RESIDUAL AND EIGENVECTOR
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))
1032 XNORM = XNORM + ABS(X(I))
1033 RNORM = RNORM + ABS(E(I)*X(I-1) + (D(I)-EVAL)*X(I) &
1037 ESTPI1 = RNORM / (EPSLON(SIZE)*ANORM*XNORM)
1040 !-----------------------------------------------------------------------------
1041 !*MODULE EIGEN *DECK ETRBK3
1042 subroutine ETRBK3(NM,N,NV,A,M,Z)
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)
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.
1058 !* THE CALCULATION IS CARRIED OUT BY FORMING THE MATRIX PRODUCT
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.
1075 !* MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
1076 !* ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
1077 !* DIMENSION STATEMENT.
1079 !* THE ORDER OF THE MATRIX A.
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.
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.
1094 !* Z - W.P. REAL (NM,M)
1095 !* CONTAINS THE TRANSFORMED EIGENVECTORS
1096 !* IN ITS FIRST M COLUMNS.
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.
1105 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
1106 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
1108 INTEGER :: I,II,IM1,IZ,J,M,N,NM,NV
1110 real(kind=8) :: A(NV),Z(NM,M)
1111 real(kind=8) :: H,S !,DDOT
1113 real(kind=8),PARAMETER :: ZERO = 0.0D+00
1115 !-----------------------------------------------------------------------
1117 IF (M .EQ. 0) RETURN
1118 IF (N .LE. 2) RETURN
1125 IF (H .EQ. ZERO) GO TO 140
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)
1133 end subroutine ETRBK3
1134 !-----------------------------------------------------------------------------
1135 !*MODULE EIGEN *DECK ETRED3
1136 subroutine ETRED3(N,NV,A,D,E,E2)
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
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.
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
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.
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.
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.
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.
1189 !* THE ORDER OF THE MATRIX.
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.
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
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.
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
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
1224 !* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
1225 !* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
1227 INTEGER :: I,IIA,IZ0,L,N,NV
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
1234 real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00
1236 !-----------------------------------------------------------------------
1238 IF (N .LE. 2) GO TO 310
1240 AIIMAX = ABS(A(IZ0))
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
1249 ! THIS ROW IS ALREADY IN TRI-DIAGONAL FORM
1252 IF (AIIMAX+D(I) .EQ. AIIMAX) D(I) = ZERO
1254 IF (AIIMAX+E(I) .EQ. AIIMAX) E(I) = ZERO
1261 SCALEI = ONE / SCALE
1262 CALL DSCAL(L,SCALEI,A(IZ0+1),1)
1263 HROOT = DNRM2(L,A(IZ0+1),1)
1269 H = HROOT*HROOT - F * G
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)
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)
1294 !* AUTHOR: S. T. ELBERT, AMES LABORATORY-USDOE, JUNE 1985
1297 !* FINDS (ALL) EIGENVALUES AND (SOME OR ALL) EIGENVECTORS
1299 !* OF A REAL SYMMETRIC PACKED MATRIX.
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
1313 !* THESE ROUTINES ARE MODIFICATIONS OF THE EISPACK 3
1314 !* ROUTINES TRED3, TQLRAT, TINVIT AND TRBAK3
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)
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.
1328 !* ORDER OF MATRIX A.
1330 !* NUMBER OF VECTORS DESIRED. 0 .LE. NVECT .LE. N.
1332 !* DIMENSION OF A IN CALLING ROUTINE. MUST NOT BE LESS
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.
1345 !* ROOT ORDERING FLAG.
1346 !* = 0, ROOTS WILL BE PUT IN ASCENDING ORDER.
1347 !* = 2, ROOTS WILL BE PUT IN DESCENDING ORDER.
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).
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.)
1365 !el LOGICAL :: GOPARR,DSKWRK,MASWRK
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)
1372 real(kind=8) :: VECT(NV,*)
1376 !el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
1377 real(kind=8) :: DSKW
1379 !el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
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')
1396 integer :: LMSGFL,I,J,L,JSV,KLIM,K
1398 !-----------------------------------------------------------------------
1401 IF (MSGFL .EQ. 0) LMSGFL=6
1403 IF (N .LE. 0) GO TO 800
1405 IF ( (N*N+N)/2 .GT. LENA) GO TO 810
1407 ! REDUCE REAL SYMMETRIC MATRIX A TO TRIDIAGONAL FORM
1409 CALL ETRED3(N,LENA,A,B(1,1),B(1,2),B(1,3))
1411 ! FIND ALL EIGENVALUES OF TRIDIAGONAL MATRIX
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
1416 ! CHECK THE DESIRED ORDER OF THE EIGENVALUES
1419 IF (IORDER .EQ. 0) GO TO 300
1420 IF (IORDER .NE. 2) GO TO 850
1422 ! ORDER ROOTS IN DESCENDING ORDER (LARGEST FIRST)...
1423 ! TURN ROOT AND IND ARRAYS END FOR END
1435 ! FIND I AND J MARKING THE START AND END OF A SEQUENCE
1436 ! OF DEGENERATE ROOTS
1441 IF (I .GT. N) GO TO 300
1443 IF (ROOT(J) .NE. ROOT(I)) GO TO 240
1448 IF (J .EQ. I) GO TO 220
1450 ! TURN AROUND IND BETWEEN I AND J
1466 IF (NVECT .LE. 0) RETURN
1467 IF (NV .LT. N) GO TO 830
1469 ! FIND EIGENVECTORS OF TRI-DIAGONAL MATRIX VIA INVERSE ITERATION
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
1476 ! FIND EIGENVECTORS OF SYMMETRIC MATRIX VIA BACK TRANSFORMATION
1479 CALL ETRBK3(NV,N,LENA,A,NVECT,VECT)
1482 ! ERROR MESSAGE SECTION
1484 800 IF (LMSGFL .LT. 0) RETURN
1485 IF (MASWRK) WRITE(LMSGFL,906)
1488 810 IF (LMSGFL .LT. 0) RETURN
1489 IF (MASWRK) WRITE(LMSGFL,901)
1492 820 IF (LMSGFL .LT. 0) RETURN
1493 IF (MASWRK) WRITE(LMSGFL,902) IERR
1496 830 IF (LMSGFL .LT. 0) RETURN
1497 IF (MASWRK) WRITE(LMSGFL,903)
1501 IF ((LMSGFL .GT. 0).AND.MASWRK) WRITE(LMSGFL,904) -IERR
1505 IF (LMSGFL .LT. 0) RETURN
1506 IF (MASWRK) WRITE(LMSGFL,905)
1510 IF (MASWRK) WRITE(LMSGFL,900) N,NVECT,LENA,NV,IORDER,IERR
1512 end subroutine EVVRSP
1513 !-----------------------------------------------------------------------------
1514 !*MODULE EIGEN *DECK FREDA
1515 subroutine FREDA(L,D,A,E)
1518 real(kind=8) :: A(*)
1519 real(kind=8) :: D(L)
1520 real(kind=8) :: E(*)!el E(L)
1526 ! .......... FORM REDUCED A ..........
1533 A(JK) = A(JK) - F * E(K) - G * D(K)
1539 end subroutine FREDA
1540 !-----------------------------------------------------------------------------
1541 !*MODULE EIGEN *DECK GIVEIS
1542 subroutine GIVEIS(N,NVECT,NV,A,B,INDB,ROOT,VECT,IERR)
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
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.
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
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.)
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.
1582 ! IF TINVIT FAILS TO CONVERGE, TQL2 IS CALLED
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.
1589 !el DATA ONE, ZERO /1.0D+00, 0.0D+00/
1590 real(kind=8) :: ZERO = 0.0D+00, ONE = 1.0D+00
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
1598 ! TO REORDER ROOTS...
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
1613 ! IF INVERSE ITERATION GIVES AN ERROR IN DETERMINING THE
1614 ! EIGENVECTORS, TRY THE QL ALGORITHM IF ALL THE EIGENVECTORS
1617 IF (NVECT .NE. N) RETURN
1624 CALL TQL2 (NV,N,B(1,1),B(1,2),VECT,IERR)
1628 IF (IERR .NE. 0) RETURN
1629 160 CALL TRBK3B(NV,N,(N*N+N)/2,A,NVECT,VECT)
1631 end subroutine GIVEIS
1632 !-----------------------------------------------------------------------------
1633 !*MODULE EIGEN *DECK GLDIAG
1634 subroutine GLDIAG(LDVECT,NVECT,N,H,WRK,EIG,VECTOR,IERR,IWRK)
1636 ! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1641 !el LOGICAL :: GOPARR,DSKWRK,MASWRK
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
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
1657 integer :: LENH,KORDER
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
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
1679 ! ----- USE STEVE ELBERT'S ROUTINE -----
1681 IF(KDIAG.LE.1 .OR. KDIAG.GT.3) THEN
1684 CALL EVVRSP(IW,N,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR, &
1688 ! ----- USE MODIFIED EISPAK ROUTINE -----
1691 CALL GIVEIS(N,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
1693 ! ----- USE JACOBI ROTATION ROUTINE -----
1697 CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,N)
1699 IF (MASWRK) WRITE(IW,9000) N,NVECT,LDVECT
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)
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
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).
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.
1733 ! N IS THE ORDER OF THE MATRIX,
1735 ! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
1737 ! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
1738 ! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
1740 ! E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
1741 ! E2(1) IS ARBITRARY.
1745 ! D AND E ARE UNALTERED,
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,
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,
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.,
1763 ! ZERO FOR NORMAL RETURN,
1764 ! J IF THE J-TH EIGENVALUE HAS NOT BEEN
1765 ! DETERMINED AFTER 30 ITERATIONS,
1767 ! RV1 IS A TEMPORARY STORAGE ARRAY.
1769 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
1770 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
1772 ! ------------------------------------------------------------------
1774 ! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
1775 ! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
1778 MACHEP = 2.0D+00**(-50)
1786 IF (I .NE. 1) RV1(I-1) = E(I)
1794 ! ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
1796 IF (M .EQ. N) GO TO 160
1797 IF (ABS(RV1(M)) .LE. MACHEP * (ABS(W(M)) + ABS(W(M+1)))) GO TO &
1799 ! ********** GUARD AGAINST UNDERFLOWED ELEMENT OF E2 **********
1800 IF (E2(M+1) .EQ. 0.0D+00) GO TO 180
1803 160 IF (M .LE. K) GO TO 200
1804 IF (M .NE. N) E2(M+1) = 0.0D+00
1808 IF (M .EQ. L) GO TO 280
1809 IF (J .EQ. 30) GO TO 380
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))
1819 ! ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
1824 IF (ABS(F) .LT. ABS(G)) GO TO 220
1826 R = SQRT(C*C+1.0D+00)
1832 R = SQRT(S*S+1.0D+00)
1837 R = (W(I) - G) * S + 2.0D+00 * C * B
1847 ! ********** ORDER EIGENVALUES **********
1848 280 IF (L .EQ. 1) GO TO 320
1849 ! ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
1852 IF (P .GE. W(I-1)) GO TO 340
1863 ! ********** SET ERROR -- NO CONVERGENCE TO AN
1864 ! EIGENVALUE AFTER 30 ITERATIONS **********
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)
1873 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
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
1881 real(kind=8),PARAMETER :: ONE=1.0D+00
1882 integer :: i,NB1,NB2,NMIN,NMAX
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.
1890 CALL VCLR(VEC,1,LDVEC*N)
1896 NB2 = (NB1*NB1+NB1)/2
1900 CALL JACDIA(A,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1903 EIG(I) = A((I*I+I)/2)
1906 CALL JACORD(VEC,EIG,NB1,LDVEC)
1908 end subroutine JACDG
1909 !-----------------------------------------------------------------------------
1910 !*MODULE EIGEN *DECK JACDIA
1911 subroutine JACDIA(F,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
1913 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
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
1922 !el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
1923 !el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
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,&
1933 real(kind=8) :: T,TT,EPS,SD,TEST,DIF,CX,SX,T2X2,T2X25,T1,T2
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
1940 ! THE ROTATIONS AMONG THE LAST NB1-NMAX BASIS FUNCTIONS ARE NOT
1946 EPS = 64.0D+00*EPSLON(ONE)
1948 ! LOOP OVER COLUMNS (K) OF TRIANGULAR MATRIX TO DETERMINE
1949 ! LARGEST OFF-DIAGONAL ELEMENTS IN ROW(I).
1954 IF(I.LT.NMIN .OR. I.EQ.1) GO TO 20
1958 IF(ABS(BIG(I)).GE.ABS(F(II+K))) GO TO 10
1964 ! ----- 2X2 JACOBI ITERATIONS BEGIN HERE -----
1966 MAXIT=MAX(NB2*20,500)
1971 ! FIND SMALLEST DIAGONAL ELEMENT
1977 SD= MIN(SD,ABS(F(JJ)))
1979 TEST = MAX(EPS, C2*MAX(SD,C6))
1981 ! FIND LARGEST OFF-DIAGONAL ELEMENT
1987 IF(T.GE.ABS(BIG(I))) GO TO 50
1992 ! TEST FOR CONVERGENCE, THEN DETERMINE ROTATION.
1994 IF(T.LT.TEST) RETURN
1997 IF(ITER.GT.MAXIT) THEN
1999 WRITE(6,*) 'JACOBI DIAGONALIZATION FAILS, DIMENSION=',NB1
2000 WRITE(6,9020) ITER,T,TEST,SD
2009 DIF=F(IAA+IA)-F(IBB+IB)
2010 IF(ABS(DIF).GT.C3*T) GO TO 70
2016 IF(T2X25 .GT. C4) GO TO 80
2020 80 IF(T2X25 .GT. C5) GO TO 90
2021 SX=T2X2*(ONE-D1500*T2X25)
2024 90 IF(T2X25 .GT. C6) GO TO 100
2025 CX=ONE+T2X25*(T2X25*D1375 - D0500)
2026 SX= T2X2*(ONE + T2X25*(T2X25*D3875 - D1500))
2028 100 T=D0250 / SQRT(D0250 + T2X25)
2030 SX= SIGN( SQRT(D0500 - T),T2X2)
2036 F(IEAR)=F(IEAR)*CX+F(IEBR)*SX
2037 F(IEBR)=T-F(IEBR)*CX
2038 IF(IR-IA) 220,120,130
2044 IF(JBIG(IR)) 200,220,200
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
2054 160 IF( ABS(T) .GE. ABS(F(IEBR))) GO TO 170
2055 IF(IB.GT.NMAX) GO TO 170
2059 180 IF( ABS(T) .LT. ABS(BIG(IR))) GO TO 190
2063 190 IF(IA .NE. JBIG(IR) .AND. IB .NE. JBIG(IR)) GO TO 220
2069 IF(ABS(BIG(IR)) .GE. ABS(F(K))) GO TO 210
2077 T1=VEC(I,IA)*CX + VEC(I,IB)*SX
2078 T2=VEC(I,IA)*SX - VEC(I,IB)*CX
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)
2090 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2092 real(kind=8),DIMENSION(LDVEC,N) :: VEC
2093 real(kind=8),DIMENSION(N) :: EIG
2097 ! ---- SORT EIGENDATA INTO ASCENDING ORDER -----
2102 IF (EIG(J) .LT. EIG(JJ)) JJ = J
2104 IF (JJ .EQ. I) GO TO 290
2110 VEC(J,JJ) = VEC(J,I)
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)
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
2132 ! ------------------------------------------------------------------
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).
2138 ! THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
2139 ! SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
2140 ! USING INVERSE ITERATION.
2144 ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2145 ! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2146 ! DIMENSION STATEMENT,
2148 ! N IS THE ORDER OF THE MATRIX,
2150 ! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2152 ! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2153 ! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
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,
2165 ! M IS THE NUMBER OF SPECIFIED EIGENVALUES,
2167 ! W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER,
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.
2176 ! ALL INPUT ARRAYS ARE UNALTERED,
2178 ! Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
2179 ! ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO,
2182 ! ZERO FOR NORMAL RETURN,
2183 ! -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
2184 ! EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS,
2186 ! RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
2188 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2189 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2191 ! ------------------------------------------------------------------
2193 ! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2194 ! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2197 MACHEP = 2.0D+00**(-50)
2200 IF (M .EQ. 0) GO TO 680
2202 ORDER = 1.0D+00 - E2(1)
2212 ! ********** ESTABLISH AND PROCESS NEXT SUBMATRIX **********
2217 IF (Q .EQ. N) GO TO 140
2218 IF (E2(Q+1) .EQ. 0.0D+00) GO TO 140
2220 ! ********** FIND VECTORS BY INVERSE ITERATION **********
2226 IF (IND(R) .NE. TAG) GO TO 660
2229 IF (S .NE. 0) GO TO 220
2230 ! ********** CHECK FOR ISOLATED ROOT **********
2232 IF (P .NE. Q) GO TO 160
2235 160 NORM = ABS(D(P))
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
2247 UK = EPS4 / SQRT(UK)
2251 ! ********** LOOK FOR CLOSE OR COINCIDENT ROOTS **********
2252 220 IF (ABS(X1-X0) .GE. EPS2) GO TO 200
2254 IF (ORDER * (X1 - X0) .LE. 0.0D+00) X1 = X0 + ORDER * EPS3
2255 ! ********** ELIMINATION WITH INTERCHANGES AND
2256 ! INITIALIZATION OF VECTOR **********
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 **********
2268 RV2(I-1) = D(I) - X1
2270 IF (I .NE. Q) RV3(I-1) = E(I+1)
2271 U = V - XU * RV2(I-1)
2279 280 U = D(I) - X1 - XU * V
2280 IF (I .NE. Q) V = E(I+1)
2283 IF (U .EQ. 0.0D+00) U = EPS3
2287 ! ********** BACK SUBSTITUTION
2288 ! FOR I=Q STEP -1 UNTIL P DO -- **********
2289 320 DO 340 II = P, Q
2291 RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
2295 ! ********** ORTHOGONALIZE WITH RESPECT TO PREVIOUS
2296 ! MEMBERS OF GROUP **********
2297 IF (GROUP .EQ. 0) GO TO 400
2300 DO 380 JJ = 1, GROUP
2302 IF (IND(J) .NE. TAG) GO TO 360
2303 XU = DDOT(IQMP,RV6(P),1,Z(P,J),1)
2305 CALL DAXPY(IQMP,-XU,Z(P,J),1,RV6(P),1)
2312 420 NORM = NORM + ABS(RV6(I))
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
2322 440 XU = EPS4 / NORM
2325 460 RV6(I) = RV6(I) * XU
2326 ! ********** ELIMINATION OPERATIONS ON NEXT VECTOR
2327 ! ITERATE **********
2328 480 DO 520 I = IP, Q
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
2336 500 RV6(I) = U - RV4(I) * RV6(I-1)
2341 ! ********** SET ERROR -- NON-CONVERGED EIGENVECTOR **********
2345 ! ********** NORMALIZE SO THAT SUM OF SQUARES IS
2346 ! 1 AND EXPAND TO FULL ORDER **********
2350 RV6(I) = RV6(I) / NORM
2351 580 U = U + RV6(I)**2
2353 XU = 1.0D+00 / SQRT(U)
2356 620 Z(I,R) = 0.0D+00
2359 640 Z(I,R) = RV6(I) * XU
2364 IF (Q .LT. N) GO TO 100
2366 ! ********** LAST CARD OF TINVIT **********
2367 end subroutine TINVTB
2368 !-----------------------------------------------------------------------------
2369 !*MODULE EIGEN *DECK TQL2
2370 subroutine TQL2(NM,N,D,E,Z,IERR)
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
2381 ! THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
2382 ! NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
2384 ! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
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.
2394 ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2395 ! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2396 ! DIMENSION STATEMENT,
2398 ! N IS THE ORDER OF THE MATRIX,
2400 ! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
2402 ! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
2403 ! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
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.
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,
2416 ! E HAS BEEN DESTROYED,
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
2424 ! ZERO FOR NORMAL RETURN,
2425 ! J IF THE J-TH EIGENVALUE HAS NOT BEEN
2426 ! DETERMINED AFTER 30 ITERATIONS.
2428 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2429 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2431 ! ------------------------------------------------------------------
2433 ! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
2434 ! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
2437 MACHEP = 2.0D+00**(-50)
2440 IF (N .EQ. 1) GO TO 400
2451 H = MACHEP * (ABS(D(L)) + ABS(E(L)))
2453 ! ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
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 **********
2460 140 IF (M .EQ. L) GO TO 280
2461 160 IF (J .EQ. 30) GO TO 380
2463 ! ********** FORM SHIFT **********
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))
2475 ! ********** QL TRANSFORMATION **********
2480 ! ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
2485 IF (ABS(P) .LT. ABS(E(I))) GO TO 200
2487 R = SQRT(C*C+1.0D+00)
2493 R = SQRT(C*C+1.0D+00)
2494 E(I+1) = S * E(I) * R
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)
2506 IF (ABS(E(L)) .GT. B) GO TO 160
2509 ! ********** ORDER EIGENVALUES AND EIGENVECTORS **********
2516 IF (D(J) .GE. P) GO TO 320
2521 IF (K .EQ. I) GO TO 360
2525 CALL DSWAP(N,Z(1,I),1,Z(1,K),1)
2530 ! ********** SET ERROR -- NO CONVERGENCE TO AN
2531 ! EIGENVALUE AFTER 30 ITERATIONS **********
2534 ! ********** LAST CARD OF TQL2 **********
2536 !-----------------------------------------------------------------------------
2537 !*MODULE EIGEN *DECK TRBK3B
2538 subroutine TRBK3B(NM,N,NV,A,M,Z)
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
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).
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.
2557 ! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
2558 ! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
2559 ! DIMENSION STATEMENT,
2561 ! N IS THE ORDER OF THE MATRIX,
2563 ! NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2564 ! AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
2566 ! A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS
2567 ! USED IN THE REDUCTION BY TRED3B IN ITS FIRST
2568 ! N*(N+1)/2 POSITIONS,
2570 ! M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
2572 ! Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
2573 ! IN ITS FIRST M COLUMNS.
2577 ! Z CONTAINS THE TRANSFORMED EIGENVECTORS
2578 ! IN ITS FIRST M COLUMNS.
2580 ! NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
2582 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2583 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2585 ! ------------------------------------------------------------------
2587 IF (M .EQ. 0) GO TO 140
2588 IF (N .EQ. 1) GO TO 140
2595 IF (H .EQ. 0.0D+00) GO TO 120
2598 S = -DDOT(L,A(IZ+1),1,Z(1,J),1)
2600 ! ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
2603 CALL DAXPY(L,S,A(IZ+1),1,Z(1,J),1)
2610 ! ********** LAST CARD OF TRBAK3 **********
2611 end subroutine TRBK3B
2612 !-----------------------------------------------------------------------------
2613 !*MODULE EIGEN *DECK TRED3B
2614 subroutine TRED3B(N,NV,A,D,E,E2)
2616 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
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
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).
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.
2634 ! N IS THE ORDER OF THE MATRIX,
2636 ! NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
2637 ! AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
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.
2645 ! A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
2646 ! TRANSFORMATIONS USED IN THE REDUCTION,
2648 ! D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
2650 ! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
2651 ! MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO,
2653 ! E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
2654 ! E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
2656 ! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
2657 ! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
2659 ! ------------------------------------------------------------------
2661 ! ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
2668 IF (L .LT. 1) GO TO 120
2669 ! ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
2673 SCALE = SCALE + ABS(D(K))
2676 IF (SCALE .NE. 0.0D+00) GO TO 140
2686 E2(I) = SCALE * SCALE * H
2688 G = -SIGN(SQRT(H),F)
2692 A(IZ) = SCALE * D(L)
2693 IF (L .EQ. 1) GO TO 280
2701 ! ********** FORM ELEMENT OF A*U **********
2702 IF (JM1 .EQ. 0) GO TO 200
2704 E(K) = E(K) + DT * A(JK)
2705 G = G + D(K) * A(JK)
2708 200 E(J) = G + A(JK) * DT
2710 ! ********** FORM ELEMENT OF P **********
2720 ! ********** FORM REDUCED A **********
2728 A(JK) = A(JK) - F * E(K) - G * D(K)
2732 A(IZ+1) = SCALE * SQRT(H)
2736 ! ********** LAST CARD OF TRED3 **********
2737 end subroutine TRED3B
2738 !-----------------------------------------------------------------------------
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
2748 ! BASIC LINEAR ALGEBRA SUBPROGRAMS (BLAS) FROM LINPACK (LEVEL 1)
2750 ! THIS MODULE SHOULD BE COMPILED ONLY IF SPECIALLY CODED
2751 ! VERSIONS OF THESE ROUTINES ARE NOT AVAILABLE ON THE TARGET MACHINE
2753 !*MODULE BLAS1 *DECK DASUM
2754 real(kind=8) function DASUM(N,DX,INCX)
2756 ! TAKES THE SUM OF THE ABSOLUTE VALUES.
2757 ! JACK DONGARRA, LINPACK, 3/11/78.
2759 real(kind=8) :: DX(1),DTEMP
2760 INTEGER :: I,INCX,M,MP1,N,NINCX
2765 IF(INCX.EQ.1)GO TO 20
2767 ! CODE FOR INCREMENT NOT EQUAL TO 1
2770 DO 10 I = 1,NINCX,INCX
2771 DTEMP = DTEMP + ABS(DX(I))
2776 ! CODE FOR INCREMENT EQUAL TO 1
2782 IF( M .EQ. 0 ) GO TO 40
2784 DTEMP = DTEMP + ABS(DX(I))
2786 IF( N .LT. 6 ) GO TO 60
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))
2795 !-----------------------------------------------------------------------------
2796 !*MODULE BLAS1 *DECK DAXPY
2797 subroutine DAXPY(N,DA,DX,INCX,DY,INCY)
2799 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2800 integer :: N,INCX,INCY
2801 real(kind=8),DIMENSION(1) :: DX,DY
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.
2809 integer :: ix,iy,i,m,mp1
2811 IF (DA .EQ. 0.0D+00) RETURN
2812 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2814 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2819 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2820 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2822 DY(IY) = DY(IY) + DA*DX(IX)
2828 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
2834 IF( M .EQ. 0 ) GO TO 40
2836 DY(I) = DY(I) + DA*DX(I)
2838 IF( N .LT. 4 ) RETURN
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)
2847 end subroutine DAXPY
2848 !-----------------------------------------------------------------------------
2849 !*MODULE BLAS1 *DECK DCOPY
2850 subroutine DCOPY(N,DX,INCX,DY,INCY)
2852 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2853 integer :: N,INCX,INCY
2854 real(kind=8),DIMENSION(*) :: DX,DY
2858 ! USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
2859 ! JACK DONGARRA, LINPACK, 3/11/78.
2861 integer :: ix,iy,m,i,mp1
2863 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2865 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2870 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2871 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2879 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
2885 IF( M .EQ. 0 ) GO TO 40
2889 IF( N .LT. 7 ) RETURN
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)
2901 end subroutine DCOPY
2902 !-----------------------------------------------------------------------------
2903 !*MODULE BLAS1 *DECK DDOT
2904 real(kind=8) function DDOT(N,DX,INCX,DY,INCY)
2906 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
2907 integer :: N,INCX,INCY
2908 real(kind=8),DIMENSION(1) :: DX,DY
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.
2915 integer ::ix,iy,m,mp1,i
2916 real(kind=8) :: DTEMP
2920 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
2922 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
2927 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
2928 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
2930 DTEMP = DTEMP + DX(IX)*DY(IY)
2937 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
2943 IF( M .EQ. 0 ) GO TO 40
2945 DTEMP = DTEMP + DX(I)*DY(I)
2947 IF( N .LT. 5 ) GO TO 60
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)
2956 !-----------------------------------------------------------------------------
2957 !*MODULE BLAS1 *DECK DNRM2
2958 real(kind=8) function DNRM2(N,DX,INCX)
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/
2966 ! EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
2968 ! IF N .LE. 0 RETURN WITH RESULT = 0.
2969 ! IF N .GE. 1 THEN INCX MUST BE .GE. 1
2971 ! C.L.LAWSON, 1978 JAN 08
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.
2978 ! EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
2979 ! U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT)
2980 ! V = LARGEST NO. (OVERFLOW LIMIT)
2982 ! BRIEF OUTLINE OF ALGORITHM..
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.
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 /
3006 IF(N .GT. 0) GO TO 10
3010 10 ASSIGN 30 TO NEXT
3015 20 GO TO NEXT,(30, 50, 70, 110)
3016 30 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
3020 ! PHASE 1. SUM IS ZERO
3022 50 IF( DX(I) .EQ. ZERO) GO TO 200
3023 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
3025 ! PREPARE FOR PHASE 2.
3029 ! PREPARE FOR PHASE 4.
3033 SUM = (SUM / DX(I)) / DX(I)
3034 105 XMAX = ABS(DX(I))
3037 ! PHASE 2. SUM IS SMALL.
3038 ! SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
3040 70 IF( ABS(DX(I)) .GT. CUTLO ) GO TO 75
3042 ! COMMON CODE FOR PHASES 2 AND 4.
3043 ! IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
3045 110 IF( ABS(DX(I)) .LE. XMAX ) GO TO 115
3046 SUM = ONE + SUM * (XMAX / DX(I))**2
3050 115 SUM = SUM + (DX(I)/XMAX)**2
3054 ! PREPARE FOR PHASE 3.
3056 75 SUM = (SUM * XMAX) * XMAX
3059 ! FOR REAL OR D.P. SET HITEST = CUTHI/N
3060 ! FOR COMPLEX SET HITEST = CUTHI/(2*N)
3064 ! PHASE 3. SUM IS MID-RANGE. NO SCALING.
3067 IF(ABS(DX(J)) .GE. HITEST) GO TO 100
3068 95 SUM = SUM + DX(J)**2
3074 IF ( I .LE. NN ) GO TO 20
3078 ! COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
3080 DNRM2 = XMAX * SQRT(SUM)
3084 !-----------------------------------------------------------------------------
3085 !*MODULE BLAS1 *DECK DROT
3086 subroutine DROT(N,DX,INCX,DY,INCY,C,S)
3088 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3089 integer :: N,INCX,INCY
3090 real(kind=8),DIMENSION(1) :: DX,DY
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.
3099 real(kind=8) :: DTEMP
3101 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
3103 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
3108 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
3109 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
3111 DTEMP = C*DX(IX) + S*DY(IY)
3112 DY(IY) = C*DY(IY) - S*DX(IX)
3119 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
3122 DTEMP = C*DX(I) + S*DY(I)
3123 DY(I) = C*DY(I) - S*DX(I)
3128 !-----------------------------------------------------------------------------
3129 !*MODULE BLAS1 *DECK DROTG
3130 subroutine DROTG(DA,DB,C,S)
3132 ! CONSTRUCT GIVENS PLANE ROTATION.
3133 ! JACK DONGARRA, LINPACK, 3/11/78.
3135 real(kind=8) :: DA,DB,C,S,ROE,SCALE,R,Z
3137 real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00
3139 !-----------------------------------------------------------------------
3143 IF( ABS(DA) .GT. ABS(DB) ) ROE = DA
3144 SCALE = ABS(DA) + ABS(DB)
3145 IF( SCALE .NE. ZERO ) GO TO 10
3151 10 R = SCALE*SQRT((DA/SCALE)**2 + (DB/SCALE)**2)
3156 IF( ABS(DA) .GT. ABS(DB) ) Z = S
3157 IF( ABS(DB) .GE. ABS(DA) .AND. C .NE. ZERO ) Z = ONE/C
3161 end subroutine DROTG
3162 !-----------------------------------------------------------------------------
3163 !*MODULE BLAS1 *DECK DSCAL
3164 subroutine DSCAL(N,DA,DX,INCX)
3166 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3168 real(kind=8),DIMENSION(1) :: DX
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.
3176 integer :: NINCX,m,mp1,i
3178 IF(INCX.EQ.1)GO TO 20
3180 ! CODE FOR INCREMENT NOT EQUAL TO 1
3183 DO 10 I = 1,NINCX,INCX
3188 ! CODE FOR INCREMENT EQUAL TO 1
3194 IF( M .EQ. 0 ) GO TO 40
3198 IF( N .LT. 5 ) RETURN
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)
3208 end subroutine DSCAL
3209 !-----------------------------------------------------------------------------
3210 !*MODULE BLAS1 *DECK DSWAP
3211 subroutine DSWAP(N,DX,INCX,DY,INCY)
3213 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3214 integer :: N,INCX,INCY
3215 real(kind=8),DIMENSION(1) :: DX,DY
3217 ! INTERCHANGES TWO VECTORS.
3219 ! USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
3220 ! JACK DONGARRA, LINPACK, 3/11/78.
3222 integer :: ix,iy,i,m,mp1
3223 real(kind=8) :: DTEMP
3225 IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
3227 ! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
3232 IF(INCX.LT.0)IX = (-N+1)*INCX + 1
3233 IF(INCY.LT.0)IY = (-N+1)*INCY + 1
3243 ! CODE FOR BOTH INCREMENTS EQUAL TO 1
3249 IF( M .EQ. 0 ) GO TO 40
3255 IF( N .LT. 3 ) RETURN
3262 DX(I + 1) = DY(I + 1)
3265 DX(I + 2) = DY(I + 2)
3269 end subroutine DSWAP
3270 !-----------------------------------------------------------------------------
3271 !*MODULE BLAS1 *DECK IDAMAX
3272 integer function IDAMAX(N,DX,INCX)
3274 ! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
3276 real(kind=8),DIMENSION(1) :: DX
3278 ! FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
3279 ! JACK DONGARRA, LINPACK, 3/11/78.
3282 real(kind=8) :: RMAX
3284 IF( N .LT. 1 ) RETURN
3287 IF(INCX.EQ.1)GO TO 20
3289 ! CODE FOR INCREMENT NOT EQUAL TO 1
3295 IF(ABS(DX(IX)).LE.RMAX) GO TO 5
3302 ! CODE FOR INCREMENT EQUAL TO 1
3304 20 RMAX = ABS(DX(1))
3306 IF(ABS(DX(I)).LE.RMAX) GO TO 30
3312 !-----------------------------------------------------------------------------
3313 !*MODULE BLAS *DECK DGEMV
3314 subroutine DGEMV(FORMA,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
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
3325 ! CLONE OF -DGEMV- WRITTEN BY MIKE SCHMIDT
3328 IF(FORMA.EQ.'T') GO TO 200
3330 ! Y = ALPHA * A * X + BETA * Y
3332 IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN
3334 Y(LOCY) = DDOT(N,A(I,1),LDA,X,INCX)
3339 Y(LOCY) = ALPHA*DDOT(N,A(I,1),LDA,X,INCX) + BETA*Y(LOCY)
3345 ! Y = ALPHA * A-TRANSPOSE * X + BETA * Y
3348 IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN
3350 Y(LOCY) = DDOT(M,A(1,I),1,X,INCX)
3355 Y(LOCY) = ALPHA*DDOT(M,A(1,I),1,X,INCX) + BETA*Y(LOCY)
3360 end subroutine DGEMV
3361 !-----------------------------------------------------------------------------
3362 !-----------------------------------------------------------------------------