MODULE AS60 ! Calculate the eigenvalues and eigenvectors of a real symmetric matrix. IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(12, 60) REAL (dp), PARAMETER :: zero = 0.0_dp, one = 1.0_dp CONTAINS SUBROUTINE tdiag (n, tol, a, d, e, z, maxdim, ifault) ! Code converted using TO_F90 by Alan Miller ! Date: 2003-12-03 Time: 12:16:20 INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: tol INTEGER, INTENT(IN OUT) :: maxdim REAL (dp), INTENT(IN) :: a(maxdim,maxdim) REAL (dp), INTENT(OUT) :: d(maxdim) REAL (dp), INTENT(OUT) :: e(maxdim) REAL (dp), INTENT(OUT) :: z(maxdim,maxdim) INTEGER, INTENT(OUT) :: ifault REAL (dp), PARAMETER :: eta = 1.0D-37, eps = 1.0D-14 ! Algorithm as 60.1 appl.statist. (1973) vol.22 no.2 ! Reduces real symmetric matrix to tridiagonal form ! tol is a machine dependent constant , tol = eta/eps , where ! eta is the smallest positive number representable in the ! computer and eps is the smallest positive number for which ! 1+eps.ne.1. ! eta=eps*tol ! eps=0.7105427358e-14 ! tol=0.3131513063e-293 ! precis=1.0e-14 ! NB ! Real constants must be <= 15 decimal digits ! the range of a real constant is from 1.0e-293 to 1.0e+322 INTEGER :: i, i1, j, j1, k, l REAL (dp) :: f, g, h, hh tol=eta/eps ifault=1 IF (n <= 1) RETURN ifault=0 DO i=1,n z(i,1:i)=a(i,1:i) END DO i=n DO i1=2,n l=i-2 f=z(i,i-1) g=zero IF (l < 1) GO TO 30 g = SUM( z(i,1:l)**2 ) 30 h=g + f*f ! If g is too small for orthogonality to be guaranteed, the ! transformation is skipped IF (g > tol) GO TO 40 e(i)=f d(i)=zero GO TO 100 40 l=l+1 g=SQRT(h) IF (f >= zero) g=-g e(i)=g h=h - f*g z(i,i-1)=f-g f=zero DO j=1,l z(j,i)=z(i,j)/h g=zero ! Form element of a * u DO k=1,j g=g + z(j,k)*z(i,k) END DO IF (j >= l) GO TO 70 j1=j+1 DO k=j1,l g=g + z(k,j)*z(i,k) END DO ! Form element of p 70 e(j)=g/h f=f + g*z(j,i) END DO ! Form k hh=f/(h+h) ! Form reduced a DO j=1,l f=z(i,j) g=e(j)-hh*f e(j)=g DO k=1,j z(j,k)=z(j,k) - f*e(k) - g*z(i,k) END DO END DO d(i)=h 100 i=i-1 END DO d(1)=zero e(1)=zero ! Accumulation of transformation matrices DO i=1,n l=i-1 IF (d(i) == zero .OR. l == 0) GO TO 140 DO j=1,l g=zero DO k=1,l g=g + z(i,k)*z(k,j) END DO DO k=1,l z(k,j)=z(k,j) - g*z(k,i) END DO END DO 140 d(i)=z(i,i) z(i,i)=one IF (l == 0) CYCLE DO j=1,l z(i,j)=zero z(j,i)=zero END DO END DO RETURN END SUBROUTINE tdiag SUBROUTINE lrvt (n, precis, d, e, z, ifault, maxdim) INTEGER, INTENT(IN) :: n REAL (dp), INTENT(OUT) :: precis INTEGER, INTENT(IN) :: maxdim REAL (dp), INTENT(IN OUT) :: d(maxdim) REAL (dp), INTENT(OUT) :: e(maxdim) REAL (dp), INTENT(IN OUT) :: z(maxdim,maxdim) INTEGER, INTENT(OUT) :: ifault ! Algorithm AS 60.2 appl.statist. (1973) vol.22, no.2 ! Finds latent roots and vectors of tridiagonal matrix INTEGER, PARAMETER :: mits = 30 REAL (dp), PARAMETER :: two = 2.0D0 INTEGER :: i, i1, j, jj, k, l, m, m1, n1 REAL (dp) :: b, c, f, g, h, p, pr, r, s precis=1.0D-14 ifault=2 IF (n <= 1) RETURN ifault=1 n1=n-1 DO i=2,n e(i-1)=e(i) END DO e(n)=zero b=zero f=zero DO l=1,n jj=0 h=precis*(ABS(d(l))+ABS(e(l))) IF (b < h) b=h ! Look for small sub-diagonal element DO m1=l,n m=m1 IF (ABS(e(m)) <= b) EXIT END DO IF (m == l) GO TO 90 40 IF (jj == mits) RETURN jj=jj+1 ! Form shift p=(d(l+1)-d(l))/(two*e(l)) r=SQRT(p*p + one) pr=p+r IF (p < zero) pr=p-r h=d(l)-e(l)/pr DO i=l,n d(i)=d(i)-h END DO f=f+h ! QL transformation p=d(m) c=one s=zero m1=m-1 i=m DO i1=l,m1 j=i i=i-1 g=c*e(i) h=c*p IF (ABS(p) >= ABS(e(i))) GO TO 60 c=p/e(i) r=SQRT(c*c + one) e(j)=s*e(i)*r s=one/r c=c/r GO TO 70 60 c=e(i)/p r=SQRT(c*c + one) e(j)=s*p*r s=c/r c=one/r 70 p=c*d(i) - s*g d(j)=h + s*(c*g + s*d(i)) ! Form vector DO k=1,n h=z(k,j) z(k,j)=s*z(k,i) + c*h z(k,i)=c*z(k,i) - s*h END DO END DO e(l)=s*p d(l)=c*p IF (ABS(e(l)) > b) GO TO 40 90 d(l)=d(l) + f END DO ! Order latent roots and vectors DO i=1,n1 k=i p=d(i) i1=i+1 DO j=i1,n IF (d(j) <= p) CYCLE k=j p=d(j) END DO IF (k == i) CYCLE d(k)=d(i) d(i)=p DO j=1,n p=z(j,i) z(j,i)=z(j,k) z(j,k)=p END DO END DO ifault=0 RETURN END SUBROUTINE lrvt END MODULE AS60