+++ /dev/null
- module MD_calc
-!-----------------------------------------------------------------------------
- use io_units
- use MD_data, only:D_ban,IP
- use geometry_data
-! use prng ! prng.f90 or prng_32.f90
- implicit none
-!
-!-----------------------------------------------------------------------------
- contains
-!-----------------------------------------------------------------------------
-! add.f
-!-----------------------------------------------------------------------------
- subroutine ABRT
- STOP 'IN ABRT'
- end subroutine ABRT
-!-----------------------------------------------------------------------------
-!*MODULE MTHLIB *DECK VCLR
- subroutine VCLR(A,INCA,N)
-!
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
-!
- real(kind=8),DIMENSION(*) :: A
-!
- real(kind=8),PARAMETER :: ZERO=0.0D+00
- integer :: INCA,N
- integer :: l,la
-!
-! ----- ZERO OUT VECTOR -A-, USING INCREMENT -INCA- -----
-!
- IF (INCA .NE. 1) GO TO 200
- DO 110 L=1,N
- A(L) = ZERO
- 110 CONTINUE
- RETURN
-!
- 200 CONTINUE
- LA=1-INCA
- DO 210 L=1,N
- LA=LA+INCA
- A(LA) = ZERO
- 210 CONTINUE
- return
- end subroutine VCLR
-!-----------------------------------------------------------------------------
-! banach.f
-!-----------------------------------------------------------------------------
- subroutine BANACH(N,NMAX,A,X,osob)
-!**********************
-! Banachiewicz
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- integer :: N,NMAX
- real(kind=8),DIMENSION(NMAX,NMAX) :: A
- real(kind=8),DIMENSION(NMAX) :: X
-!el real(kind=8),DIMENSION(6*nres) :: D !(MAXRES6) maxres6=6*maxres
-!el COMMON /BANII/ D
- logical :: osob
- real(kind=8) :: xx,aij,aijd
- integer :: i,j,k,jjjj
-
-!el allocate(D_ban(6*nres))
-
- osob=.false.
- if (dabs(a(1,1)).lt.1.0d-15) then
- osob=.true.
- return
- endif
- D_ban(1)=1./A(1,1)
- DO 80 I=2,N
- A(I,1)=A(1,I)
- DO 81 J=2,I-1
- XX=A(J,I)
- DO 82 K=1,J-1
- XX=XX-A(I,K)*A(J,K)
- 82 CONTINUE
- A(I,J)=XX
- 81 CONTINUE
- XX=A(I,I)
- JJJJ=I-1
- DO 83 J=1,JJJJ
- AIJ=A(I,J)
- AIJD=AIJ*D_ban(J)
- A(I,J)=AIJD
- XX=XX-AIJ*AIJD
- 83 CONTINUE
- if (dabs(xx).lt.1.0d-15) then
- osob=.true.
- return
- endif
- D_ban(I)=1./XX
- 80 CONTINUE
-!
- CALL BANAII(N,NMAX,A,X)
- return
- end subroutine BANACH
-!-----------------------------------------------------------------------------
- subroutine BANAII(N,NMAX,A,X)
-!************************
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- integer :: N,NMAX
- real(kind=8),DIMENSION(NMAX,NMAX) :: A
- real(kind=8),DIMENSION(NMAX) :: X
-!el real(kind=8),DIMENSION(6*nres) :: D !(MAXRES6) maxres6=6*maxres
-!el COMMON /BANII/ D ---> D_ban
- real(kind=8) :: Z
- integer :: i,j,jjjj
- DO 90 I=1,N
- Z=X(I)
- JJJJ=I-1
- DO 91 J=JJJJ,1,-1
- Z=Z-A(I,J)*X(J)
- 91 CONTINUE
- X(I)=Z
- 90 CONTINUE
- DO 92 I=N,1,-1
- Z=X(I)*D_ban(I)
- JJJJ=I+1
- DO 93 J=JJJJ,N
- Z=Z-A(J,I)*X(J)
- 93 CONTINUE
- X(I)=Z
- 92 CONTINUE
- return
- end subroutine BANAII
-!-----------------------------------------------------------------------------
- subroutine MATINVERT(N,NMAX,A,A1,osob)
-
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- integer :: N,NMAX
- real(kind=8),DIMENSION(NMAX,NMAX) :: A,A1
-!el real(kind=8),DIMENSION(6*nres) :: D !(MAXRES6) maxres6=6*maxres
-!el COMMON /BANII/ D
- real(kind=8),DIMENSION(NMAX) :: X
- logical :: osob
- integer :: i,j
- DO I=1,N
- X(I)=0.0
- ENDDO
- X(1)=1.0
- CALL BANACH(N,NMAX,A,X,osob)
- if (osob) return
- DO I=1,N
- A1(I,1)=X(I)
- ENDDO
- DO I=2,N
- DO J=1,N
- X(J)=0.0
- ENDDO
- X(I)=1.0
- CALL BANAII(N,NMAX,A,X)
- DO J=1,N
- A1(J,I)=X(J)
- ENDDO
- ENDDO
- return
- end subroutine MATINVERT
-!-----------------------------------------------------------------------------
-! bond_move.f
-!-----------------------------------------------------------------------------
- subroutine bond_move(nbond,nstart,psi,lprint,error)
-
- use mcm_data, only:print_mc
- use geometry, only:alpha,beta,refsys,matmult
-! Move NBOND fragment starting from the CA(nstart) by angle PSI.
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
- integer :: nbond,nstart
- real(kind=8) :: psi
- logical :: fail,error,lprint
-! include 'COMMON.GEO'
-! include 'COMMON.CHAIN'
-! include 'COMMON.VAR'
-! include 'COMMON.IOUNITS'
-! include 'COMMON.MCM'
- real(kind=8),dimension(3) :: x,e1,e2,e3
- real(kind=8),dimension(3,3) :: e,rot,trans
- real(kind=8) :: cospsi,sinpsi,rij
- integer :: i,j,nend,i2,i3,i4,k
- error=.false.
- nend=nstart+nbond
- if (print_mc.gt.2) then
- write (iout,*) 'nstart=',nstart,' nend=',nend,' nbond=',nbond
- write (iout,*) 'psi=',psi
- write (iout,'(a)') 'Original coordinates of the fragment'
- do i=nstart,nend
- write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
- enddo
- endif
- if (nstart.lt.1 .or. nend .gt.nres .or. nbond.lt.2 .or. &
- nbond.ge.nres-1) then
- write (iout,'(a)') 'Bad data in BOND_MOVE.'
- error=.true.
- return
- endif
-! Generate the reference system.
- i2=nend
- i3=nstart
- i4=nstart+1
- call refsys(i2,i3,i4,e1,e2,e3,error)
-! Return, if couldn't define the reference system.
- if (error) return
-! Compute the transformation matrix.
- cospsi=dcos(psi)
- sinpsi=dsin(psi)
- rot(1,1)=1.0D0
- rot(1,2)=0.0D0
- rot(1,3)=0.0D0
- rot(2,1)=0.0D0
- rot(2,2)=cospsi
- rot(2,3)=-sinpsi
- rot(3,1)=0.0D0
- rot(3,2)=sinpsi
- rot(3,3)=cospsi
- do i=1,3
- e(1,i)=e1(i)
- e(2,i)=e2(i)
- e(3,i)=e3(i)
- enddo
-
- if (print_mc.gt.2) then
- write (iout,'(a)') 'Reference system and matrix r:'
- do i=1,3
- write(iout,'(i5,2(3f10.5,5x))')i,(e(i,j),j=1,3),(rot(i,j),j=1,3)
- enddo
- endif
-
- call matmult(rot,e,trans)
- do i=1,3
- do j=1,3
- e(i,1)=e1(i)
- e(i,2)=e2(i)
- e(i,3)=e3(i)
- enddo
- enddo
- call matmult(e,trans,trans)
-
- if (lprint) then
- write (iout,'(a)') 'The trans matrix:'
- do i=1,3
- write (iout,'(i5,3f10.5)') i,(trans(i,j),j=1,3)
- enddo
- endif
-
- do i=nstart,nend
- do j=1,3
- rij=c(j,nstart)
- do k=1,3
- rij=rij+trans(j,k)*(c(k,i)-c(k,nstart))
- enddo
- x(j)=rij
- enddo
- do j=1,3
- c(j,i)=x(j)
- enddo
- enddo
-
- if (lprint) then
- write (iout,'(a)') 'Rotated coordinates of the fragment'
- do i=nstart,nend
- write (iout,'(i5,3f10.5)') i,(c(j,i),j=1,3)
- enddo
- endif
-
-! call int_from_cart(.false.,lprint)
- if (nstart.gt.1) then
- theta(nstart+1)=alpha(nstart-1,nstart,nstart+1)
- phi(nstart+2)=beta(nstart-1,nstart,nstart+1,nstart+2)
- if (nstart.gt.2) phi(nstart+1)= &
- beta(nstart-2,nstart-1,nstart,nstart+1)
- endif
- if (nend.lt.nres) then
- theta(nend+1)=alpha(nend-1,nend,nend+1)
- phi(nend+1)=beta(nend-2,nend-1,nend,nend+1)
- if (nend.lt.nres-1) phi(nend+2)= &
- beta(nend-1,nend,nend+1,nend+2)
- endif
- if (print_mc.gt.2) then
- write (iout,'(/a,i3,a,i3,a/)') &
- 'Moved internal coordinates of the ',nstart,'-',nend,&
- ' fragment:'
- do i=nstart+1,nstart+2
- write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
- enddo
- do i=nend+1,nend+2
- write (iout,'(i5,2f10.5)') i,rad2deg*theta(i),rad2deg*phi(i)
- enddo
- endif
- return
- end subroutine bond_move
-!-----------------------------------------------------------------------------
-! eigen.f
-!-----------------------------------------------------------------------------
-! 10 AUG 94 - MWS - INCREASE NUMBER OF DAF RECORDS
-! 31 MAR 94 - MWS - ADD A VARIABLE TO END OF MACHSW COMMON
-! 26 JUN 93 - MWS - ETRED3: ADD RETURN FOR SPECIAL CASE N=1
-! 4 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
-! 30 AUG 91 - MWS - JACDIA: LIMIT ITERATIONS, USE EPSLON IN TEST.
-! 14 JUL 91 - MWS - JACOBI DIAGONALIZATION ALLOWS FOR LDVEC.NE.N
-! 29 JAN 91 - TLW - GLDIAG: CHANGED COMMON DIAGSW TO MACHSW
-! 29 OCT 90 - STE - FIX JACDIA UNDEFINED VARIABLE BUG
-! 14 SEP 90 - MK - NEW JACOBI DIAGONALIZATION (KDIAG=3)
-! 27 MAR 88 - MWS - ALLOW FOR VECTOR ROUTINE IN GLDIAG
-! 11 AUG 87 - MWS - SANITIZE CONSTANTS IN EQLRAT
-! 15 FEB 87 - STE - FIX EINVIT SUB-MATRIX LOOP LIMIT
-! SCRATCH ARRAYS ARE N*8 REAL AND N INTEGER
-! 8 DEC 86 - STE - USE PERF INDEX FROM ESTPI1 TO JUDGE EINVIT FAILURE
-! 30 NOV 86 - STE - DELETE LIGENB, MAKE EVVRSP DEFAULT
-! (GIVEIS FAILS ON CRAY FOR BENCHMC AND BENCHCI)
-! 7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
-! 11 OCT 85 - STE - LIGENB,TQL2: USE DROT,DSWAP; TINVTB: SCALE VECTOR
-! BEFORE NORMALIZING; GENERIC FUNCTIONS
-! 24 FEB 84 - STE - INITIALIZE INDEX ARRAY FOR LIGENB IN GLDIAG
-! 1 DEC 83 - STE - CHANGE MACHEP FROM 2**-54 TO 2**-50
-! 28 SEP 82 - MWS - CONVERT TO IBM
-!
-!*MODULE EIGEN *DECK EINVIT
- subroutine EINVIT(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
-!*
-!* AUTHORS-
-!* THIS IS A MODIFICATION OF TINVIT FROM EISPACK EDITION 3
-!* DATED AUGUST 1983.
-!* TINVIT IS A TRANSLATION OF THE INVERSE ITERATION TECHNIQUE
-!* IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
-!* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
-!* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
-!*
-!* PURPOSE -
-!* THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
-!* SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES.
-!*
-!* METHOD -
-!* INVERSE ITERATION.
-!*
-!* ON ENTRY -
-!* NM - INTEGER
-!* MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
-!* ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
-!* DIMENSION STATEMENT.
-!* N - INTEGER
-!* D - W.P. REAL (N)
-!* CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
-!* E - W.P. REAL (N)
-!* CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
-!* IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY.
-!* E2 - W.P. REAL (N)
-!* CONTAINS THE SQUARES OF CORRESPONDING ELEMENTS OF E,
-!* WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
-!* E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
-!* THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE
-!* SUM OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST
-!* CONTAIN 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER,
-!* OR 2.0 IF THE EIGENVALUES ARE IN DESCENDING ORDER.
-!* IF TQLRAT, BISECT, TRIDIB, OR IMTQLV
-!* HAS BEEN USED TO FIND THE EIGENVALUES, THEIR
-!* OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE.
-!* M - INTEGER
-!* THE NUMBER OF SPECIFIED EIGENVECTORS.
-!* W - W.P. REAL (M)
-!* CONTAINS THE M EIGENVALUES IN ASCENDING
-!* OR DESCENDING ORDER.
-!* IND - INTEGER (M)
-!* CONTAINS IN FIRST M POSITIONS THE SUBMATRIX INDICES
-!* ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
-!* 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX
-!* FROM THE TOP, 2 FOR THOSE BELONGING TO THE SECOND
-!* SUBMATRIX, ETC.
-!* IERR - INTEGER (LOGICAL UNIT NUMBER)
-!* LOGICAL UNIT FOR ERROR MESSAGES
-!*
-!* ON EXIT -
-!* ALL INPUT ARRAYS ARE UNALTERED.
-!* Z - W.P. REAL (NM,M)
-!* CONTAINS THE ASSOCIATED SET OF ORTHONORMAL
-!* EIGENVECTORS. ANY VECTOR WHICH WHICH FAILS TO CONVERGE
-!* IS LEFT AS IS (BUT NORMALIZED) WHEN ITERATING STOPPED.
-!* IERR - INTEGER
-!* SET TO
-!* ZERO FOR NORMAL RETURN,
-!* -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
-!* EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS.
-!* (ONLY LAST FAILURE TO CONVERGE IS REPORTED)
-!*
-!* RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
-!*
-!* RV1 - W.P. REAL (N)
-!* DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
-!* RV2 - W.P. REAL (N)
-!* SUPER(1)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
-!* RV3 - W.P. REAL (N)
-!* SUPER(2)-DIAGONAL ELEMENTS OF U FROM LU DECOMPOSITION
-!* RV4 - W.P. REAL (N)
-!* ELEMENTS DEFINING L IN LU DECOMPOSITION
-!* RV6 - W.P. REAL (N)
-!* APPROXIMATE EIGENVECTOR
-!*
-!* DIFFERENCES FROM EISPACK 3 -
-!* EPS3 IS SCALED BY EPSCAL (ENHANCES CONVERGENCE, BUT
-!* LOWERS ACCURACY)!
-!* ONE MORE ITERATION (MINIMUM 2) IS PERFORMED AFTER CONVERGENCE
-!* (ENHANCES ACCURACY)!
-!* REPLACE LOOP WITH PYTHAG WITH SINGLE CALL TO DNRM2!
-!* IF NOT CONVERGED, USE PERFORMANCE INDEX TO DECIDE ON ERROR
-!* VALUE SETTING, BUT DO NOT STOP!
-!* L.U. FOR ERROR MESSAGES PASSED THROUGH IERR
-!* USE PARAMETER STATEMENTS AND GENERIC INTRINSIC FUNCTIONS
-!* USE LEVEL 1 BLAS
-!* USE IF-THEN-ELSE TO CLARIFY LOGIC
-!* LOOP OVER SUBSPACES MADE INTO DO LOOP.
-!* LOOP OVER INVERSE ITERATIONS MADE INTO DO LOOP
-!* ZERO ONLY REQUIRED PORTIONS OF OUTPUT VECTOR
-!*
-!* NOTE -
-!* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
-!* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
-!*
-!
- use comm_par
- LOGICAL :: CONVGD !el,GOPARR,DSKWRK,MASWRK
-!
- INTEGER :: GROUP,I,IERR,ITS,J,JJ,M,N,NM,P,Q,R,S,SUBMAT,TAG
- INTEGER :: IND(M)
-!
- real(kind=8),dimension(N) :: D,E2
- real(kind=8) :: E(*)!el E(L)
- real(kind=8) :: W(M),Z(NM,M)
- real(kind=8),dimension(N) :: RV1,RV2,RV3,RV4,RV6
- real(kind=8) :: ANORM,EPS2,EPS3,EPS4,NORM,ORDER,RHO,U,UK,V
- real(kind=8) :: X0,X1,XU
-! real(kind=8) :: ESTPI1 !, DASUM, DDOT, DNRM2 EPSLON,
-!
-!el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
-!el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
-!
- real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00, GRPTOL = 0.001D+00
- real(kind=8),PARAMETER :: EPSCAL = 0.5D+00, HUNDRD = 100.0D+00, TEN = 10.0D+00
-!
- 001 FORMAT(' EIGENVECTOR ROUTINE EINVIT DID NOT CONVERGE FOR VECTOR' &
- ,I5,'. NORM =',1P,E10.2,' PERFORMANCE INDEX =',E10.2/ &
- ' (AN ERROR HALT WILL OCCUR IF THE PI IS GREATER THAN 100)')
- integer :: LUEMSG
-!
-!-----------------------------------------------------------------------
-!
- LUEMSG = IERR
- IERR = 0
- X0 = ZERO
- UK = ZERO
- NORM = ZERO
- EPS2 = ZERO
- EPS3 = ZERO
- EPS4 = ZERO
- GROUP = 0
- TAG = 0
- ORDER = ONE - E2(1)
- Q = 0
- DO 930 SUBMAT = 1, N
- P = Q + 1
-!
-! .......... ESTABLISH AND PROCESS NEXT SUBMATRIX ..........
-!
- DO 120 Q = P, N-1
- IF (E2(Q+1) .EQ. ZERO) GO TO 140
- 120 CONTINUE
- Q = N
-!
-! .......... FIND VECTORS BY INVERSE ITERATION ..........
-!
- 140 CONTINUE
- TAG = TAG + 1
- ANORM = ZERO
- S = 0
-!
- DO 920 R = 1, M
- IF (IND(R) .NE. TAG) GO TO 920
- ITS = 1
- X1 = W(R)
- IF (S .NE. 0) GO TO 510
-!
-! .......... CHECK FOR ISOLATED ROOT ..........
-!
- XU = ONE
- IF (P .EQ. Q) THEN
- RV6(P) = ONE
- CONVGD = .TRUE.
- GO TO 860
-!
- END IF
- NORM = ABS(D(P))
- DO 500 I = P+1, Q
- NORM = MAX( NORM, ABS(D(I)) + ABS(E(I)) )
- 500 CONTINUE
-!
-! .......... EPS2 IS THE CRITERION FOR GROUPING,
-! EPS3 REPLACES ZERO PIVOTS AND EQUAL
-! ROOTS ARE MODIFIED BY EPS3,
-! EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW .........
-!
- EPS2 = GRPTOL * NORM
- EPS3 = EPSCAL * EPSLON(NORM)
- UK = Q - P + 1
- EPS4 = UK * EPS3
- UK = EPS4 / SQRT(UK)
- S = P
- GROUP = 0
- GO TO 520
-!
-! .......... LOOK FOR CLOSE OR COINCIDENT ROOTS ..........
-!
- 510 IF (ABS(X1-X0) .GE. EPS2) THEN
-!
-! ROOTS ARE SEPERATE
-!
- GROUP = 0
- ELSE
-!
-! ROOTS ARE CLOSE
-!
- GROUP = GROUP + 1
- IF (ORDER * (X1 - X0) .LE. EPS3) X1 = X0 + ORDER * EPS3
- END IF
-!
-! .......... ELIMINATION WITH INTERCHANGES AND
-! INITIALIZATION OF VECTOR ..........
-!
- 520 CONTINUE
-!
- U = D(P) - X1
- V = E(P+1)
- RV6(P) = UK
- DO 550 I = P+1, Q
- RV6(I) = UK
- IF (ABS(E(I)) .GT. ABS(U)) THEN
-!
-! EXCHANGE ROWS BEFORE ELIMINATION
-!
-! *** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
-! E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY .......
-!
- XU = U / E(I)
- RV4(I) = XU
- RV1(I-1) = E(I)
- RV2(I-1) = D(I) - X1
- RV3(I-1) = E(I+1)
- U = V - XU * RV2(I-1)
- V = -XU * RV3(I-1)
-!
- ELSE
-!
-! STRAIGHT ELIMINATION
-!
- XU = E(I) / U
- RV4(I) = XU
- RV1(I-1) = U
- RV2(I-1) = V
- RV3(I-1) = ZERO
- U = D(I) - X1 - XU * V
- V = E(I+1)
- END IF
- 550 CONTINUE
-!
- IF (ABS(U) .LE. EPS3) U = EPS3
- RV1(Q) = U
- RV2(Q) = ZERO
- RV3(Q) = ZERO
-!
-! DO INVERSE ITERATIONS
-!
- CONVGD = .FALSE.
- DO 800 ITS = 1, 5
- IF (ITS .EQ. 1) GO TO 600
-!
-! .......... FORWARD SUBSTITUTION ..........
-!
- IF (NORM .EQ. ZERO) THEN
- RV6(S) = EPS4
- S = S + 1
- IF (S .GT. Q) S = P
- ELSE
- XU = EPS4 / NORM
- CALL DSCAL (Q-P+1, XU, RV6(P), 1)
- END IF
-!
-! ... ELIMINATION OPERATIONS ON NEXT VECTOR
-!
- DO 590 I = P+1, Q
- U = RV6(I)
-!
-! IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
-! WAS PERFORMED EARLIER IN THE
-! TRIANGULARIZATION PROCESS ..........
-!
- IF (RV1(I-1) .EQ. E(I)) THEN
- U = RV6(I-1)
- RV6(I-1) = RV6(I)
- ELSE
- U = RV6(I)
- END IF
- RV6(I) = U - RV4(I) * RV6(I-1)
- 590 CONTINUE
- 600 CONTINUE
-!
-! .......... BACK SUBSTITUTION
-!
- RV6(Q) = RV6(Q) / RV1(Q)
- V = U
- U = RV6(Q)
- NORM = ABS(U)
- DO 620 I = Q-1, P, -1
- RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
- V = U
- U = RV6(I)
- NORM = NORM + ABS(U)
- 620 CONTINUE
- IF (GROUP .EQ. 0) GO TO 700
-!
-! ....... ORTHOGONALIZE WITH RESPECT TO PREVIOUS
-! MEMBERS OF GROUP ..........
-!
- J = R
- DO 680 JJ = 1, GROUP
- 630 J = J - 1
- IF (IND(J) .NE. TAG) GO TO 630
- CALL DAXPY(Q-P+1, -DDOT(Q-P+1,RV6(P),1,Z(P,J),1),&
- Z(P,J),1,RV6(P),1)
- 680 CONTINUE
- NORM = DASUM(Q-P+1, RV6(P), 1)
- 700 CONTINUE
-!
- IF (CONVGD) GO TO 840
- IF (NORM .GE. ONE) CONVGD = .TRUE.
- 800 CONTINUE
-!
-! .......... NORMALIZE SO THAT SUM OF SQUARES IS
-! 1 AND EXPAND TO FULL ORDER ..........
-!
- 840 CONTINUE
-!
- XU = ONE / DNRM2(Q-P+1,RV6(P),1)
-!
- 860 CONTINUE
- DO 870 I = 1, P-1
- Z(I,R) = ZERO
- 870 CONTINUE
- DO 890 I = P,Q
- Z(I,R) = RV6(I) * XU
- 890 CONTINUE
- DO 900 I = Q+1, N
- Z(I,R) = ZERO
- 900 CONTINUE
-!
- IF (.NOT.CONVGD) THEN
- RHO = ESTPI1(Q-P+1,X1,D(P),E(P),Z(P,R),ANORM)
- IF (RHO .GE. TEN .AND. LUEMSG .GT. 0 .AND. MASWRK) &
- WRITE(LUEMSG,001) R,NORM,RHO
-!
-! *** SET ERROR -- NON-CONVERGED EIGENVECTOR ..........
-!
- IF (RHO .GT. HUNDRD) IERR = -R
- END IF
-!
- X0 = X1
- 920 CONTINUE
-!
- IF (Q .EQ. N) GO TO 940
- 930 CONTINUE
- 940 CONTINUE
- return
- end subroutine EINVIT
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK ELAUM
- subroutine ELAU(HINV,L,D,A,E)
-!
- integer :: L,JL,JK,J,JM1,K
- real(kind=8) :: A(*)
- real(kind=8) :: D(L)
-!el real(kind=8) :: E(L)
- real(kind=8) :: E(*)!el E(L)
- real(kind=8) :: F
- real(kind=8) :: G
- real(kind=8) :: HH
- real(kind=8) :: HINV
-!
- real(kind=8),PARAMETER :: ZERO = 0.0D+00, HALF = 0.5D+00
-!
- JL = L
- E(1) = A(1) * D(1)
- JK = 2
- DO 210 J = 2, JL
- F = D(J)
- G = ZERO
- JM1 = J - 1
-!
- DO 200 K = 1, JM1
- G = G + A(JK) * D(K)
- E(K) = E(K) + A(JK) * F
- JK = JK + 1
- 200 CONTINUE
-!
- E(J) = G + A(JK) * F
- JK = JK + 1
- 210 CONTINUE
-!
-! .......... FORM P ..........
-!
- F = ZERO
- DO 245 J = 1, L
- E(J) = E(J) * HINV
- F = F + E(J) * D(J)
- 245 CONTINUE
-!
-! .......... FORM Q ..........
-!
- HH = F * HALF * HINV
- DO 250 J = 1, L
- 250 E(J) = E(J) - HH * D(J)
-!
- return
- end subroutine ELAU
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK EPSLON
- real(kind=8) function EPSLON(X)
-!*
-!* AUTHORS -
-!* THIS ROUTINE WAS TAKEN FROM EISPACK EDITION 3 DATED 4/6/83
-!* THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE NOV 1986
-!*
-!* PURPOSE -
-!* ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X.
-!*
-!* ON ENTRY -
-!* X - WORKING PRECISION REAL
-!* VALUES TO FIND EPSLON FOR
-!*
-!* ON EXIT -
-!* EPSLON - WORKING PRECISION REAL
-!* SMALLEST POSITIVE VALUE SUCH THAT X+EPSLON .NE. ZERO
-!*
-!* QUALIFICATIONS -
-!* THIS ROUTINE SHOULD PERFORM PROPERLY ON ALL SYSTEMS
-!* SATISFYING THE FOLLOWING TWO ASSUMPTIONS,
-!* 1. THE BASE USED IN REPRESENTING FLOATING POINT
-!* NUMBERS IS NOT A POWER OF THREE.
-!* 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO
-!* THE ACCURACY USED IN FLOATING POINT VARIABLES
-!* THAT ARE STORED IN MEMORY.
-!* THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO
-!* FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING
-!* ASSUMPTION 2.
-!* UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT,
-!* A IS NOT EXACTLY EQUAL TO FOUR-THIRDS,
-!* B HAS A ZERO FOR ITS LAST BIT OR DIGIT,
-!* C IS NOT EXACTLY EQUAL TO ONE,
-!* EPS MEASURES THE SEPARATION OF 1.0 FROM
-!* THE NEXT LARGER FLOATING POINT NUMBER.
-!* THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED
-!* ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD.
-!*
-!* DIFFERENCES FROM EISPACK 3 -
-!* USE IS MADE OF PARAMETER STATEMENTS AND INTRINSIC FUNCTIONS
-!* --NO EXECUTEABLE CODE CHANGES--
-!*
-!* NOTE -
-!* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
-!* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
-!
- real(kind=8) :: A,B,C,EPS,X
-!
- real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00, THREE=3.0D+00, FOUR=4.0D+00
-!
-!-----------------------------------------------------------------------
-!
- A = FOUR/THREE
- 10 B = A - ONE
- C = B + B + B
- EPS = ABS(C - ONE)
- IF (EPS .EQ. ZERO) GO TO 10
- EPSLON = EPS*ABS(X)
- return
- end function EPSLON
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK EQLRAT
- subroutine EQLRAT(N,DIAG,E,E2IN,D,IND,IERR,E2)
-!*
-!* AUTHORS -
-!* THIS IS A MODIFICATION OF ROUTINE EQLRAT FROM EISPACK EDITION 3
-!* DATED AUGUST 1983.
-!* TQLRAT IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT,
-!* ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH.
-!* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
-!*
-!* PURPOSE -
-!* THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC
-!* TRIDIAGONAL MATRIX
-!*
-!* METHOD -
-!* RATIONAL QL
-!*
-!* ON ENTRY -
-!* N - INTEGER
-!* THE ORDER OF THE MATRIX.
-!* D - W.P. REAL (N)
-!* CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
-!* E2 - W.P. REAL (N)
-!* CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF
-!* THE INPUT MATRIX IN ITS LAST N-1 POSITIONS.
-!* E2(1) IS ARBITRARY.
-!*
-!* ON EXIT -
-!* D - W.P. REAL (N)
-!* CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
-!* ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
-!* ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
-!* THE SMALLEST EIGENVALUES.
-!* E2 - W.P. REAL (N)
-!* DESTROYED.
-!* IERR - INTEGER
-!* SET TO
-!* ZERO FOR NORMAL RETURN,
-!* J IF THE J-TH EIGENVALUE HAS NOT BEEN
-!* DETERMINED AFTER 30 ITERATIONS.
-!*
-!* DIFFERENCES FROM EISPACK 3 -
-!* G=G+B INSTEAD OF IF(G.EQ.0) G=B ; B=B/4
-!* F77 BACKWARD LOOPS INSTEAD OF F66 CONSTRUCT
-!* GENERIC INTRINSIC FUNCTIONS
-!* ARRARY IND ADDED FOR USE BY EINVIT
-!*
-!* NOTE -
-!* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
-!* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
-!
- INTEGER :: I,J,L,M,N,II,L1,IERR
- INTEGER,dimension(N) :: IND
-!
- real(kind=8),dimension(N) :: D,E,E2,DIAG,E2IN
- real(kind=8) :: B,C,F,G,H,P,R,S,T !,EPSLON
-!
- real(kind=8),PARAMETER :: ZERO = 0.0D+00, SCALE= 1.0D+00/64.0D+00, ONE = 1.0D+00
-!
- integer :: K,ITAG
-!-----------------------------------------------------------------------
- IERR = 0
- D(1)=DIAG(1)
- IND(1) = 1
- K = 0
- ITAG = 0
- IF (N .EQ. 1) GO TO 1001
-!
- DO 100 I = 2, N
- D(I)=DIAG(I)
- 100 E2(I-1) = E2IN(I)
-!
- F = ZERO
- T = ZERO
- B = EPSLON(ONE)
- C = B *B
- B = B * SCALE
- E2(N) = ZERO
-!
- DO 290 L = 1, N
- H = ABS(D(L)) + ABS(E(L))
- IF (T .GE. H) GO TO 105
- T = H
- B = EPSLON(T)
- C = B * B
- B = B * SCALE
- 105 CONTINUE
-! .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT ..........
- M = L - 1
- 110 M = M + 1
- IF (E2(M) .GT. C) GO TO 110
-! .......... E2(N) IS ALWAYS ZERO, SO THERE IS AN EXIT
-! FROM THE LOOP ..........
-!
- IF (M .LE. K) GO TO 125
- IF (M .NE. N) E2IN(M+1) = ZERO
- K = M
- ITAG = ITAG + 1
- 125 CONTINUE
- IF (M .EQ. L) GO TO 210
-!
-! ITERATE
-!
- DO 205 J = 1, 30
-! .......... FORM SHIFT ..........
- L1 = L + 1
- S = SQRT(E2(L))
- G = D(L)
- P = (D(L1) - G) / (2.0D+00 * S)
- R = SQRT(P*P+1.0D+00)
- D(L) = S / (P + SIGN(R,P))
- H = G - D(L)
-!
- DO 140 I = L1, N
- 140 D(I) = D(I) - H
-!
- F = F + H
-! .......... RATIONAL QL TRANSFORMATION ..........
- G = D(M) + B
- H = G
- S = ZERO
- DO 200 I = M-1,L,-1
- P = G * H
- R = P + E2(I)
- E2(I+1) = S * R
- S = E2(I) / R
- D(I+1) = H + S * (H + D(I))
- G = D(I) - E2(I) / G + B
- H = G * P / R
- 200 CONTINUE
-!
- E2(L) = S * G
- D(L) = H
-! .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST
- IF (H .EQ. ZERO) GO TO 210
- IF (ABS(E2(L)) .LE. ABS(C/H)) GO TO 210
- E2(L) = H * E2(L)
- IF (E2(L) .EQ. ZERO) GO TO 210
- 205 CONTINUE
-! .......... SET ERROR -- NO CONVERGENCE TO AN
-! EIGENVALUE AFTER 30 ITERATIONS ..........
- IERR = L
- GO TO 1001
-!
-! CONVERGED
-!
- 210 P = D(L) + F
-! .......... ORDER EIGENVALUES ..........
- I = 1
- IF (L .EQ. 1) GO TO 250
- IF (P .LT. D(1)) GO TO 230
- I = L
-! .......... LOOP TO FIND ORDERED POSITION
- 220 I = I - 1
- IF (P .LT. D(I)) GO TO 220
-!
- I = I + 1
- IF (I .EQ. L) GO TO 250
- 230 CONTINUE
- DO 240 II = L, I+1, -1
- D(II) = D(II-1)
- IND(II) = IND(II-1)
- 240 CONTINUE
-!
- 250 CONTINUE
- D(I) = P
- IND(I) = ITAG
- 290 CONTINUE
-!
- 1001 return
- end subroutine EQLRAT
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK ESTPI1
- real(kind=8) function ESTPI1(N,EVAL,D,E,X,ANORM)
-!*
-!* AUTHOR -
-!* STEPHEN T. ELBERT (AMES LABORATORY-USDOE) DATE: 5 DEC 1986
-!*
-!* PURPOSE -
-!* EVALUATE SYMMETRIC TRIDIAGONAL MATRIX PERFORMANCE INDEX
-!* * * * * *
-!* FOR 1 EIGENVECTOR
-!* *
-!*
-!* METHOD -
-!* THIS ROUTINE FORMS THE 1-NORM OF THE RESIDUAL MATRIX A*X-X*EVAL
-!* WHERE A IS A SYMMETRIC TRIDIAGONAL MATRIX STORED
-!* IN THE DIAGONAL (D) AND SUB-DIAGONAL (E) VECTORS, EVAL IS THE
-!* EIGENVALUE OF AN EIGENVECTOR OF A, NAMELY X.
-!* THIS NORM IS SCALED BY MACHINE ACCURACY FOR THE PROBLEM SIZE.
-!* ALL NORMS APPEARING IN THE COMMENTS BELOW ARE 1-NORMS.
-!*
-!* ON ENTRY -
-!* N - INTEGER
-!* THE ORDER OF THE MATRIX A.
-!* EVAL - W.P. REAL
-!* THE EIGENVALUE CORRESPONDING TO VECTOR X.
-!* D - W.P. REAL (N)
-!* THE DIAGONAL VECTOR OF A.
-!* E - W.P. REAL (N)
-!* THE SUB-DIAGONAL VECTOR OF A.
-!* X - W.P. REAL (N)
-!* AN EIGENVECTOR OF A.
-!* ANORM - W.P. REAL
-!* THE NORM OF A IF IT HAS BEEN PREVIOUSLY COMPUTED.
-!*
-!* ON EXIT -
-!* ANORM - W.P. REAL
-!* THE NORM OF A, COMPUTED IF INITIALLY ZERO.
-!* ESTPI1 - W.P. REAL
-!* !!A*X-X*EVAL!! / (EPSLON(10*N)*!!A!!*!!X!!);
-!* WHERE EPSLON(X) IS THE SMALLEST NUMBER SUCH THAT
-!* X + EPSLON(X) .NE. X
-!*
-!* ESTPI1 .LT. 1 == SATISFACTORY PERFORMANCE
-!* .GE. 1 AND .LE. 100 == MARGINAL PERFORMANCE
-!* .GT. 100 == POOR PERFORMANCE
-!* (SEE LECT. NOTES IN COMP. SCI. VOL.6 PP 124-125)
-!
- integer :: N,I
- real(kind=8) :: ANORM,EVAL,RNORM,SIZE,XNORM
- real(kind=8),dimension(N) :: D,X
- real(kind=8) :: E(*)!el E(L)
-! real(kind=8) :: EPSLON
-!
- real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00
-!
-!-----------------------------------------------------------------------
-!
- ESTPI1 = ZERO
- IF( N .LE. 1 ) RETURN
- SIZE = 10 * N
- IF (ANORM .EQ. ZERO) THEN
-!
-! COMPUTE NORM OF A
-!
- ANORM = MAX( ABS(D(1)) + ABS(E(2)), &
- ABS(D(N)) + ABS(E(N)))
- DO 110 I = 2, N-1
- ANORM = MAX( ANORM, ABS(E(I))+ABS(D(I))+ABS(E(I+1)))
- 110 CONTINUE
- IF(ANORM .EQ. ZERO) ANORM = ONE
- END IF
-!
-! COMPUTE NORMS OF RESIDUAL AND EIGENVECTOR
-!
- XNORM = ABS(X(1)) + ABS(X(N))
- RNORM = ABS( (D(1)-EVAL)*X(1) + E(2)*X(2)) &
- +ABS( (D(N)-EVAL)*X(N) + E(N)*X(N-1))
- DO 120 I = 2, N-1
- XNORM = XNORM + ABS(X(I))
- RNORM = RNORM + ABS(E(I)*X(I-1) + (D(I)-EVAL)*X(I) &
- + E(I+1)*X(I+1))
- 120 CONTINUE
-!
- ESTPI1 = RNORM / (EPSLON(SIZE)*ANORM*XNORM)
- return
- end function ESTPI1
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK ETRBK3
- subroutine ETRBK3(NM,N,NV,A,M,Z)
-!*
-!* AUTHORS-
-!* THIS IS A MODIFICATION OF ROUTINE TRBAK3 FROM EISPACK EDITION 3
-!* DATED AUGUST 1983.
-!* EISPACK TRBAK3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
-!* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
-!* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
-!* THIS VERSION IS BY S. T. ELBERT (AMES LABORATORY-USDOE)
-!*
-!* PURPOSE -
-!* THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
-!* MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
-!* SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY ETRED3.
-!*
-!* METHOD -
-!* THE CALCULATION IS CARRIED OUT BY FORMING THE MATRIX PRODUCT
-!* Q*Z
-!* WHERE Q IS A PRODUCT OF THE ORTHOGONAL SYMMETRIC MATRICES
-!* Q = PROD(I)[1 - U(I)*.TRANSPOSE.U(I)*H(I)]
-!* U IS THE AUGMENTED SUB-DIAGONAL ROWS OF A AND
-!* Z IS THE SET OF EIGENVECTORS OF THE TRIDIAGONAL
-!* MATRIX F WHICH WAS FORMED FROM THE ORIGINAL SYMMETRIC
-!* MATRIX C BY THE SIMILARITY TRANSFORMATION
-!* F = Q(TRANSPOSE) C Q
-!* NOTE THAT ETRBK3 PRESERVES VECTOR EUCLIDEAN NORMS.
-!*
-!*
-!* COMPLEXITY -
-!* M*N**2
-!*
-!* ON ENTRY-
-!* NM - INTEGER
-!* MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
-!* ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
-!* DIMENSION STATEMENT.
-!* N - INTEGER
-!* THE ORDER OF THE MATRIX A.
-!* NV - INTEGER
-!* MUST BE SET TO THE DIMENSION OF THE ARRAY A AS
-!* DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT.
-!* A - W.P. REAL (NV)
-!* CONTAINS INFORMATION ABOUT THE ORTHOGONAL
-!* TRANSFORMATIONS USED IN THE REDUCTION BY ETRED3 IN
-!* ITS FIRST NV = N*(N+1)/2 POSITIONS.
-!* M - INTEGER
-!* THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
-!* Z - W.P REAL (NM,M)
-!* CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
-!* IN ITS FIRST M COLUMNS.
-!*
-!* ON EXIT-
-!* Z - W.P. REAL (NM,M)
-!* CONTAINS THE TRANSFORMED EIGENVECTORS
-!* IN ITS FIRST M COLUMNS.
-!*
-!* DIFFERENCES WITH EISPACK 3 -
-!* THE TWO INNER LOOPS ARE REPLACED BY DDOT AND DAXPY.
-!* MULTIPLICATION USED INSTEAD OF DIVISION TO FIND S.
-!* OUTER LOOP RANGE CHANGED FROM 2,N TO 3,N.
-!* ADDRESS POINTERS FOR A SIMPLIFIED.
-!*
-!* NOTE -
-!* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
-!* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
-!
- INTEGER :: I,II,IM1,IZ,J,M,N,NM,NV
-!
- real(kind=8) :: A(NV),Z(NM,M)
- real(kind=8) :: H,S !,DDOT
-!
- real(kind=8),PARAMETER :: ZERO = 0.0D+00
-!
-!-----------------------------------------------------------------------
-!
- IF (M .EQ. 0) RETURN
- IF (N .LE. 2) RETURN
-!
- II=3
- DO 140 I = 3, N
- IZ=II+1
- II=II+I
- H = A(II)
- IF (H .EQ. ZERO) GO TO 140
- IM1 = I - 1
- DO 130 J = 1, M
- S = -( DDOT(IM1,A(IZ),1,Z(1,J),1) * H) * H
- CALL DAXPY(IM1,S,A(IZ),1,Z(1,J),1)
- 130 CONTINUE
- 140 CONTINUE
- return
- end subroutine ETRBK3
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK ETRED3
- subroutine ETRED3(N,NV,A,D,E,E2)
-!*
-!* AUTHORS -
-!* THIS IS A MODIFICATION OF ROUTINE TRED3 FROM EISPACK EDITION 3
-!* DATED AUGUST 1983.
-!* EISPACK TRED3 IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
-!* NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
-!* HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
-!* THIS VERSION IS BY S. T. ELBERT, AMES LABORATORY-USDOE JUN 1986
-!*
-!* PURPOSE -
-!* THIS ROUTINE REDUCES A REAL SYMMETRIC (PACKED) MATRIX, STORED
-!* AS A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
-!* USING ORTHOGONAL SIMILARITY TRANSFORMATIONS, PRESERVING THE
-!* INFORMATION ABOUT THE TRANSFORMATIONS IN A.
-!*
-!* METHOD -
-!* THE TRIDIAGONAL REDUCTION IS PERFORMED IN THE FOLLOWING WAY.
-!* STARTING WITH J=N, THE ELEMENTS IN THE J-TH ROW TO THE
-!* LEFT OF THE DIAGONAL ARE FIRST SCALED, TO AVOID POSSIBLE
-!* UNDERFLOW IN THE TRANSFORMATION THAT MIGHT RESULT IN SEVERE
-!* DEPARTURE FROM ORTHOGONALITY. THE SUM OF SQUARES SIGMA OF
-!* THESE SCALED ELEMENTS IS NEXT FORMED. THEN, A VECTOR U AND
-!* A SCALAR
-!* H = U(TRANSPOSE) * U / 2
-!* DEFINE A REFLECTION OPERATOR
-!* P = I - U * U(TRANSPOSE) / H
-!* WHICH IS ORTHOGONAL AND SYMMETRIC AND FOR WHICH THE
-!* SIMILIARITY TRANSFORMATION PAP ELIMINATES THE ELEMENTS IN
-!* THE J-TH ROW OF A TO THE LEFT OF THE SUBDIAGONAL AND THE
-!* SYMMETRICAL ELEMENTS IN THE J-TH COLUMN.
-!*
-!* THE NON-ZERO COMPONENTS OF U ARE THE ELEMENTS OF THE J-TH
-!* ROW TO THE LEFT OF THE DIAGONAL WITH THE LAST OF THEM
-!* AUGMENTED BY THE SQUARE ROOT OF SIGMA PREFIXED BY THE SIGN
-!* OF THE SUBDIAGONAL ELEMENT. BY STORING THE TRANSFORMED SUB-
-!* DIAGONAL ELEMENT IN E(J) AND NOT OVERWRITING THE ROW
-!* ELEMENTS ELIMINATED IN THE TRANSFORMATION, FULL INFORMATION
-!* ABOUT P IS SAVE FOR LATER USE IN ETRBK3.
-!*
-!* THE TRANSFORMATION SETS E2(J) EQUAL TO SIGMA AND E(J)
-!* EQUAL TO THE SQUARE ROOT OF SIGMA PREFIXED BY THE SIGN
-!* OF THE REPLACED SUBDIAGONAL ELEMENT.
-!*
-!* THE ABOVE STEPS ARE REPEATED ON FURTHER ROWS OF THE
-!* TRANSFORMED A IN REVERSE ORDER UNTIL A IS REDUCED TO TRI-
-!* DIAGONAL FORM, THAT IS, REPEATED FOR J = N-1,N-2,...,3.
-!*
-!* COMPLEXITY -
-!* 2/3 N**3
-!*
-!* ON ENTRY-
-!* N - INTEGER
-!* THE ORDER OF THE MATRIX.
-!* NV - INTEGER
-!* MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
-!* AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT
-!* A - W.P. REAL (NV)
-!* CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
-!* INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
-!* ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
-!*
-!* ON EXIT-
-!* A - W.P. REAL (NV)
-!* CONTAINS INFORMATION ABOUT THE ORTHOGONAL
-!* TRANSFORMATIONS USED IN THE REDUCTION.
-!* D - W.P. REAL (N)
-!* CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL
-!* MATRIX.
-!* E - W.P. REAL (N)
-!* CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
-!* MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO
-!* E2 - W.P. REAL (N)
-!* CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF
-!* E. MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
-!*
-!* DIFFERENCES FROM EISPACK 3 -
-!* OUTER LOOP CHANGED FROM II=1,N TO I=N,3,-1
-!* PARAMETER STATEMENT AND GENERIC INTRINSIC FUNCTIONS USED
-!* SCALE.NE.0 TEST NOW SPOTS TRI-DIAGONAL FORM
-!* VALUES LESS THAN EPSLON CLEARED TO ZERO
-!* USE BLAS(1)
-!* U NOT COPIED TO D, LEFT IN A
-!* E2 COMPUTED FROM E
-!* INNER LOOPS SPLIT INTO ROUTINES ELAU AND FREDA
-!* INVERSE OF H STORED INSTEAD OF H
-!*
-!* NOTE -
-!* QUESTIONS AND COMMENTS CONCERNING EISPACK SHOULD BE DIRECTED TO
-!* B. S. GARBOW, APPLIED MATH. DIVISION, ARGONNE NATIONAL LAB.
-!
- INTEGER :: I,IIA,IZ0,L,N,NV
-!
- real(kind=8) :: A(NV),D(N),E2(N)
- real(kind=8) :: E(*)!el E(L)
- real(kind=8) :: AIIMAX,F,G,H,HROOT,SCALE,SCALEI
-! real(kind=8) :: DASUM, DNRM2
-!
- real(kind=8),PARAMETER :: ZERO = 0.0D+00, ONE = 1.0D+00
-!
-!-----------------------------------------------------------------------
-!
- IF (N .LE. 2) GO TO 310
- IZ0 = (N*N+N)/2
- AIIMAX = ABS(A(IZ0))
- DO 300 I = N, 3, -1
- L = I - 1
- IIA = IZ0
- IZ0 = IZ0 - I
- AIIMAX = MAX(AIIMAX, ABS(A(IIA)))
- SCALE = DASUM (L, A(IZ0+1), 1)
- IF(SCALE .EQ. ABS(A(IIA-1)) .OR. AIIMAX+SCALE .EQ. AIIMAX) THEN
-!
-! THIS ROW IS ALREADY IN TRI-DIAGONAL FORM
-!
- D(I) = A(IIA)
- IF (AIIMAX+D(I) .EQ. AIIMAX) D(I) = ZERO
- E(I) = A(IIA-1)
- IF (AIIMAX+E(I) .EQ. AIIMAX) E(I) = ZERO
- E2(I) = E(I)*E(I)
- A(IIA) = ZERO
- GO TO 300
-!
- END IF
-!
- SCALEI = ONE / SCALE
- CALL DSCAL(L,SCALEI,A(IZ0+1),1)
- HROOT = DNRM2(L,A(IZ0+1),1)
-!
- F = A(IZ0+L)
- G = -SIGN(HROOT,F)
- E(I) = SCALE * G
- E2(I) = E(I)*E(I)
- H = HROOT*HROOT - F * G
- A(IZ0+L) = F - G
- D(I) = A(IIA)
- A(IIA) = ONE / SQRT(H)
-! .......... FORM P THEN Q IN E(1:L) ..........
- CALL ELAU(ONE/H,L,A(IZ0+1),A,E)
-! .......... FORM REDUCED A ..........
- CALL FREDA(L,A(IZ0+1),A,E)
-!
- 300 CONTINUE
- 310 CONTINUE
- E(1) = ZERO
- E2(1)= ZERO
- D(1) = A(1)
- IF(N.EQ.1) RETURN
-!
- E(2) = A(2)
- E2(2)= A(2)*A(2)
- D(2) = A(3)
- return
- end subroutine ETRED3
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK EVVRSP
- subroutine EVVRSP(MSGFL,N,NVECT,LENA,NV,A,B,IND,ROOT,VECT,IORDER,IERR)
-!*
-!* AUTHOR: S. T. ELBERT, AMES LABORATORY-USDOE, JUNE 1985
-!*
-!* PURPOSE -
-!* FINDS (ALL) EIGENVALUES AND (SOME OR ALL) EIGENVECTORS
-!* * * *
-!* OF A REAL SYMMETRIC PACKED MATRIX.
-!* * * *
-!*
-!* METHOD -
-!* THE METHOD AS PRESENTED IN THIS ROUTINE CONSISTS OF FOUR STEPS:
-!* FIRST, THE INPUT MATRIX IS REDUCED TO TRIDIAGONAL FORM BY THE
-!* HOUSEHOLDER TECHNIQUE (ORTHOGONAL SIMILARITY TRANSFORMATIONS).
-!* SECOND, THE ROOTS ARE LOCATED USING THE RATIONAL QL METHOD.
-!* THIRD, THE VECTORS OF THE TRIDIAGONAL FORM ARE EVALUATED BY THE
-!* INVERSE ITERATION TECHNIQUE. VECTORS FOR DEGENERATE OR NEAR-
-!* DEGENERATE ROOTS ARE FORCED TO BE ORTHOGONAL.
-!* FOURTH, THE TRIDIAGONAL VECTORS ARE ROTATED TO VECTORS OF THE
-!* ORIGINAL ARRAY.
-!*
-!* THESE ROUTINES ARE MODIFICATIONS OF THE EISPACK 3
-!* ROUTINES TRED3, TQLRAT, TINVIT AND TRBAK3
-!*
-!* FOR FURTHER DETAILS, SEE EISPACK USERS GUIDE, B. T. SMITH
-!* ET AL, SPRINGER-VERLAG, LECTURE NOTES IN COMPUTER SCIENCE,
-!* VOL. 6, 2-ND EDITION, 1976. ANOTHER GOOD REFERENCE IS
-!* THE SYMMETRIC EIGENVALUE PROBLEM BY B. N. PARLETT
-!* PUBLISHED BY PRENTICE-HALL, INC., ENGLEWOOD CLIFFS, N.J. (1980)
-!*
-!* ON ENTRY -
-!* MSGFL - INTEGER (LOGICAL UNIT NO.)
-!* FILE WHERE ERROR MESSAGES WILL BE PRINTED.
-!* IF MSGFL IS 0, ERROR MESSAGES WILL BE PRINTED ON LU 6.
-!* IF MSGFL IS NEGATIVE, NO ERROR MESSAGES PRINTED.
-!* N - INTEGER
-!* ORDER OF MATRIX A.
-!* NVECT - INTEGER
-!* NUMBER OF VECTORS DESIRED. 0 .LE. NVECT .LE. N.
-!* LENA - INTEGER
-!* DIMENSION OF A IN CALLING ROUTINE. MUST NOT BE LESS
-!* THAN (N*N+N)/2.
-!* NV - INTEGER
-!* ROW DIMENSION OF VECT IN CALLING ROUTINE. N .LE. NV.
-!* A - WORKING PRECISION REAL (LENA)
-!* INPUT MATRIX, ROWS OF THE LOWER TRIANGLE PACKED INTO
-!* LINEAR ARRAY OF DIMENSION N*(N+1)/2. THE PACKED ORDER
-!* IS A(1,1), A(2,1), A(2,2), A(3,1), A(3,2), ...
-!* B - WORKING PRECISION REAL (N,8)
-!* SCRATCH ARRAY, 8*N ELEMENTS
-!* IND - INTEGER (N)
-!* SCRATCH ARRAY OF LENGTH N.
-!* IORDER - INTEGER
-!* ROOT ORDERING FLAG.
-!* = 0, ROOTS WILL BE PUT IN ASCENDING ORDER.
-!* = 2, ROOTS WILL BE PUT IN DESCENDING ORDER.
-!*
-!* ON EXIT -
-!* A - DESTORYED. NOW HOLDS REFLECTION OPERATORS.
-!* ROOT - WORKING PRECISION REAL (N)
-!* ALL EIGENVALUES IN ASCENDING OR DESCENDING ORDER.
-!* IF IORDER = 0, ROOT(1) .LE. ... .LE. ROOT(N)
-!* IF IORDER = 2, ROOT(1) .GE. ... .GE. ROOT(N)
-!* VECT - WORKING PRECISION REAL (NV,NVECT)
-!* EIGENVECTORS FOR ROOT(1), ..., ROOT(NVECT).
-!* IERR - INTEGER
-!* = 0 IF NO ERROR DETECTED,
-!* = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
-!* = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
-!* (FAILURES SHOULD BE VERY RARE. CONTACT C. MOLER.)
-!*
-!
- use comm_par
-!el LOGICAL :: GOPARR,DSKWRK,MASWRK
-!
- integer :: MSGFL,N,NVECT,LENA,NV,IORDER,IERR
- real(kind=8) :: A(LENA)
- real(kind=8) :: B(N,8)
- real(kind=8) :: ROOT(N)
- real(kind=8) :: T
- real(kind=8) :: VECT(NV,*)
-!
- INTEGER :: IND(N)
-!
-!el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
- real(kind=8) :: DSKW
-
-!el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
-!
- 900 FORMAT(26H0*** EVVRSP PARAMETERS ***/ &
- 14H *** N = ,I8,4H ***/ &
- 14H *** NVECT = ,I8,4H ***/ &
- 14H *** LENA = ,I8,4H ***/ &
- 14H *** NV = ,I8,4H ***/ &
- 14H *** IORDER = ,I8,4H ***/ &
- 14H *** IERR = ,I8,4H ***)
- 901 FORMAT(37H VALUE OF LENA IS LESS THAN (N*N+N)/2)
- 902 FORMAT(39H EQLRAT HAS FAILED TO CONVERGE FOR ROOT,I5)
- 903 FORMAT(18H NV IS LESS THAN N)
- 904 FORMAT(41H EINVIT HAS FAILED TO CONVERGE FOR VECTOR,I5)
- 905 FORMAT(51H VALUE OF IORDER MUST BE 0 (SMALLEST ROOT FIRST) OR &
- ,23H 2 (LARGEST ROOT FIRST))
- 906 FORMAT(' VALUE OF N IS LESS THAN OR EQUAL ZERO')
-
- integer :: LMSGFL,I,J,L,JSV,KLIM,K
-!
-!-----------------------------------------------------------------------
-!
- LMSGFL=MSGFL
- IF (MSGFL .EQ. 0) LMSGFL=6
- IERR = N - 1
- IF (N .LE. 0) GO TO 800
- IERR = N + 1
- IF ( (N*N+N)/2 .GT. LENA) GO TO 810
-!
-! REDUCE REAL SYMMETRIC MATRIX A TO TRIDIAGONAL FORM
-!
- CALL ETRED3(N,LENA,A,B(1,1),B(1,2),B(1,3))
-!
-! FIND ALL EIGENVALUES OF TRIDIAGONAL MATRIX
-!
- CALL EQLRAT(N,B(1,1),B(1,2),B(1,3),ROOT,IND,IERR,B(1,4))
- IF (IERR .NE. 0) GO TO 820
-!
-! CHECK THE DESIRED ORDER OF THE EIGENVALUES
-!
- B(1,3) = IORDER
- IF (IORDER .EQ. 0) GO TO 300
- IF (IORDER .NE. 2) GO TO 850
-!
-! ORDER ROOTS IN DESCENDING ORDER (LARGEST FIRST)...
-! TURN ROOT AND IND ARRAYS END FOR END
-!
- DO 210 I = 1, N/2
- J = N+1-I
- T = ROOT(I)
- ROOT(I) = ROOT(J)
- ROOT(J) = T
- L = IND(I)
- IND(I) = IND(J)
- IND(J) = L
- 210 CONTINUE
-!
-! FIND I AND J MARKING THE START AND END OF A SEQUENCE
-! OF DEGENERATE ROOTS
-!
- I=0
- 220 CONTINUE
- I = I+1
- IF (I .GT. N) GO TO 300
- DO 230 J=I,N
- IF (ROOT(J) .NE. ROOT(I)) GO TO 240
- 230 CONTINUE
- J = N+1
- 240 CONTINUE
- J = J-1
- IF (J .EQ. I) GO TO 220
-!
-! TURN AROUND IND BETWEEN I AND J
-!
- JSV = J
- KLIM = (J-I+1)/2
- DO 250 K=1,KLIM
- L = IND(J)
- IND(J) = IND(I)
- IND(I) = L
- I = I+1
- J = J-1
- 250 CONTINUE
- I = JSV
- GO TO 220
-!
- 300 CONTINUE
-!
- IF (NVECT .LE. 0) RETURN
- IF (NV .LT. N) GO TO 830
-!
-! FIND EIGENVECTORS OF TRI-DIAGONAL MATRIX VIA INVERSE ITERATION
-!
- IERR = LMSGFL
- CALL EINVIT(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,IND,&
- VECT,IERR,B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
- IF (IERR .NE. 0) GO TO 840
-!
-! FIND EIGENVECTORS OF SYMMETRIC MATRIX VIA BACK TRANSFORMATION
-!
- 400 CONTINUE
- CALL ETRBK3(NV,N,LENA,A,NVECT,VECT)
- RETURN
-!
-! ERROR MESSAGE SECTION
-!
- 800 IF (LMSGFL .LT. 0) RETURN
- IF (MASWRK) WRITE(LMSGFL,906)
- GO TO 890
-!
- 810 IF (LMSGFL .LT. 0) RETURN
- IF (MASWRK) WRITE(LMSGFL,901)
- GO TO 890
-!
- 820 IF (LMSGFL .LT. 0) RETURN
- IF (MASWRK) WRITE(LMSGFL,902) IERR
- GO TO 890
-!
- 830 IF (LMSGFL .LT. 0) RETURN
- IF (MASWRK) WRITE(LMSGFL,903)
- GO TO 890
-!
- 840 CONTINUE
- IF ((LMSGFL .GT. 0).AND.MASWRK) WRITE(LMSGFL,904) -IERR
- GO TO 400
-!
- 850 IERR=-1
- IF (LMSGFL .LT. 0) RETURN
- IF (MASWRK) WRITE(LMSGFL,905)
- GO TO 890
-!
- 890 CONTINUE
- IF (MASWRK) WRITE(LMSGFL,900) N,NVECT,LENA,NV,IORDER,IERR
- return
- end subroutine EVVRSP
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK FREDA
- subroutine FREDA(L,D,A,E)
-!
- integer :: l,jk,j,k
- real(kind=8) :: A(*)
- real(kind=8) :: D(L)
- real(kind=8) :: E(*)!el E(L)
- real(kind=8) :: F
- real(kind=8) :: G
-!
- JK = 1
-!
-! .......... FORM REDUCED A ..........
-!
- DO 280 J = 1, L
- F = D(J)
- G = E(J)
-!
- DO 260 K = 1, J
- A(JK) = A(JK) - F * E(K) - G * D(K)
- JK = JK + 1
- 260 CONTINUE
-!
- 280 CONTINUE
- return
- end subroutine FREDA
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK GIVEIS
- subroutine GIVEIS(N,NVECT,NV,A,B,INDB,ROOT,VECT,IERR)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,NVECT,NV,IERR
- real(kind=8),DIMENSION(*) :: A
- real(kind=8),DIMENSION(N,8) :: B
- integer,DIMENSION(N) :: INDB
- real(kind=8),DIMENSION(N) :: ROOT
- real(kind=8),DIMENSION(NV,NVECT) :: VECT
-!
-! EISPACK-BASED SUBSTITUTE FOR QCPE ROUTINE GIVENS.
-! FINDS ALL EIGENVALUES AND SOME EIGENVECTORS OF A REAL SYMMETRIC
-! MATRIX. AUTHOR.. C. MOLER AND D. SPANGLER, N.R.C.C., 4/1/79.
-!
-! INPUT..
-! N = ORDER OF MATRIX .
-! NVECT = NUMBER OF VECTORS DESIRED. 0 .LE. NVECT .LE. N .
-! NV = LEADING DIMENSION OF VECT .
-! A = INPUT MATRIX, COLUMNS OF THE UPPER TRIANGLE PACKED INTO
-! LINEAR ARRAY OF DIMENSION N*(N+1)/2 .
-! B = SCRATCH ARRAY, 8*N ELEMENTS (NOTE THIS IS MORE THAN
-! PREVIOUS VERSIONS OF GIVENS.)
-! IND = INDEX ARRAY OF N ELEMENTS
-!
-! OUTPUT..
-! A DESTROYED .
-! ROOT = ALL EIGENVALUES, ROOT(1) .LE. ... .LE. ROOT(N) .
-! (FOR OTHER ORDERINGS, SEE BELOW.)
-! VECT = EIGENVECTORS FOR ROOT(1),..., ROOT(NVECT) .
-! IERR = 0 IF NO ERROR DETECTED,
-! = K IF ITERATION FOR K-TH EIGENVALUE FAILED,
-! = -K IF ITERATION FOR K-TH EIGENVECTOR FAILED.
-! (FAILURES SHOULD BE VERY RARE. CONTACT MOLER.)
-!
-! CALLS MODIFIED EISPACK ROUTINES TRED3B, IMTQLV, TINVTB, AND
-! TRBK3B. THE ROUTINES TRED3B, TINVTB, AND TRBK3B.
-! THE ORIGINAL EISPACK ROUTINES TRED3, TINVIT, AND TRBAK3
-! WERE MODIFIED BY THE INTRODUCTION OF TWO ROUTINES FROM THE
-! BLAS LIBRARY - DDOT AND DAXPY.
-!
-! IF TINVIT FAILS TO CONVERGE, TQL2 IS CALLED
-!
-! SEE EISPACK USERS GUIDE, B. T. SMITH ET AL, SPRINGER-VERLAG
-! LECTURE NOTES IN COMPUTER SCIENCE, VOL. 6, 2-ND EDITION, 1976 .
-! NOTE THAT IMTQLV AND TINVTB HAVE INTERNAL MACHINE
-! DEPENDENT CONSTANTS.
-!
-!el DATA ONE, ZERO /1.0D+00, 0.0D+00/
- real(kind=8) :: ZERO = 0.0D+00, ONE = 1.0D+00
-
- integer :: i,j
-
- CALL TRED3B(N,(N*N+N)/2,A,B(1,1),B(1,2),B(1,3))
- CALL IMTQLV(N,B(1,1),B(1,2),B(1,3),ROOT,INDB,IERR,B(1,4))
- IF (IERR .NE. 0) RETURN
-!
-! TO REORDER ROOTS...
-! K = N/2
-! B(1,3) = 2.0D+00
-! DO 50 I = 1, K
-! J = N+1-I
-! T = ROOT(I)
-! ROOT(I) = ROOT(J)
-! ROOT(J) = T
-! 50 CONTINUE
-!
- IF (NVECT .LE. 0) RETURN
- CALL TINVTB(NV,N,B(1,1),B(1,2),B(1,3),NVECT,ROOT,INDB,VECT,IERR,&
- B(1,4),B(1,5),B(1,6),B(1,7),B(1,8))
- IF (IERR .EQ. 0) GO TO 160
-!
-! IF INVERSE ITERATION GIVES AN ERROR IN DETERMINING THE
-! EIGENVECTORS, TRY THE QL ALGORITHM IF ALL THE EIGENVECTORS
-! ARE DESIRED.
-!
- IF (NVECT .NE. N) RETURN
- DO 120 I = 1, NVECT
- DO 100 J = 1, N
- VECT(I,J) = ZERO
- 100 CONTINUE
- VECT(I,I) = ONE
- 120 CONTINUE
- CALL TQL2 (NV,N,B(1,1),B(1,2),VECT,IERR)
- DO 140 I = 1, NVECT
- ROOT(I) = B(I,1)
- 140 CONTINUE
- IF (IERR .NE. 0) RETURN
- 160 CALL TRBK3B(NV,N,(N*N+N)/2,A,NVECT,VECT)
- return
- end subroutine GIVEIS
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK GLDIAG
- subroutine GLDIAG(LDVECT,NVECT,N,H,WRK,EIG,VECTOR,IERR,IWRK)
-!
-! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-!
- use comm_iofile
- use comm_machsw
- use comm_par
-!el LOGICAL :: GOPARR,DSKWRK,MASWRK
-!
- integer :: LDVECT,NVECT,N,IERR
- real(kind=8),DIMENSION(*) :: H
- real(kind=8),DIMENSION(N,8) :: WRK
- real(kind=8),DIMENSION(N) :: EIG
- integer,DIMENSION(N) :: IWRK
- real(kind=8),DIMENSION(LDVECT,NVECT) :: VECTOR
-!
-!el integer :: IR,IW,IP,IJK,IPK,IDAF,NAV,IODA(400)
-!el integer :: KDIAG,ICORFL,IXDR
-!el COMMON /IOFILE/ IR,IW,IP,IJK,IPK,IDAF,NAV,IODA
-!el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
-!el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
-!el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
-!
- integer :: LENH,KORDER
-
-! ----- GENERAL ROUTINE TO DIAGONALIZE A SYMMETRIC MATRIX -----
-! IF KDIAG = 0, USE A ROUTINE FROM THE VECTOR LIBRARY,
-! IF AVAILABLE (SEE THE SUBROUTINE 'GLDIAG'
-! IN VECTOR.SRC), OR EVVRSP OTHERWISE
-! = 1, USE EVVRSP
-! = 2, USE GIVEIS
-! = 3, USE JACOBI
-!
-! N = DIMENSION (ORDER) OF MATRIX TO BE SOLVED
-! LDVECT = LEADING DIMENSION OF VECTOR
-! NVECT = NUMBER OF VECTORS DESIRED
-! H = MATRIX TO BE DIAGONALIZED
-! WRK = N*8 W.P. REAL WORDS OF SCRATCH SPACE
-! EIG = EIGENVECTORS (OUTPUT)
-! VECTOR = EIGENVECTORS (OUTPUT)
-! IERR = ERROR FLAG (OUTPUT)
-! IWRK = N INTEGER WORDS OF SCRATCH SPACE
-!
- IERR = 0
-!
-! ----- USE STEVE ELBERT'S ROUTINE -----
-!
- IF(KDIAG.LE.1 .OR. KDIAG.GT.3) THEN
- LENH = (N*N+N)/2
- KORDER =0
- CALL EVVRSP(IW,N,NVECT,LENH,LDVECT,H,WRK,IWRK,EIG,VECTOR, &
- KORDER,IERR)
- END IF
-!
-! ----- USE MODIFIED EISPAK ROUTINE -----
-!
- IF(KDIAG.EQ.2) &
- CALL GIVEIS(N,NVECT,LDVECT,H,WRK,IWRK,EIG,VECTOR,IERR)
-!
-! ----- USE JACOBI ROTATION ROUTINE -----
-!
- IF(KDIAG.EQ.3) THEN
- IF(NVECT.EQ.N) THEN
- CALL JACDG(H,VECTOR,EIG,IWRK,WRK,LDVECT,N)
- ELSE
- IF (MASWRK) WRITE(IW,9000) N,NVECT,LDVECT
- CALL ABRT
- END IF
- END IF
- RETURN
-!
- 9000 FORMAT(1X,'IN -GLDIAG-, N,NVECT,LDVECT=',3I8/ &
- 1X,'THE JACOBI CODE CANNOT COPE WITH N.NE.NVECT!'/ &
- 1X,'SO THIS RUN DOES NOT PERMIT KDIAG=3.')
- end subroutine GLDIAG
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK IMTQLV
- subroutine IMTQLV(N,D,E,E2,W,IND,IERR,RV1)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- INTEGER :: N,TAG,IERR
- real(kind=8) :: MACHEP
- real(kind=8),DIMENSION(N) :: D,E2,W,RV1
- real(kind=8) :: E(*)!el E(L)
- integer,DIMENSION(N) :: IND
- integer :: k,i,l,j,m,mml,ii
- real(kind=8) :: c,p,s,f,b,g,r
-!
-! THIS ROUTINE IS A VARIANT OF IMTQL1 WHICH IS A TRANSLATION OF
-! ALGOL PROCEDURE IMTQL1, NUM. MATH. 12, 377-383(1968) BY MARTIN AND
-! WILKINSON, AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE.
-! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971).
-!
-! THIS ROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC TRIDIAGONAL
-! MATRIX BY THE IMPLICIT QL METHOD AND ASSOCIATES WITH THEM
-! THEIR CORRESPONDING SUBMATRIX INDICES.
-!
-! ON INPUT-
-!
-! N IS THE ORDER OF THE MATRIX,
-!
-! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
-!
-! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
-! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
-!
-! E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
-! E2(1) IS ARBITRARY.
-!
-! ON OUTPUT-
-!
-! D AND E ARE UNALTERED,
-!
-! ELEMENTS OF E2, CORRESPONDING TO ELEMENTS OF E REGARDED
-! AS NEGLIGIBLE, HAVE BEEN REPLACED BY ZERO CAUSING THE
-! MATRIX TO SPLIT INTO A DIRECT SUM OF SUBMATRICES.
-! E2(1) IS ALSO SET TO ZERO,
-!
-! W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
-! ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND
-! ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE
-! THE SMALLEST EIGENVALUES,
-!
-! IND CONTAINS THE SUBMATRIX INDICES ASSOCIATED WITH THE
-! CORRESPONDING EIGENVALUES IN W -- 1 FOR EIGENVALUES
-! BELONGING TO THE FIRST SUBMATRIX FROM THE TOP,
-! 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.,
-!
-! IERR IS SET TO
-! ZERO FOR NORMAL RETURN,
-! J IF THE J-TH EIGENVALUE HAS NOT BEEN
-! DETERMINED AFTER 30 ITERATIONS,
-!
-! RV1 IS A TEMPORARY STORAGE ARRAY.
-!
-! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
-! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
-!
-! ------------------------------------------------------------------
-!
-! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
-! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
-!
-! **********
- MACHEP = 2.0D+00**(-50)
-!
- IERR = 0
- K = 0
- TAG = 0
-!
- DO 100 I = 1, N
- W(I) = D(I)
- IF (I .NE. 1) RV1(I-1) = E(I)
- 100 CONTINUE
-!
- E2(1) = 0.0D+00
- RV1(N) = 0.0D+00
-!
- DO 360 L = 1, N
- J = 0
-! ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
- 120 DO 140 M = L, N
- IF (M .EQ. N) GO TO 160
- IF (ABS(RV1(M)) .LE. MACHEP * (ABS(W(M)) + ABS(W(M+1)))) GO TO &
- 160
-! ********** GUARD AGAINST UNDERFLOWED ELEMENT OF E2 **********
- IF (E2(M+1) .EQ. 0.0D+00) GO TO 180
- 140 CONTINUE
-!
- 160 IF (M .LE. K) GO TO 200
- IF (M .NE. N) E2(M+1) = 0.0D+00
- 180 K = M
- TAG = TAG + 1
- 200 P = W(L)
- IF (M .EQ. L) GO TO 280
- IF (J .EQ. 30) GO TO 380
- J = J + 1
-! ********** FORM SHIFT **********
- G = (W(L+1) - P) / (2.0D+00 * RV1(L))
- R = SQRT(G*G+1.0D+00)
- G = W(M) - P + RV1(L) / (G + SIGN(R,G))
- S = 1.0D+00
- C = 1.0D+00
- P = 0.0D+00
- MML = M - L
-! ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
- DO 260 II = 1, MML
- I = M - II
- F = S * RV1(I)
- B = C * RV1(I)
- IF (ABS(F) .LT. ABS(G)) GO TO 220
- C = G / F
- R = SQRT(C*C+1.0D+00)
- RV1(I+1) = F * R
- S = 1.0D+00 / R
- C = C * S
- GO TO 240
- 220 S = F / G
- R = SQRT(S*S+1.0D+00)
- RV1(I+1) = G * R
- C = 1.0D+00 / R
- S = S * C
- 240 G = W(I+1) - P
- R = (W(I) - G) * S + 2.0D+00 * C * B
- P = S * R
- W(I+1) = G + P
- G = C * R - B
- 260 CONTINUE
-!
- W(L) = W(L) - P
- RV1(L) = G
- RV1(M) = 0.0D+00
- GO TO 120
-! ********** ORDER EIGENVALUES **********
- 280 IF (L .EQ. 1) GO TO 320
-! ********** FOR I=L STEP -1 UNTIL 2 DO -- **********
- DO 300 II = 2, L
- I = L + 2 - II
- IF (P .GE. W(I-1)) GO TO 340
- W(I) = W(I-1)
- IND(I) = IND(I-1)
- 300 CONTINUE
-!
- 320 I = 1
- 340 W(I) = P
- IND(I) = TAG
- 360 CONTINUE
-!
- GO TO 400
-! ********** SET ERROR -- NO CONVERGENCE TO AN
-! EIGENVALUE AFTER 30 ITERATIONS **********
- 380 IERR = L
- 400 return
-! ********** LAST CARD OF IMTQLV **********
- end subroutine IMTQLV
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK JACDG
- subroutine JACDG(A,VEC,EIG,JBIG,BIG,LDVEC,N)
-!
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
-!
- integer :: LDVEC,N
- real(kind=8),DIMENSION(*) :: A
- real(kind=8),DIMENSION(LDVEC,N) :: VEC
- real(kind=8),DIMENSION(N) :: EIG,BIG
- integer,DIMENSION(N) :: JBIG
-!
- real(kind=8),PARAMETER :: ONE=1.0D+00
- integer :: i,NB1,NB2,NMIN,NMAX
-!
-! ----- JACOBI DIAGONALIZATION OF SYMMETRIC MATRIX -----
-! SYMMETRIC MATRIX -A- OF DIMENSION -N- IS DESTROYED ON EXIT.
-! ALL EIGENVECTORS ARE FOUND, SO -VEC- MUST BE SQUARE,
-! UNLESS SOMEONE TAKES THE TROUBLE TO LOOK AT -NMAX- BELOW.
-! -BIG- AND -JBIG- ARE SCRATCH WORK ARRAYS.
-!
- CALL VCLR(VEC,1,LDVEC*N)
- DO 20 I = 1,N
- VEC(I,I) = ONE
- 20 CONTINUE
-!
- NB1 = N
- NB2 = (NB1*NB1+NB1)/2
- NMIN = 1
- NMAX = NB1
-!
- CALL JACDIA(A,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
-!
- DO 30 I=1,N
- EIG(I) = A((I*I+I)/2)
- 30 CONTINUE
-!
- CALL JACORD(VEC,EIG,NB1,LDVEC)
- return
- end subroutine JACDG
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK JACDIA
- subroutine JACDIA(F,VEC,NB1,NB2,LDVEC,NMIN,NMAX,BIG,JBIG)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- use comm_par
- integer :: NB1,NB2,LDVEC,NMIN,NMAX
-!el LOGICAL :: GOPARR,DSKWRK,MASWRK
- real(kind=8),DIMENSION(NB2) :: F
- real(kind=8),DIMENSION(LDVEC,NB1) :: VEC
- real(kind=8),DIMENSION(NB1) :: BIG
- integer,DIMENSION(NB1) :: JBIG
-!
-!el integer :: ME,MASTER,NPROC,IBTYP,IPTIM
-!el COMMON /PAR / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
-!
- real(kind=8),PARAMETER :: ROOT2=0.707106781186548D+00
- real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00, D1050=1.05D+00,&
- D1500=1.5D+00, D3875=3.875D+00,&
- D0500=0.5D+00, D1375=1.375D+00, D0250=0.25D+00
- real(kind=8),PARAMETER :: C2=1.0D-12, C3=4.0D-16,&
- C4=2.0D-16, C5=8.0D-09, C6=3.0D-06
- integer :: i,ii,j,k,jj,IEAA,IEAB,MAXIT,ITER,I1,IA,IB,IAA,IBB,IEAR,&
- IEBR,IR,IT,KQ,IR1
- real(kind=8) :: T,TT,EPS,SD,TEST,DIF,CX,SX,T2X2,T2X25,T1,T2
-!
-! F IS THE MATRIX TO BE DIAGONALIZED, F IS STORED TRIANGULAR
-! VEC IS THE ARRAY OF EIGENVECTORS, DIMENSION NB1*NB1
-! BIG AND JBIG ARE TEMPORARY SCRATCH AREAS OF DIMENSION NB1
-! THE ROTATIONS AMONG THE FIRST NMIN BASIS FUNCTIONS ARE NOT
-! ACCOUNTED FOR.
-! THE ROTATIONS AMONG THE LAST NB1-NMAX BASIS FUNCTIONS ARE NOT
-! ACCOUNTED FOR.
-!
- IEAA=0
- IEAB=0
- TT=ZERO
- EPS = 64.0D+00*EPSLON(ONE)
-!
-! LOOP OVER COLUMNS (K) OF TRIANGULAR MATRIX TO DETERMINE
-! LARGEST OFF-DIAGONAL ELEMENTS IN ROW(I).
-!
- DO 20 I=1,NB1
- BIG(I)=ZERO
- JBIG(I)=0
- IF(I.LT.NMIN .OR. I.EQ.1) GO TO 20
- II = (I*I-I)/2
- J=MIN(I-1,NMAX)
- DO 10 K=1,J
- IF(ABS(BIG(I)).GE.ABS(F(II+K))) GO TO 10
- BIG(I)=F(II+K)
- JBIG(I)=K
- 10 CONTINUE
- 20 CONTINUE
-!
-! ----- 2X2 JACOBI ITERATIONS BEGIN HERE -----
-!
- MAXIT=MAX(NB2*20,500)
- ITER=0
- 30 CONTINUE
- ITER=ITER+1
-!
-! FIND SMALLEST DIAGONAL ELEMENT
-!
- SD=D1050
- JJ=0
- DO 40 J=1,NB1
- JJ=JJ+J
- SD= MIN(SD,ABS(F(JJ)))
- 40 CONTINUE
- TEST = MAX(EPS, C2*MAX(SD,C6))
-!
-! FIND LARGEST OFF-DIAGONAL ELEMENT
-!
- T=ZERO
- I1=MAX(2,NMIN)
- IB = I1
- DO 50 I=I1,NB1
- IF(T.GE.ABS(BIG(I))) GO TO 50
- T= ABS(BIG(I))
- IB=I
- 50 CONTINUE
-!
-! TEST FOR CONVERGENCE, THEN DETERMINE ROTATION.
-!
- IF(T.LT.TEST) RETURN
-! ******
-!
- IF(ITER.GT.MAXIT) THEN
- IF (MASWRK) THEN
- WRITE(6,*) 'JACOBI DIAGONALIZATION FAILS, DIMENSION=',NB1
- WRITE(6,9020) ITER,T,TEST,SD
- ENDIF
- CALL ABRT
- STOP
- END IF
-!
- IA=JBIG(IB)
- IAA=IA*(IA-1)/2
- IBB=IB*(IB-1)/2
- DIF=F(IAA+IA)-F(IBB+IB)
- IF(ABS(DIF).GT.C3*T) GO TO 70
- SX=ROOT2
- CX=ROOT2
- GO TO 110
- 70 T2X2=BIG(IB)/DIF
- T2X25=T2X2*T2X2
- IF(T2X25 .GT. C4) GO TO 80
- CX=ONE
- SX=T2X2
- GO TO 110
- 80 IF(T2X25 .GT. C5) GO TO 90
- SX=T2X2*(ONE-D1500*T2X25)
- CX=ONE-D0500*T2X25
- GO TO 110
- 90 IF(T2X25 .GT. C6) GO TO 100
- CX=ONE+T2X25*(T2X25*D1375 - D0500)
- SX= T2X2*(ONE + T2X25*(T2X25*D3875 - D1500))
- GO TO 110
- 100 T=D0250 / SQRT(D0250 + T2X25)
- CX= SQRT(D0500 + T)
- SX= SIGN( SQRT(D0500 - T),T2X2)
- 110 IEAR=IAA+1
- IEBR=IBB+1
-!
- DO 230 IR=1,NB1
- T=F(IEAR)*SX
- F(IEAR)=F(IEAR)*CX+F(IEBR)*SX
- F(IEBR)=T-F(IEBR)*CX
- IF(IR-IA) 220,120,130
- 120 TT=F(IEBR)
- IEAA=IEAR
- IEAB=IEBR
- F(IEBR)=BIG(IB)
- IEAR=IEAR+IR-1
- IF(JBIG(IR)) 200,220,200
- 130 T=F(IEAR)
- IT=IA
- IEAR=IEAR+IR-1
- IF(IR-IB) 180,150,160
- 150 F(IEAA)=F(IEAA)*CX+F(IEAB)*SX
- F(IEAB)=TT*CX+F(IEBR)*SX
- F(IEBR)=TT*SX-F(IEBR)*CX
- IEBR=IEBR+IR-1
- GO TO 200
- 160 IF( ABS(T) .GE. ABS(F(IEBR))) GO TO 170
- IF(IB.GT.NMAX) GO TO 170
- T=F(IEBR)
- IT=IB
- 170 IEBR=IEBR+IR-1
- 180 IF( ABS(T) .LT. ABS(BIG(IR))) GO TO 190
- BIG(IR) = T
- JBIG(IR) = IT
- GO TO 220
- 190 IF(IA .NE. JBIG(IR) .AND. IB .NE. JBIG(IR)) GO TO 220
- 200 KQ=IEAR-IR-IA+1
- BIG(IR)=ZERO
- IR1=MIN(IR-1,NMAX)
- DO 210 I=1,IR1
- K=KQ+I
- IF(ABS(BIG(IR)) .GE. ABS(F(K))) GO TO 210
- BIG(IR) = F(K)
- JBIG(IR)=I
- 210 CONTINUE
- 220 IEAR=IEAR+1
- 230 IEBR=IEBR+1
-!
- DO 240 I=1,NB1
- T1=VEC(I,IA)*CX + VEC(I,IB)*SX
- T2=VEC(I,IA)*SX - VEC(I,IB)*CX
- VEC(I,IA)=T1
- VEC(I,IB)=T2
- 240 CONTINUE
- GO TO 30
-!
- 9020 FORMAT(1X,'ITER=',I6,' T,TEST,SD=',1P,3E20.10)
- end subroutine JACDIA
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK JACORD
- subroutine JACORD(VEC,EIG,N,LDVEC)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,LDVEC
- real(kind=8),DIMENSION(LDVEC,N) :: VEC
- real(kind=8),DIMENSION(N) :: EIG
- integer :: i,jj,j
- real(kind=8) :: T
-!
-! ---- SORT EIGENDATA INTO ASCENDING ORDER -----
-!
- DO 290 I = 1, N
- JJ = I
- DO 270 J = I, N
- IF (EIG(J) .LT. EIG(JJ)) JJ = J
- 270 CONTINUE
- IF (JJ .EQ. I) GO TO 290
- T = EIG(JJ)
- EIG(JJ) = EIG(I)
- EIG(I) = T
- DO 280 J = 1, N
- T = VEC(J,JJ)
- VEC(J,JJ) = VEC(J,I)
- VEC(J,I) = T
- 280 CONTINUE
- 290 CONTINUE
- return
- end subroutine JACORD
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK TINVTB
- subroutine TINVTB(NM,N,D,E,E2,M,W,IND,Z,IERR,RV1,RV2,RV3,RV4,RV6)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: NM,N,M,IERR
- real(kind=8),DIMENSION(N) :: D,E,E2
- real(kind=8),DIMENSION(M) :: W
- real(kind=8),DIMENSION(NM,M) :: Z
- real(kind=8),DIMENSION(N) :: RV1,RV2,RV3,RV4,RV6
- integer,DIMENSION(M) :: IND
- real(kind=8) :: MACHEP,NORM
- INTEGER :: P,Q,R,S,TAG,GROUP
- integer :: ii,j,jj,i,iqmp,its
- real(kind=8) :: ORDER,XU,UK,X0,U,EPS2,EPS3,EPS4,x1,v
-
-! ------------------------------------------------------------------
-!
-! THIS ROUTINE IS A TRANSLATION OF THE INVERSE ITERATION TECH-
-! NIQUE IN THE ALGOL PROCEDURE TRISTURM BY PETERS AND WILKINSON.
-! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 418-439(1971).
-!
-! THIS ROUTINE FINDS THOSE EIGENVECTORS OF A TRIDIAGONAL
-! SYMMETRIC MATRIX CORRESPONDING TO SPECIFIED EIGENVALUES,
-! USING INVERSE ITERATION.
-!
-! ON INPUT-
-!
-! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
-! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
-! DIMENSION STATEMENT,
-!
-! N IS THE ORDER OF THE MATRIX,
-!
-! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
-!
-! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
-! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
-!
-! E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E,
-! WITH ZEROS CORRESPONDING TO NEGLIGIBLE ELEMENTS OF E.
-! E(I) IS CONSIDERED NEGLIGIBLE IF IT IS NOT LARGER THAN
-! THE PRODUCT OF THE RELATIVE MACHINE PRECISION AND THE SUM
-! OF THE MAGNITUDES OF D(I) AND D(I-1). E2(1) MUST CONTAIN
-! 0.0 IF THE EIGENVALUES ARE IN ASCENDING ORDER, OR 2.0
-! IF THE EIGENVALUES ARE IN DESCENDING ORDER. IF BISECT,
-! TRIDIB, OR IMTQLV HAS BEEN USED TO FIND THE EIGENVALUES,
-! THEIR OUTPUT E2 ARRAY IS EXACTLY WHAT IS EXPECTED HERE,
-!
-! M IS THE NUMBER OF SPECIFIED EIGENVALUES,
-!
-! W CONTAINS THE M EIGENVALUES IN ASCENDING OR DESCENDING ORDER,
-!
-! IND CONTAINS IN ITS FIRST M POSITIONS THE SUBMATRIX INDICES
-! ASSOCIATED WITH THE CORRESPONDING EIGENVALUES IN W --
-! 1 FOR EIGENVALUES BELONGING TO THE FIRST SUBMATRIX FROM
-! THE TOP, 2 FOR THOSE BELONGING TO THE SECOND SUBMATRIX, ETC.
-!
-! ON OUTPUT-
-!
-! ALL INPUT ARRAYS ARE UNALTERED,
-!
-! Z CONTAINS THE ASSOCIATED SET OF ORTHONORMAL EIGENVECTORS.
-! ANY VECTOR WHICH FAILS TO CONVERGE IS SET TO ZERO,
-!
-! IERR IS SET TO
-! ZERO FOR NORMAL RETURN,
-! -R IF THE EIGENVECTOR CORRESPONDING TO THE R-TH
-! EIGENVALUE FAILS TO CONVERGE IN 5 ITERATIONS,
-!
-! RV1, RV2, RV3, RV4, AND RV6 ARE TEMPORARY STORAGE ARRAYS.
-!
-! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
-! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
-!
-! ------------------------------------------------------------------
-!
-! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
-! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
-!
-! **********
- MACHEP = 2.0D+00**(-50)
-!
- IERR = 0
- IF (M .EQ. 0) GO TO 680
- TAG = 0
- ORDER = 1.0D+00 - E2(1)
- XU = 0.0D+00
- UK = 0.0D+00
- X0 = 0.0D+00
- U = 0.0D+00
- EPS2 = 0.0D+00
- EPS3 = 0.0D+00
- EPS4 = 0.0D+00
- GROUP = 0
- Q = 0
-! ********** ESTABLISH AND PROCESS NEXT SUBMATRIX **********
- 100 P = Q + 1
- IP = P + 1
-!
- DO 120 Q = P, N
- IF (Q .EQ. N) GO TO 140
- IF (E2(Q+1) .EQ. 0.0D+00) GO TO 140
- 120 CONTINUE
-! ********** FIND VECTORS BY INVERSE ITERATION **********
- 140 TAG = TAG + 1
- IQMP = Q - P + 1
- S = 0
-!
- DO 660 R = 1, M
- IF (IND(R) .NE. TAG) GO TO 660
- ITS = 1
- X1 = W(R)
- IF (S .NE. 0) GO TO 220
-! ********** CHECK FOR ISOLATED ROOT **********
- XU = 1.0D+00
- IF (P .NE. Q) GO TO 160
- RV6(P) = 1.0D+00
- GO TO 600
- 160 NORM = ABS(D(P))
-!
- DO 180 I = IP, Q
- 180 NORM = NORM + ABS(D(I)) + ABS(E(I))
-! ********** EPS2 IS THE CRITERION FOR GROUPING,
-! EPS3 REPLACES ZERO PIVOTS AND EQUAL
-! ROOTS ARE MODIFIED BY EPS3,
-! EPS4 IS TAKEN VERY SMALL TO AVOID OVERFLOW **********
- EPS2 = 1.0D-03 * NORM
- EPS3 = MACHEP * NORM
- UK = IQMP
- EPS4 = UK * EPS3
- UK = EPS4 / SQRT(UK)
- S = P
- 200 GROUP = 0
- GO TO 240
-! ********** LOOK FOR CLOSE OR COINCIDENT ROOTS **********
- 220 IF (ABS(X1-X0) .GE. EPS2) GO TO 200
- GROUP = GROUP + 1
- IF (ORDER * (X1 - X0) .LE. 0.0D+00) X1 = X0 + ORDER * EPS3
-! ********** ELIMINATION WITH INTERCHANGES AND
-! INITIALIZATION OF VECTOR **********
- 240 V = 0.0D+00
-!
- DO 300 I = P, Q
- RV6(I) = UK
- IF (I .EQ. P) GO TO 280
- IF (ABS(E(I)) .LT. ABS(U)) GO TO 260
-! ********** WARNING -- A DIVIDE CHECK MAY OCCUR HERE IF
-! E2 ARRAY HAS NOT BEEN SPECIFIED CORRECTLY **********
- XU = U / E(I)
- RV4(I) = XU
- RV1(I-1) = E(I)
- RV2(I-1) = D(I) - X1
- RV3(I-1) = 0.0D+00
- IF (I .NE. Q) RV3(I-1) = E(I+1)
- U = V - XU * RV2(I-1)
- V = -XU * RV3(I-1)
- GO TO 300
- 260 XU = E(I) / U
- RV4(I) = XU
- RV1(I-1) = U
- RV2(I-1) = V
- RV3(I-1) = 0.0D+00
- 280 U = D(I) - X1 - XU * V
- IF (I .NE. Q) V = E(I+1)
- 300 CONTINUE
-!
- IF (U .EQ. 0.0D+00) U = EPS3
- RV1(Q) = U
- RV2(Q) = 0.0D+00
- RV3(Q) = 0.0D+00
-! ********** BACK SUBSTITUTION
-! FOR I=Q STEP -1 UNTIL P DO -- **********
- 320 DO 340 II = P, Q
- I = P + Q - II
- RV6(I) = (RV6(I) - U * RV2(I) - V * RV3(I)) / RV1(I)
- V = U
- U = RV6(I)
- 340 CONTINUE
-! ********** ORTHOGONALIZE WITH RESPECT TO PREVIOUS
-! MEMBERS OF GROUP **********
- IF (GROUP .EQ. 0) GO TO 400
- J = R
-!
- DO 380 JJ = 1, GROUP
- 360 J = J - 1
- IF (IND(J) .NE. TAG) GO TO 360
- XU = DDOT(IQMP,RV6(P),1,Z(P,J),1)
-!
- CALL DAXPY(IQMP,-XU,Z(P,J),1,RV6(P),1)
-!
- 380 CONTINUE
-!
- 400 NORM = 0.0D+00
-!
- DO 420 I = P, Q
- 420 NORM = NORM + ABS(RV6(I))
-!
- IF (NORM .GE. 1.0D+00) GO TO 560
-! ********** FORWARD SUBSTITUTION **********
- IF (ITS .EQ. 5) GO TO 540
- IF (NORM .NE. 0.0D+00) GO TO 440
- RV6(S) = EPS4
- S = S + 1
- IF (S .GT. Q) S = P
- GO TO 480
- 440 XU = EPS4 / NORM
-!
- DO 460 I = P, Q
- 460 RV6(I) = RV6(I) * XU
-! ********** ELIMINATION OPERATIONS ON NEXT VECTOR
-! ITERATE **********
- 480 DO 520 I = IP, Q
- U = RV6(I)
-! ********** IF RV1(I-1) .EQ. E(I), A ROW INTERCHANGE
-! WAS PERFORMED EARLIER IN THE
-! TRIANGULARIZATION PROCESS **********
- IF (RV1(I-1) .NE. E(I)) GO TO 500
- U = RV6(I-1)
- RV6(I-1) = RV6(I)
- 500 RV6(I) = U - RV4(I) * RV6(I-1)
- 520 CONTINUE
-!
- ITS = ITS + 1
- GO TO 320
-! ********** SET ERROR -- NON-CONVERGED EIGENVECTOR **********
- 540 IERR = -R
- XU = 0.0D+00
- GO TO 600
-! ********** NORMALIZE SO THAT SUM OF SQUARES IS
-! 1 AND EXPAND TO FULL ORDER **********
- 560 U = 0.0D+00
-!
- DO 580 I = P, Q
- RV6(I) = RV6(I) / NORM
- 580 U = U + RV6(I)**2
-!
- XU = 1.0D+00 / SQRT(U)
-!
- 600 DO 620 I = 1, N
- 620 Z(I,R) = 0.0D+00
-!
- DO 640 I = P, Q
- 640 Z(I,R) = RV6(I) * XU
-!
- X0 = X1
- 660 CONTINUE
-!
- IF (Q .LT. N) GO TO 100
- 680 return
-! ********** LAST CARD OF TINVIT **********
- end subroutine TINVTB
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK TQL2
- subroutine TQL2(NM,N,D,E,Z,IERR)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: NM,N,IERR
- real(kind=8) :: MACHEP
- real(kind=8),DIMENSION(N) :: D!,E
- real(kind=8) :: E(*)!el E(L)
- real(kind=8),DIMENSION(NM,N) :: Z
- integer :: ii,i,j,mml,m,l1,k,l
- real(kind=8) :: c,f,b,h,g,p,r,s
-!
-! THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2,
-! NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND
-! WILKINSON.
-! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971).
-!
-! THIS ROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
-! OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE QL METHOD.
-! THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO
-! BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS
-! FULL MATRIX TO TRIDIAGONAL FORM.
-!
-! ON INPUT-
-!
-! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
-! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
-! DIMENSION STATEMENT,
-!
-! N IS THE ORDER OF THE MATRIX,
-!
-! D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX,
-!
-! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
-! IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY,
-!
-! Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
-! REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS
-! OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN
-! THE IDENTITY MATRIX.
-!
-! ON OUTPUT-
-!
-! D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN
-! ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT
-! UNORDERED FOR INDICES 1,2,...,IERR-1,
-!
-! E HAS BEEN DESTROYED,
-!
-! Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC
-! TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE,
-! Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED
-! EIGENVALUES,
-!
-! IERR IS SET TO
-! ZERO FOR NORMAL RETURN,
-! J IF THE J-TH EIGENVALUE HAS NOT BEEN
-! DETERMINED AFTER 30 ITERATIONS.
-!
-! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
-! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
-!
-! ------------------------------------------------------------------
-!
-! ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
-! THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.
-!
-! **********
- MACHEP = 2.0D+00**(-50)
-!
- IERR = 0
- IF (N .EQ. 1) GO TO 400
-!
- DO 100 I = 2, N
- 100 E(I-1) = E(I)
-!
- F = 0.0D+00
- B = 0.0D+00
- E(N) = 0.0D+00
-!
- DO 300 L = 1, N
- J = 0
- H = MACHEP * (ABS(D(L)) + ABS(E(L)))
- IF (B .LT. H) B = H
-! ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT **********
- DO 120 M = L, N
- IF (ABS(E(M)) .LE. B) GO TO 140
-! ********** E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT
-! THROUGH THE BOTTOM OF THE LOOP **********
- 120 CONTINUE
-!
- 140 IF (M .EQ. L) GO TO 280
- 160 IF (J .EQ. 30) GO TO 380
- J = J + 1
-! ********** FORM SHIFT **********
- L1 = L + 1
- G = D(L)
- P = (D(L1) - G) / (2.0D+00 * E(L))
- R = SQRT(P*P+1.0D+00)
- D(L) = E(L) / (P + SIGN(R,P))
- H = G - D(L)
-!
- DO 180 I = L1, N
- 180 D(I) = D(I) - H
-!
- F = F + H
-! ********** QL TRANSFORMATION **********
- P = D(M)
- C = 1.0D+00
- S = 0.0D+00
- MML = M - L
-! ********** FOR I=M-1 STEP -1 UNTIL L DO -- **********
- DO 260 II = 1, MML
- I = M - II
- G = C * E(I)
- H = C * P
- IF (ABS(P) .LT. ABS(E(I))) GO TO 200
- C = E(I) / P
- R = SQRT(C*C+1.0D+00)
- E(I+1) = S * P * R
- S = C / R
- C = 1.0D+00 / R
- GO TO 220
- 200 C = P / E(I)
- R = SQRT(C*C+1.0D+00)
- E(I+1) = S * E(I) * R
- S = 1.0D+00 / R
- C = C * S
- 220 P = C * D(I) - S * G
- D(I+1) = H + S * (C * G + S * D(I))
-! ********** FORM VECTOR **********
- CALL DROT(N,Z(1,I+1),1,Z(1,I),1,C,S)
-!
- 260 CONTINUE
-!
- E(L) = S * P
- D(L) = C * P
- IF (ABS(E(L)) .GT. B) GO TO 160
- 280 D(L) = D(L) + F
- 300 CONTINUE
-! ********** ORDER EIGENVALUES AND EIGENVECTORS **********
- DO 360 II = 2, N
- I = II - 1
- K = I
- P = D(I)
-!
- DO 320 J = II, N
- IF (D(J) .GE. P) GO TO 320
- K = J
- P = D(J)
- 320 CONTINUE
-!
- IF (K .EQ. I) GO TO 360
- D(K) = D(I)
- D(I) = P
-!
- CALL DSWAP(N,Z(1,I),1,Z(1,K),1)
-!
- 360 CONTINUE
-!
- GO TO 400
-! ********** SET ERROR -- NO CONVERGENCE TO AN
-! EIGENVALUE AFTER 30 ITERATIONS **********
- 380 IERR = L
- 400 return
-! ********** LAST CARD OF TQL2 **********
- end subroutine TQL2
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK TRBK3B
- subroutine TRBK3B(NM,N,NV,A,M,Z)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: NM,N,NV,M
- real(kind=8),DIMENSION(NV) :: A
- real(kind=8),DIMENSION(NM,M) :: Z
- integer :: i,l,iz,ik,j
- real(kind=8) :: h,s
-!
-! THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRBAK3,
-! NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
-! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
-!
-! THIS ROUTINE FORMS THE EIGENVECTORS OF A REAL SYMMETRIC
-! MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
-! SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY TRED3B.
-!
-! ON INPUT-
-!
-! NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL
-! ARRAY PARAMETERS AS DECLARED IN THE CALLING ROUTINE
-! DIMENSION STATEMENT,
-!
-! N IS THE ORDER OF THE MATRIX,
-!
-! NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
-! AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
-!
-! A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANSFORMATIONS
-! USED IN THE REDUCTION BY TRED3B IN ITS FIRST
-! N*(N+1)/2 POSITIONS,
-!
-! M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED,
-!
-! Z CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
-! IN ITS FIRST M COLUMNS.
-!
-! ON OUTPUT-
-!
-! Z CONTAINS THE TRANSFORMED EIGENVECTORS
-! IN ITS FIRST M COLUMNS.
-!
-! NOTE THAT TRBAK3 PRESERVES VECTOR EUCLIDEAN NORMS.
-!
-! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
-! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
-!
-! ------------------------------------------------------------------
-!
- IF (M .EQ. 0) GO TO 140
- IF (N .EQ. 1) GO TO 140
-!
- DO 120 I = 2, N
- L = I - 1
- IZ = (I * L) / 2
- IK = IZ + I
- H = A(IK)
- IF (H .EQ. 0.0D+00) GO TO 120
-!
- DO 100 J = 1, M
- S = -DDOT(L,A(IZ+1),1,Z(1,J),1)
-!
-! ********** DOUBLE DIVISION AVOIDS POSSIBLE UNDERFLOW **********
- S = (S / H) / H
-!
- CALL DAXPY(L,S,A(IZ+1),1,Z(1,J),1)
-!
- 100 CONTINUE
-!
- 120 CONTINUE
-!
- 140 return
-! ********** LAST CARD OF TRBAK3 **********
- end subroutine TRBK3B
-!-----------------------------------------------------------------------------
-!*MODULE EIGEN *DECK TRED3B
- subroutine TRED3B(N,NV,A,D,E,E2)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,NV
- real(kind=8),DIMENSION(NV) :: A
- real(kind=8),DIMENSION(N) :: D,E2
- real(kind=8) :: E(*)!el E(L)
- integer :: ii,i,l,iz,k,jk,j,jm1
- real(kind=8) :: h,f,g,scale,dt,hh
-!
-! THIS ROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED3,
-! NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON.
-! HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
-!
-! THIS ROUTINE REDUCES A REAL SYMMETRIC MATRIX, STORED AS
-! A ONE-DIMENSIONAL ARRAY, TO A SYMMETRIC TRIDIAGONAL MATRIX
-! USING ORTHOGONAL SIMILARITY TRANSFORMATIONS.
-!
-! ON INPUT-
-!
-! N IS THE ORDER OF THE MATRIX,
-!
-! NV MUST BE SET TO THE DIMENSION OF THE ARRAY PARAMETER A
-! AS DECLARED IN THE CALLING ROUTINE DIMENSION STATEMENT,
-!
-! A CONTAINS THE LOWER TRIANGLE OF THE REAL SYMMETRIC
-! INPUT MATRIX, STORED ROW-WISE AS A ONE-DIMENSIONAL
-! ARRAY, IN ITS FIRST N*(N+1)/2 POSITIONS.
-!
-! ON OUTPUT-
-!
-! A CONTAINS INFORMATION ABOUT THE ORTHOGONAL
-! TRANSFORMATIONS USED IN THE REDUCTION,
-!
-! D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX,
-!
-! E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL
-! MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO,
-!
-! E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E.
-! E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED.
-!
-! QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,
-! APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY
-!
-! ------------------------------------------------------------------
-!
-! ********** FOR I=N STEP -1 UNTIL 1 DO -- **********
- DO 300 II = 1, N
- I = N + 1 - II
- L = I - 1
- IZ = (I * L) / 2
- H = 0.0D+00
- SCALE = 0.0D+00
- IF (L .LT. 1) GO TO 120
-! ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) **********
- DO 100 K = 1, L
- IZ = IZ + 1
- D(K) = A(IZ)
- SCALE = SCALE + ABS(D(K))
- 100 CONTINUE
-!
- IF (SCALE .NE. 0.0D+00) GO TO 140
- 120 E(I) = 0.0D+00
- E2(I) = 0.0D+00
- GO TO 280
-!
- 140 DO 160 K = 1, L
- D(K) = D(K) / SCALE
- H = H + D(K) * D(K)
- 160 CONTINUE
-!
- E2(I) = SCALE * SCALE * H
- F = D(L)
- G = -SIGN(SQRT(H),F)
- E(I) = SCALE * G
- H = H - F * G
- D(L) = F - G
- A(IZ) = SCALE * D(L)
- IF (L .EQ. 1) GO TO 280
- F = 0.0D+00
-!
- JK = 1
- DO 220 J = 1, L
- JM1 = J - 1
- DT = D(J)
- G = 0.0D+00
-! ********** FORM ELEMENT OF A*U **********
- IF (JM1 .EQ. 0) GO TO 200
- DO 180 K = 1, JM1
- E(K) = E(K) + DT * A(JK)
- G = G + D(K) * A(JK)
- JK = JK + 1
- 180 CONTINUE
- 200 E(J) = G + A(JK) * DT
- JK = JK + 1
-! ********** FORM ELEMENT OF P **********
- 220 CONTINUE
- F = 0.0D+00
- DO 240 J = 1, L
- E(J) = E(J) / H
- F = F + E(J) * D(J)
- 240 CONTINUE
-!
- HH = F / (H + H)
- JK = 0
-! ********** FORM REDUCED A **********
- DO 260 J = 1, L
- F = D(J)
- G = E(J) - HH * F
- E(J) = G
-!
- DO 260 K = 1, J
- JK = JK + 1
- A(JK) = A(JK) - F * E(K) - G * D(K)
- 260 CONTINUE
-!
- 280 D(I) = A(IZ+1)
- A(IZ+1) = SCALE * SQRT(H)
- 300 CONTINUE
-!
- return
-! ********** LAST CARD OF TRED3 **********
- end subroutine TRED3B
-!-----------------------------------------------------------------------------
-! blas.f
-!-----------------------------------------------------------------------------
-! 10 NOV 94 - MWS - DNRM2: REMOVE FTNCHECK WARNINGS
-! 11 JUN 94 - MWS - INCLUDE A COPY OF DGEMV (LEVEL TWO ROUTINE)
-! 11 AUG 87 - MWS - SANITIZE FLOATING POINT CONSTANTS IN DNRM2
-! 26 MAR 87 - MWS - USE GENERIC SIGN IN DROTG
-! 28 NOV 86 - STE - SUPPLY ALL LEVEL ONE BLAS
-! 7 JUL 86 - JAB - SANITIZE FLOATING POINT CONSTANTS
-!
-! BASIC LINEAR ALGEBRA SUBPROGRAMS (BLAS) FROM LINPACK (LEVEL 1)
-!
-! THIS MODULE SHOULD BE COMPILED ONLY IF SPECIALLY CODED
-! VERSIONS OF THESE ROUTINES ARE NOT AVAILABLE ON THE TARGET MACHINE
-!
-!*MODULE BLAS1 *DECK DASUM
- real(kind=8) function DASUM(N,DX,INCX)
-!
-! TAKES THE SUM OF THE ABSOLUTE VALUES.
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- real(kind=8) :: DX(1),DTEMP
- INTEGER :: I,INCX,M,MP1,N,NINCX
-!
- DASUM = 0.0D+00
- DTEMP = 0.0D+00
- IF(N.LE.0)RETURN
- IF(INCX.EQ.1)GO TO 20
-!
-! CODE FOR INCREMENT NOT EQUAL TO 1
-!
- NINCX = N*INCX
- DO 10 I = 1,NINCX,INCX
- DTEMP = DTEMP + ABS(DX(I))
- 10 CONTINUE
- DASUM = DTEMP
- RETURN
-!
-! CODE FOR INCREMENT EQUAL TO 1
-!
-!
-! CLEAN-UP LOOP
-!
- 20 M = MOD(N,6)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DTEMP = DTEMP + ABS(DX(I))
- 30 CONTINUE
- IF( N .LT. 6 ) GO TO 60
- 40 MP1 = M + 1
- DO 50 I = MP1,N,6
- DTEMP = DTEMP + ABS(DX(I)) + ABS(DX(I + 1)) + ABS(DX(I + 2)) &
- + ABS(DX(I + 3)) + ABS(DX(I + 4)) + ABS(DX(I + 5))
- 50 CONTINUE
- 60 DASUM = DTEMP
- return
- end function DASUM
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK DAXPY
- subroutine DAXPY(N,DA,DX,INCX,DY,INCY)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,INCX,INCY
- real(kind=8),DIMENSION(1) :: DX,DY
- real(kind=8) :: DA
-!
-! CONSTANT TIMES A VECTOR PLUS A VECTOR.
-! DY(I) = DY(I) + DA * DX(I)
-! USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- integer :: ix,iy,i,m,mp1
- IF(N.LE.0)RETURN
- IF (DA .EQ. 0.0D+00) RETURN
- IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
-!
-! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
-! NOT EQUAL TO 1
-!
- IX = 1
- IY = 1
- IF(INCX.LT.0)IX = (-N+1)*INCX + 1
- IF(INCY.LT.0)IY = (-N+1)*INCY + 1
- DO 10 I = 1,N
- DY(IY) = DY(IY) + DA*DX(IX)
- IX = IX + INCX
- IY = IY + INCY
- 10 CONTINUE
- RETURN
-!
-! CODE FOR BOTH INCREMENTS EQUAL TO 1
-!
-!
-! CLEAN-UP LOOP
-!
- 20 M = MOD(N,4)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DY(I) = DY(I) + DA*DX(I)
- 30 CONTINUE
- IF( N .LT. 4 ) RETURN
- 40 MP1 = M + 1
- DO 50 I = MP1,N,4
- DY(I) = DY(I) + DA*DX(I)
- DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
- DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
- DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
- 50 CONTINUE
- return
- end subroutine DAXPY
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK DCOPY
- subroutine DCOPY(N,DX,INCX,DY,INCY)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,INCX,INCY
- real(kind=8),DIMENSION(*) :: DX,DY
-!
-! COPIES A VECTOR.
-! DY(I) <== DX(I)
-! USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- integer :: ix,iy,m,i,mp1
- IF(N.LE.0)RETURN
- IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
-!
-! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
-! NOT EQUAL TO 1
-!
- IX = 1
- IY = 1
- IF(INCX.LT.0)IX = (-N+1)*INCX + 1
- IF(INCY.LT.0)IY = (-N+1)*INCY + 1
- DO 10 I = 1,N
- DY(IY) = DX(IX)
- IX = IX + INCX
- IY = IY + INCY
- 10 CONTINUE
- RETURN
-!
-! CODE FOR BOTH INCREMENTS EQUAL TO 1
-!
-!
-! CLEAN-UP LOOP
-!
- 20 M = MOD(N,7)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DY(I) = DX(I)
- 30 CONTINUE
- IF( N .LT. 7 ) RETURN
- 40 MP1 = M + 1
- DO 50 I = MP1,N,7
- DY(I) = DX(I)
- DY(I + 1) = DX(I + 1)
- DY(I + 2) = DX(I + 2)
- DY(I + 3) = DX(I + 3)
- DY(I + 4) = DX(I + 4)
- DY(I + 5) = DX(I + 5)
- DY(I + 6) = DX(I + 6)
- 50 CONTINUE
- return
- end subroutine DCOPY
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK DDOT
- real(kind=8) function DDOT(N,DX,INCX,DY,INCY)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,INCX,INCY
- real(kind=8),DIMENSION(1) :: DX,DY
-!
-! FORMS THE DOT PRODUCT OF TWO VECTORS.
-! DOT = DX(I) * DY(I)
-! USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- integer ::ix,iy,m,mp1,i
- real(kind=8) :: DTEMP
- DDOT = 0.0D+00
- DTEMP = 0.0D+00
- IF(N.LE.0)RETURN
- IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
-!
-! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
-! NOT EQUAL TO 1
-!
- IX = 1
- IY = 1
- IF(INCX.LT.0)IX = (-N+1)*INCX + 1
- IF(INCY.LT.0)IY = (-N+1)*INCY + 1
- DO 10 I = 1,N
- DTEMP = DTEMP + DX(IX)*DY(IY)
- IX = IX + INCX
- IY = IY + INCY
- 10 CONTINUE
- DDOT = DTEMP
- RETURN
-!
-! CODE FOR BOTH INCREMENTS EQUAL TO 1
-!
-!
-! CLEAN-UP LOOP
-!
- 20 M = MOD(N,5)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DTEMP = DTEMP + DX(I)*DY(I)
- 30 CONTINUE
- IF( N .LT. 5 ) GO TO 60
- 40 MP1 = M + 1
- DO 50 I = MP1,N,5
- DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + &
- DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4)
- 50 CONTINUE
- 60 DDOT = DTEMP
- return
- end function DDOT
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK DNRM2
- real(kind=8) function DNRM2(N,DX,INCX)
-
- INTEGER :: NEXT,N,INCX
- real(kind=8) :: DX(1), CUTLO, CUTHI, HITEST, SUM, XMAX,ZERO,ONE
- DATA ZERO, ONE /0.0D+00, 1.0D+00/
-
- integer :: i,j,nn
-!
-! EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE
-! INCREMENT INCX .
-! IF N .LE. 0 RETURN WITH RESULT = 0.
-! IF N .GE. 1 THEN INCX MUST BE .GE. 1
-!
-! C.L.LAWSON, 1978 JAN 08
-!
-! FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE
-! HOPEFULLY APPLICABLE TO ALL MACHINES.
-! CUTLO = MAXIMUM OF SQRT(U/EPS) OVER ALL KNOWN MACHINES.
-! CUTHI = MINIMUM OF SQRT(V) OVER ALL KNOWN MACHINES.
-! WHERE
-! EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1.
-! U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT)
-! V = LARGEST NO. (OVERFLOW LIMIT)
-!
-! BRIEF OUTLINE OF ALGORITHM..
-!
-! PHASE 1 SCANS ZERO COMPONENTS.
-! MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO
-! MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO
-! MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M
-! WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX.
-!
-! VALUES FOR CUTLO AND CUTHI..
-! FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER
-! DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS..
-! CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE
-! UNIVAC AND DEC AT 2**(-103)
-! THUS CUTLO = 2**(-51) = 4.44089E-16
-! CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC.
-! THUS CUTHI = 2**(63.5) = 1.30438E19
-! CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC.
-! THUS CUTLO = 2**(-33.5) = 8.23181D-11
-! CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D+19
-! DATA CUTLO, CUTHI / 8.232D-11, 1.304D+19 /
-! DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 /
- DATA CUTLO, CUTHI / 8.232D-11, 1.304D+19 /
-!
- J=0
- IF(N .GT. 0) GO TO 10
- DNRM2 = ZERO
- GO TO 300
-!
- 10 ASSIGN 30 TO NEXT
- SUM = ZERO
- NN = N * INCX
-! BEGIN MAIN LOOP
- I = 1
- 20 GO TO NEXT,(30, 50, 70, 110)
- 30 IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
- ASSIGN 50 TO NEXT
- XMAX = ZERO
-!
-! PHASE 1. SUM IS ZERO
-!
- 50 IF( DX(I) .EQ. ZERO) GO TO 200
- IF( ABS(DX(I)) .GT. CUTLO) GO TO 85
-!
-! PREPARE FOR PHASE 2.
- ASSIGN 70 TO NEXT
- GO TO 105
-!
-! PREPARE FOR PHASE 4.
-!
- 100 I = J
- ASSIGN 110 TO NEXT
- SUM = (SUM / DX(I)) / DX(I)
- 105 XMAX = ABS(DX(I))
- GO TO 115
-!
-! PHASE 2. SUM IS SMALL.
-! SCALE TO AVOID DESTRUCTIVE UNDERFLOW.
-!
- 70 IF( ABS(DX(I)) .GT. CUTLO ) GO TO 75
-!
-! COMMON CODE FOR PHASES 2 AND 4.
-! IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW.
-!
- 110 IF( ABS(DX(I)) .LE. XMAX ) GO TO 115
- SUM = ONE + SUM * (XMAX / DX(I))**2
- XMAX = ABS(DX(I))
- GO TO 200
-!
- 115 SUM = SUM + (DX(I)/XMAX)**2
- GO TO 200
-!
-!
-! PREPARE FOR PHASE 3.
-!
- 75 SUM = (SUM * XMAX) * XMAX
-!
-!
-! FOR REAL OR D.P. SET HITEST = CUTHI/N
-! FOR COMPLEX SET HITEST = CUTHI/(2*N)
-!
- 85 HITEST = CUTHI/N
-!
-! PHASE 3. SUM IS MID-RANGE. NO SCALING.
-!
- DO 95 J =I,NN,INCX
- IF(ABS(DX(J)) .GE. HITEST) GO TO 100
- 95 SUM = SUM + DX(J)**2
- DNRM2 = SQRT( SUM )
- GO TO 300
-!
- 200 CONTINUE
- I = I + INCX
- IF ( I .LE. NN ) GO TO 20
-!
-! END OF MAIN LOOP.
-!
-! COMPUTE SQUARE ROOT AND ADJUST FOR SCALING.
-!
- DNRM2 = XMAX * SQRT(SUM)
- 300 CONTINUE
- return
- end function DNRM2
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK DROT
- subroutine DROT(N,DX,INCX,DY,INCY,C,S)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,INCX,INCY
- real(kind=8),DIMENSION(1) :: DX,DY
- real(kind=8) :: C,S
-!
-! APPLIES A PLANE ROTATION.
-! DX(I) = C*DX(I) + S*DY(I)
-! DY(I) = -S*DX(I) + C*DY(I)
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- integer :: ix,iy,i
- real(kind=8) :: DTEMP
- IF(N.LE.0)RETURN
- IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
-!
-! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
-! TO 1
-!
- IX = 1
- IY = 1
- IF(INCX.LT.0)IX = (-N+1)*INCX + 1
- IF(INCY.LT.0)IY = (-N+1)*INCY + 1
- DO 10 I = 1,N
- DTEMP = C*DX(IX) + S*DY(IY)
- DY(IY) = C*DY(IY) - S*DX(IX)
- DX(IX) = DTEMP
- IX = IX + INCX
- IY = IY + INCY
- 10 CONTINUE
- RETURN
-!
-! CODE FOR BOTH INCREMENTS EQUAL TO 1
-!
- 20 DO 30 I = 1,N
- DTEMP = C*DX(I) + S*DY(I)
- DY(I) = C*DY(I) - S*DX(I)
- DX(I) = DTEMP
- 30 CONTINUE
- return
- end subroutine DROT
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK DROTG
- subroutine DROTG(DA,DB,C,S)
-!
-! CONSTRUCT GIVENS PLANE ROTATION.
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- real(kind=8) :: DA,DB,C,S,ROE,SCALE,R,Z
-!
- real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00
-!
-!-----------------------------------------------------------------------
-!
-!
- ROE = DB
- IF( ABS(DA) .GT. ABS(DB) ) ROE = DA
- SCALE = ABS(DA) + ABS(DB)
- IF( SCALE .NE. ZERO ) GO TO 10
- C = ONE
- S = ZERO
- R = ZERO
- GO TO 20
-!
- 10 R = SCALE*SQRT((DA/SCALE)**2 + (DB/SCALE)**2)
- R = SIGN(ONE,ROE)*R
- C = DA/R
- S = DB/R
- 20 Z = ONE
- IF( ABS(DA) .GT. ABS(DB) ) Z = S
- IF( ABS(DB) .GE. ABS(DA) .AND. C .NE. ZERO ) Z = ONE/C
- DA = R
- DB = Z
- return
- end subroutine DROTG
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK DSCAL
- subroutine DSCAL(N,DA,DX,INCX)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,INCX
- real(kind=8),DIMENSION(1) :: DX
- real(kind=8) :: DA
-!
-! SCALES A VECTOR BY A CONSTANT.
-! DX(I) = DA * DX(I)
-! USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE.
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- integer :: NINCX,m,mp1,i
- IF(N.LE.0)RETURN
- IF(INCX.EQ.1)GO TO 20
-!
-! CODE FOR INCREMENT NOT EQUAL TO 1
-!
- NINCX = N*INCX
- DO 10 I = 1,NINCX,INCX
- DX(I) = DA*DX(I)
- 10 CONTINUE
- RETURN
-!
-! CODE FOR INCREMENT EQUAL TO 1
-!
-!
-! CLEAN-UP LOOP
-!
- 20 M = MOD(N,5)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DX(I) = DA*DX(I)
- 30 CONTINUE
- IF( N .LT. 5 ) RETURN
- 40 MP1 = M + 1
- DO 50 I = MP1,N,5
- DX(I) = DA*DX(I)
- DX(I + 1) = DA*DX(I + 1)
- DX(I + 2) = DA*DX(I + 2)
- DX(I + 3) = DA*DX(I + 3)
- DX(I + 4) = DA*DX(I + 4)
- 50 CONTINUE
- return
- end subroutine DSCAL
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK DSWAP
- subroutine DSWAP(N,DX,INCX,DY,INCY)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,INCX,INCY
- real(kind=8),DIMENSION(1) :: DX,DY
-!
-! INTERCHANGES TWO VECTORS.
-! DX(I) <==> DY(I)
-! USES UNROLLED LOOPS FOR INCREMENTS EQUAL ONE.
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- integer :: ix,iy,i,m,mp1
- real(kind=8) :: DTEMP
- IF(N.LE.0)RETURN
- IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
-!
-! CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL
-! TO 1
-!
- IX = 1
- IY = 1
- IF(INCX.LT.0)IX = (-N+1)*INCX + 1
- IF(INCY.LT.0)IY = (-N+1)*INCY + 1
- DO 10 I = 1,N
- DTEMP = DX(IX)
- DX(IX) = DY(IY)
- DY(IY) = DTEMP
- IX = IX + INCX
- IY = IY + INCY
- 10 CONTINUE
- RETURN
-!
-! CODE FOR BOTH INCREMENTS EQUAL TO 1
-!
-!
-! CLEAN-UP LOOP
-!
- 20 M = MOD(N,3)
- IF( M .EQ. 0 ) GO TO 40
- DO 30 I = 1,M
- DTEMP = DX(I)
- DX(I) = DY(I)
- DY(I) = DTEMP
- 30 CONTINUE
- IF( N .LT. 3 ) RETURN
- 40 MP1 = M + 1
- DO 50 I = MP1,N,3
- DTEMP = DX(I)
- DX(I) = DY(I)
- DY(I) = DTEMP
- DTEMP = DX(I + 1)
- DX(I + 1) = DY(I + 1)
- DY(I + 1) = DTEMP
- DTEMP = DX(I + 2)
- DX(I + 2) = DY(I + 2)
- DY(I + 2) = DTEMP
- 50 CONTINUE
- return
- end subroutine DSWAP
-!-----------------------------------------------------------------------------
-!*MODULE BLAS1 *DECK IDAMAX
- integer function IDAMAX(N,DX,INCX)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: N,INCX
- real(kind=8),DIMENSION(1) :: DX
-!
-! FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE.
-! JACK DONGARRA, LINPACK, 3/11/78.
-!
- integer :: ix,iy,i
- real(kind=8) :: RMAX
- IDAMAX = 0
- IF( N .LT. 1 ) RETURN
- IDAMAX = 1
- IF(N.EQ.1)RETURN
- IF(INCX.EQ.1)GO TO 20
-!
-! CODE FOR INCREMENT NOT EQUAL TO 1
-!
- IX = 1
- RMAX = ABS(DX(1))
- IX = IX + INCX
- DO 10 I = 2,N
- IF(ABS(DX(IX)).LE.RMAX) GO TO 5
- IDAMAX = I
- RMAX = ABS(DX(IX))
- 5 IX = IX + INCX
- 10 CONTINUE
- RETURN
-!
-! CODE FOR INCREMENT EQUAL TO 1
-!
- 20 RMAX = ABS(DX(1))
- DO 30 I = 2,N
- IF(ABS(DX(I)).LE.RMAX) GO TO 30
- IDAMAX = I
- RMAX = ABS(DX(I))
- 30 CONTINUE
- return
- end function IDAMAX
-!-----------------------------------------------------------------------------
-!*MODULE BLAS *DECK DGEMV
- subroutine DGEMV(FORMA,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
-
-! IMPLICIT DOUBLE PRECISION(A-H,O-Z)
- integer :: M,N,INCX,INCY,LDA
- CHARACTER(len=1) :: FORMA
- real(kind=8),DIMENSION(LDA,*) :: A
- real(kind=8),DIMENSION(*) :: X,Y
- real(kind=8),PARAMETER :: ZERO=0.0D+00, ONE=1.0D+00
- real(kind=8) :: ALPHA,BETA
- integer :: i,locy
-!
-! CLONE OF -DGEMV- WRITTEN BY MIKE SCHMIDT
-!
- LOCY = 1
- IF(FORMA.EQ.'T') GO TO 200
-!
-! Y = ALPHA * A * X + BETA * Y
-!
- IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN
- DO 110 I=1,M
- Y(LOCY) = DDOT(N,A(I,1),LDA,X,INCX)
- LOCY = LOCY+INCY
- 110 CONTINUE
- ELSE
- DO 120 I=1,M
- Y(LOCY) = ALPHA*DDOT(N,A(I,1),LDA,X,INCX) + BETA*Y(LOCY)
- LOCY = LOCY+INCY
- 120 CONTINUE
- END IF
- RETURN
-!
-! Y = ALPHA * A-TRANSPOSE * X + BETA * Y
-!
- 200 CONTINUE
- IF(ALPHA.EQ.ONE .AND. BETA.EQ.ZERO) THEN
- DO 210 I=1,N
- Y(LOCY) = DDOT(M,A(1,I),1,X,INCX)
- LOCY = LOCY+INCY
- 210 CONTINUE
- ELSE
- DO 220 I=1,N
- Y(LOCY) = ALPHA*DDOT(M,A(1,I),1,X,INCX) + BETA*Y(LOCY)
- LOCY = LOCY+INCY
- 220 CONTINUE
- END IF
- return
- end subroutine DGEMV
-!-----------------------------------------------------------------------------
-!-----------------------------------------------------------------------------
- end module MD_calc
projects / unres4.git / blobdiff