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