1! 2! Copyright (C) 2001-2013 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8! 9 10! 11! 12!---------------------------------------------------------------------- 13subroutine cgsolve_all_gamma (h_psi, cg_psi, e, d0psi, dpsi, h_diag, & 14 ndmx, ndim, ethr, ik, kter, conv_root, anorm, nbnd, npol) 15 !---------------------------------------------------------------------- 16 ! 17 ! iterative solution of the linear system: 18 ! 19 ! ( h - e + Q ) * dpsi = d0psi (1) 20 ! 21 ! where h is a complex hermitean matrix, e is a real sca 22 ! dpsi and d0psi are complex vectors 23 ! 24 ! on input: 25 ! h_psi EXTERNAL name of a subroutine: 26 ! h_psi(ndim,psi,psip) 27 ! Calculates H*psi products. 28 ! Vectors psi and psip should be dimensined 29 ! (ndmx,nvec). nvec=1 is used! 30 ! 31 ! cg_psi EXTERNAL name of a subroutine: 32 ! g_psi(ndmx,ndim,notcnv,psi,e) 33 ! which calculates (h-e)^-1 * psi, with 34 ! some approximation, e.g. (diag(h)-e) 35 ! 36 ! e real unperturbed eigenvalue. 37 ! 38 ! dpsi contains an estimate of the solution 39 ! vector. 40 ! 41 ! d0psi contains the right hand side vector 42 ! of the system. 43 ! 44 ! ndmx integer row dimension of dpsi, ecc. 45 ! 46 ! ndim integer actual row dimension of dpsi 47 ! 48 ! ethr real convergence threshold. solution 49 ! improvement is stopped when the error in 50 ! eq (1), defined as l.h.s. - r.h.s., becomes 51 ! less than ethr in norm. 52 ! 53 ! on output: dpsi contains the refined estimate of the 54 ! solution vector. 55 ! 56 ! d0psi is corrupted on exit 57 ! 58 ! revised (extensively) 6 Apr 1997 by A. Dal Corso & F. Mauri 59 ! revised (to reduce memory) 29 May 2004 by S. de Gironcoli 60 ! 61 USE kinds, ONLY : DP 62 USE mp_pools, ONLY : intra_pool_comm 63 USE mp, ONLY : mp_sum 64 USE control_flags, ONLY : gamma_only 65 USE gvect, ONLY : gstart 66 67 implicit none 68 ! 69 ! first the I/O variables 70 ! 71 integer :: ndmx, & ! input: the maximum dimension of the vectors 72 ndim, & ! input: the actual dimension of the vectors 73 kter, & ! output: counter on iterations 74 nbnd, & ! input: the number of bands 75 npol, & ! input: number of components of the wavefunctions 76 ik ! input: the k point 77 78 real(DP) :: & 79 e(nbnd), & ! input: the actual eigenvalue 80 anorm, & ! output: the norm of the error in the solution 81 h_diag(ndmx*npol,nbnd), & ! input: an estimate of ( H - \epsilon ) 82 ethr ! input: the required precision 83 84 complex(DP) :: & 85 dpsi (ndmx*npol, nbnd), & ! output: the solution of the linear syst 86 d0psi (ndmx*npol, nbnd) ! input: the known term 87 88 logical :: conv_root ! output: if true the root is converged 89 external h_psi ! input: the routine computing h_psi 90 external cg_psi ! input: the routine computing cg_psi 91 ! 92 ! here the local variables 93 ! 94 integer, parameter :: maxter = 200 95 ! the maximum number of iterations 96 integer :: iter, ibnd, lbnd 97 ! counters on iteration, bands 98 integer , allocatable :: conv (:) 99 ! if 1 the root is converged 100 101 complex(DP), allocatable :: g (:,:), t (:,:), h (:,:), hold (:,:) 102 ! the gradient of psi 103 ! the preconditioned gradient 104 ! the delta gradient 105 ! the conjugate gradient 106 ! work space 107 complex(DP) :: dcgamma, dclambda 108 ! the ratio between rho 109 ! step length 110 complex(DP), external :: zdotc 111 REAL(kind=dp), EXTERNAL :: ddot 112 ! the scalar product 113 real(DP), allocatable :: rho (:), rhoold (:), eu (:), a(:), c(:) 114 ! the residue 115 ! auxiliary for h_diag 116 real(DP) :: kter_eff 117 ! account the number of iterations with b 118 ! coefficient of quadratic form 119 ! 120 call start_clock ('cgsolve') 121 allocate ( g(ndmx*npol,nbnd), t(ndmx*npol,nbnd), h(ndmx*npol,nbnd), & 122 hold(ndmx*npol ,nbnd) ) 123 allocate (a(nbnd), c(nbnd)) 124 allocate (conv ( nbnd)) 125 allocate (rho(nbnd),rhoold(nbnd)) 126 allocate (eu ( nbnd)) 127 ! WRITE( stdout,*) g,t,h,hold 128 129 kter_eff = 0.d0 130 do ibnd = 1, nbnd 131 conv (ibnd) = 0 132 enddo 133 g=(0.d0,0.d0) 134 t=(0.d0,0.d0) 135 h=(0.d0,0.d0) 136 hold=(0.d0,0.d0) 137 do iter = 1, maxter 138 ! 139 ! compute the gradient. can reuse information from previous step 140 ! 141 if (iter == 1) then 142 call h_psi (ndim, dpsi, g, e, ik, nbnd) 143 do ibnd = 1, nbnd 144 call zaxpy (ndim, (-1.d0,0.d0), d0psi(1,ibnd), 1, g(1,ibnd), 1) 145 enddo 146 IF (npol==2) THEN 147 do ibnd = 1, nbnd 148 call zaxpy (ndim, (-1.d0,0.d0), d0psi(ndmx+1,ibnd), 1, & 149 g(ndmx+1,ibnd), 1) 150 enddo 151 END IF 152 endif 153 ! 154 ! compute preconditioned residual vector and convergence check 155 ! 156 lbnd = 0 157 do ibnd = 1, nbnd 158 if (conv (ibnd) .eq.0) then 159 lbnd = lbnd+1 160 call zcopy (ndmx*npol, g (1, ibnd), 1, h (1, ibnd), 1) 161 call cg_psi(ndmx, ndim, 1, h(1,ibnd), h_diag(1,ibnd) ) 162 163 IF (gamma_only) THEN 164 rho(lbnd)=2.0d0*ddot(2*ndmx*npol,h(1,ibnd),1,g(1,ibnd),1) 165 IF(gstart==2) THEN 166 rho(lbnd)=rho(lbnd)-DBLE(h(1,ibnd))*DBLE(g(1,ibnd)) 167 ENDIF 168 ELSE 169 rho(lbnd) = zdotc (ndmx*npol, h(1,ibnd), 1, g(1,ibnd), 1) 170 ENDIF 171 172 endif 173 enddo 174 kter_eff = kter_eff + DBLE (lbnd) / DBLE (nbnd) 175#if defined(__MPI) 176 call mp_sum( rho(1:lbnd) , intra_pool_comm ) 177#endif 178 do ibnd = nbnd, 1, -1 179 if (conv(ibnd).eq.0) then 180 rho(ibnd)=rho(lbnd) 181 lbnd = lbnd -1 182 anorm = sqrt (rho (ibnd) ) 183! write(6,*) ibnd, anorm 184 if (anorm.lt.ethr) conv (ibnd) = 1 185 endif 186 enddo 187! 188 conv_root = .true. 189 do ibnd = 1, nbnd 190 conv_root = conv_root.and. (conv (ibnd) .eq.1) 191 enddo 192 if (conv_root) goto 100 193 ! 194 ! compute the step direction h. Conjugate it to previous step 195 ! 196 lbnd = 0 197 do ibnd = 1, nbnd 198 if (conv (ibnd) .eq.0) then 199! 200! change sign to h 201! 202 call dscal (2 * ndmx * npol, - 1.d0, h (1, ibnd), 1) 203 if (iter.ne.1) then 204 dcgamma = rho (ibnd) / rhoold (ibnd) 205 call zaxpy (ndmx*npol, dcgamma, hold (1, ibnd), 1, h (1, ibnd), 1) 206 endif 207 208! 209! here hold is used as auxiliary vector in order to efficiently compute t = A*h 210! it is later set to the current (becoming old) value of h 211! 212 lbnd = lbnd+1 213 call zcopy (ndmx*npol, h (1, ibnd), 1, hold (1, lbnd), 1) 214 eu (lbnd) = e (ibnd) 215 endif 216 enddo 217 ! 218 ! compute t = A*h 219 ! 220 call h_psi (ndim, hold, t, eu, ik, lbnd) 221 ! 222 ! compute the coefficients a and c for the line minimization 223 ! compute step length lambda 224 lbnd=0 225 do ibnd = 1, nbnd 226 if (conv (ibnd) .eq.0) then 227 lbnd=lbnd+1 228 IF (gamma_only) THEN 229 a(lbnd) = 2.0d0*ddot(2*ndmx*npol,h(1,ibnd),1,g(1,ibnd),1) 230 c(lbnd) = 2.0d0*ddot(2*ndmx*npol,h(1,ibnd),1,t(1,lbnd),1) 231 IF (gstart == 2) THEN 232 a(lbnd)=a(lbnd)-DBLE(h(1,ibnd))*DBLE(g(1,ibnd)) 233 c(lbnd)=c(lbnd)-DBLE(h(1,ibnd))*DBLE(t(1,lbnd)) 234 ENDIF 235 ELSE 236 a(lbnd) = zdotc (ndmx*npol, h(1,ibnd), 1, g(1,ibnd), 1) 237 c(lbnd) = zdotc (ndmx*npol, h(1,ibnd), 1, t(1,lbnd), 1) 238 ENDIF 239 end if 240 end do 241#if defined(__MPI) 242 call mp_sum( a(1:lbnd), intra_pool_comm ) 243 call mp_sum( c(1:lbnd), intra_pool_comm ) 244#endif 245 lbnd=0 246 do ibnd = 1, nbnd 247 if (conv (ibnd) .eq.0) then 248 lbnd=lbnd+1 249 dclambda = CMPLX( - a(lbnd) / c(lbnd), 0.d0,kind=DP) 250 ! 251 ! move to new position 252 ! 253 call zaxpy (ndmx*npol, dclambda, h(1,ibnd), 1, dpsi(1,ibnd), 1) 254 ! 255 ! update to get the gradient 256 ! 257 !g=g+lam 258 call zaxpy (ndmx*npol, dclambda, t(1,lbnd), 1, g(1,ibnd), 1) 259 ! 260 ! save current (now old) h and rho for later use 261 ! 262 call zcopy (ndmx*npol, h(1,ibnd), 1, hold(1,ibnd), 1) 263 rhoold (ibnd) = rho (ibnd) 264 endif 265 enddo 266 enddo 267100 continue 268 kter = kter_eff 269 deallocate (eu) 270 deallocate (rho, rhoold) 271 deallocate (conv) 272 deallocate (a,c) 273 deallocate (g, t, h, hold) 274 275 call stop_clock ('cgsolve') 276 return 277end subroutine cgsolve_all_gamma 278