1c 2c 3c ################################################### 4c ## COPYRIGHT (C) 1996 by Jay William Ponder ## 5c ## All Rights Reserved ## 6c ################################################### 7c 8c ################################################################## 9c ## ## 10c ## subroutine deflate -- eigenvalues by method of deflation ## 11c ## ## 12c ################################################################## 13c 14c 15c "deflate" uses the power method with deflation to compute the 16c few largest eigenvalues and eigenvectors of a symmetric matrix 17c 18c n dimension of the matrix to be diagonalized 19c nv number of largest eigenvalues to be extracted 20c a input with the matrix to be diagonalized; only 21c the lower triangle and diagonal are required 22c ev returned with the eigenvalues in descending order 23c vec returned with the eigenvectors of the matrix 24c work local vector containing temporary work space 25c 26c 27 subroutine deflate (n,nv,a,ev,vec) 28 use iounit 29 implicit none 30 integer i,j,k,n,nv 31 integer iter,maxiter 32 real*8 random,eps 33 real*8 dot1,dot2,ratio 34 real*8 ev(*) 35 real*8, allocatable :: work(:) 36 real*8 a(n,*) 37 real*8 vec(n,*) 38 external random 39c 40c 41c initialize number of iterations and convergence criteria 42c 43 maxiter = 500 44 eps = 1.0d-6 45c 46c use identity vector as initial guess for eigenvectors 47c 48 do j = 1, nv 49 do i = 1, n 50 vec(i,j) = 1.0d0 51 end do 52 end do 53c 54c perform dynamic allocation of some local arrays 55c 56 allocate (work(n)) 57c 58c find the few largest eigenvalues and eigenvectors 59c 60 do k = 1, nv 61 ev(k) = 0.0d0 62 dot1 = 0.0d0 63 do i = 1, n 64 work(i) = 0.0d0 65 do j = 1, i-1 66 work(i) = work(i) + a(i,j)*vec(j,k) 67 end do 68 do j = i, n 69 work(i) = work(i) + a(j,i)*vec(j,k) 70 end do 71 dot1 = dot1 + work(i)**2 72 end do 73c 74c if in or near null space, use random guess as eigenvector 75c 76 if (dot1 .le. 100.0d0*eps*dble(n)) then 77 do i = 1, n 78 work(i) = random () 79 end do 80 end if 81c 82c find the current eigenvalue by iterating to convergence; 83c first multiply vector by matrix and compute dot products 84c 85 do iter = 1, maxiter 86 dot1 = 0.0d0 87 dot2 = 0.0d0 88 do i = 1, n 89 vec(i,k) = 0.0d0 90 do j = 1, i-1 91 vec(i,k) = vec(i,k) + a(i,j)*work(j) 92 end do 93 do j = i, n 94 vec(i,k) = vec(i,k) + a(j,i)*work(j) 95 end do 96 dot1 = dot1 + vec(i,k)**2 97 dot2 = dot2 + vec(i,k)*work(i) 98 end do 99c 100c normalize new eigenvector and substitute for old one 101c 102 ratio = abs((ev(k)-dot2) / dot2) 103 ev(k) = dot2 104 dot1 = sqrt(dot1) 105 do i = 1, n 106 vec(i,k) = vec(i,k) / dot1 107 work(i) = vec(i,k) 108 end do 109 if (ratio .lt. eps) goto 20 110 end do 111 write (iout,10) k 112 10 format (/,' DEFLATE -- Eigenvalue',i3,' not Fully Converged') 113c 114c eliminate the current eigenvalue from the matrix 115c 116 20 continue 117 do i = 1, n 118 do j = i, n 119 a(j,i) = a(j,i) - ev(k)*vec(i,k)*vec(j,k) 120 end do 121 end do 122 end do 123c 124c perform deallocation of some local arrays 125c 126 deallocate (work) 127 return 128 end 129