1! 2! Copyright (C) 2001-2004 PWSCF 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!---------------------------------------------------------------------- 10subroutine gmressolve_all (h_psi, cg_psi, e, d0psi, dpsi, h_diag, & 11 ndmx, ndim, ethr, ik, kter, conv_root, anorm, nbnd, m) 12 !---------------------------------------------------------------------- 13 ! 14 ! iterative solution of the linear system by GMRES(m) method: 15 ! 16 ! ( h - e + Q ) * dpsi = d0psi (1) 17 ! 18 ! where h is a complex hermitean matrix, e is a complex sca 19 ! dpsi and d0psi are complex vectors 20 ! 21 ! on input: 22 ! h_psi EXTERNAL name of a subroutine: 23 ! h_psi(ndim,psi,psip) 24 ! Calculates H*psi products. 25 ! Vectors psi and psip should be dimensined 26 ! (ndmx,nvec). nvec=1 is used! 27 ! 28 ! cg_psi EXTERNAL name of a subroutine: 29 ! g_psi(ndmx,ndim,notcnv,psi,e) 30 ! which calculates (h-e)^-1 * psi, with 31 ! some approximation, e.g. (diag(h)-e) 32 ! 33 ! e complex unperturbed eigenvalue plus 34 ! imaginary frequency. 35 ! 36 ! dpsi contains an estimate of the solution 37 ! vector. 38 ! 39 ! d0psi contains the right hand side vector 40 ! of the system. 41 ! 42 ! ndmx integer row dimension of dpsi, ecc. 43 ! 44 ! ndim integer actual row dimension of dpsi 45 ! 46 ! ethr real convergence threshold. solution 47 ! improvement is stopped when the error in 48 ! eq (1), defined as l.h.s. - r.h.s., becomes 49 ! less than ethr in norm. 50 ! 51 ! m integer # of basis vectors 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_bands, ONLY: intra_bgrp_comm 63 USE mp, ONLY: mp_sum 64 65 implicit none 66 ! 67 ! first the I/O variables 68 ! 69 integer :: ndmx, & ! input: the maximum dimension of the vectors 70 ndim, & ! input: the actual dimension of the vectors 71 kter, & ! output: counter on iterations 72 nbnd, & ! input: the number of bands 73 ik, & ! input: the k point 74 m ! # of basic vector 75 76 real(kind=DP) :: & 77 anorm, & ! output: the norm of the error in the solution 78 ethr ! input: the required precision 79 complex(kind=DP) :: h_diag(ndmx,nbnd) ! input: an estimate of ( H - \epsilon - iu ) 80 81 complex(kind=DP) :: & 82 e(nbnd), & ! input: the actual eigenvalue plus imaginary freq. 83 dpsi (ndmx, nbnd), & ! output: the solution of the linear syst 84 d0psi (ndmx, nbnd) ! input: the known term 85 86 logical :: conv_root ! output: if true the root is converged 87 external h_psi ! input: the routine computing h_psi 88 external cg_psi ! input: the routine computing cg_psi 89 ! 90 ! here the local variables 91 ! 92 integer, parameter :: maxter = 5000 93 ! the maximum number of iterations 94 integer :: iter, ibnd, i, j, bnd 95 ! counters on iteration, bands 96 ! control variables 97 integer , allocatable :: conv (:) 98 ! if 1 the root is converged 99 100 complex(kind=DP), allocatable :: r (:,:), v(:,:,:), w (:,:)!, zz(:,:), p(:,:), pp(:,:) 101 ! the gradient of psi 102 ! the preconditioned gradient 103 ! the delta gradient 104 ! the conjugate gradient 105 ! work space 106 complex(kind=DP) :: bk, ak 107 ! the ratio between rho 108 ! step length 109 complex(kind=DP), external :: zdotc 110 ! the scalar product 111 real(kind=DP) :: t 112 complex(kind=DP):: c, s, ei 113 real(kind=DP), allocatable :: bet (:) 114 real(kind=DP), allocatable :: res (:) 115 complex(kind=DP) :: hm (m+1,m), & ! the Hessenberg matrix 116 e1(m+1) ! unit vector 117 complex(kind=DP) :: hm4para(1) ! temp variable for hm in paralell calculation 118! real(kind=DP), allocatable :: rho (:), rhoold (:), eu (:), a(:), c(:) 119 ! the residue 120 ! auxiliary for h_diag 121 real(kind=DP) :: kter_eff 122 ! account the number of iterations with b 123 ! coefficient of quadratic form 124 ! 125 integer :: lbnd 126 ! 127 ! 128 ! 129 call start_clock ('gmres_solve') 130 ! 131 if (m .lt. 1) then 132 write(*,*) '# of basis vectors is less than 1. Stop' 133 stop 134 else if (m .gt. 30) then 135 write(*,*) '# of basis vectors is too large. Stop' 136 stop 137 endif 138 ! 139 allocate ( r(ndmx,nbnd), v(ndmx,nbnd,m+1), w(ndmx,nbnd)) 140 allocate (conv ( nbnd)) 141 allocate (bet(nbnd), res(nbnd)) 142 ! WRITE( stdout,*) g,t,h,hold 143 144 kter_eff = 0.d0 145 do ibnd = 1, nbnd 146 conv (ibnd) = 0 147 enddo 148 ! 149 do iter = 1, maxter 150 ! 151!print*, 'iter=', iter 152 do ibnd = 1, nbnd ! loop over bands 153 ! 154 if (conv(ibnd) .eq. 0) then 155 ! 156 ! preliminary step to construct the basis set 157 ! 158 ! r = H*dpsi 159 call h_psi (ndim, dpsi(1,ibnd), r(1,ibnd), e(ibnd), ik, 1) 160!print*,'dpsi',sum(dpsi),sum(d0psi) 161 ! 162 ! r = H*dpsi - d0psi 163 call zaxpy (ndim, (-1.d0,0.d0), d0psi(1,ibnd), 1, r(1,ibnd), 1) 164!print*,'r1',sum(dpsi),sum(d0psi) 165 ! change the size of r : r = d0psi - H*dpsi 166 call dscal (2 * ndim, - 1.d0, r (1, ibnd), 1) 167!print*,'r2',sum(dpsi),sum(d0psi) 168 ! compute the preconditioned r : r = M^-1*r 169 call cg_psi(ndmx, ndim, 1, r(1,ibnd), h_diag(1,ibnd), 1 ) 170!print*,'r3',sum(dpsi),sum(d0psi) 171 ! norm of pre. r : bet = |r| 172 bet(ibnd) = zdotc (ndim, r(1,ibnd), 1, r(1,ibnd), 1) 173#if defined(__MPI) 174 call mp_sum ( bet(ibnd), intra_bgrp_comm ) 175#endif 176 bet(ibnd) = sqrt( bet(ibnd) ) 177 ! 178 endif 179 ! 180 enddo 181 ! 182 ! check the convergence 183 ! 184 lbnd = 0 185 do ibnd = 1, nbnd 186 ! 187 if ( conv(ibnd) .eq. 0 ) then 188 lbnd = lbnd + 1 189!if (mod(iter,10) .eq. 0) print*, iter, bet(ibnd), ethr 190 if (bet(ibnd) .lt. ethr) conv(ibnd) = 1 191 endif 192 ! 193 enddo 194 kter_eff = kter_eff + DBLE (lbnd) / DBLE (nbnd) 195 ! 196 conv_root = .true. 197 do ibnd = 1, nbnd 198 conv_root = conv_root .and. (conv (ibnd) .eq. 1) 199 enddo 200 if (conv_root) goto 100 201 ! 202 ! 203 ! 204 do ibnd = 1, nbnd 205 ! 206 if ( conv(ibnd) .eq. 0 ) then 207 ! 208 hm (:,:) = (0.d0, 0.d0) 209 ! normalize pre. r and keep in v(1) 210 call dscal (2 * ndim, 1.d0/bet(ibnd), r (1, ibnd), 1) 211 j = 1 212 call zcopy (ndim, r (1, ibnd), 1, v (1, ibnd, j), 1) 213!print*,'v',sum(r(1:ndim,ibnd)) 214 ! 215 ! 216 ! loop to construct basis set 217 ! 218 ! 219 do j = 1, m 220 ! w = A*v 221 call h_psi (ndim, v(1,ibnd,j), w(1,ibnd), e, ik, 1) ! NEED to be checked 222!print*,'w1',sum(w(:,ibnd)) 223 ! 224 ! compute w = M^-1*A*v 225 call cg_psi(ndmx, ndim, 1, w(1,ibnd), h_diag(1,ibnd), 1 ) 226!print*,'w2',sum(w(:,ibnd)) 227!print*,'h_diag',sum(h_diag) 228 ! 229 do i = 1, j 230 ! 231 ! compute hm(i,j) 232! hm(i,j) = zdotc (ndim, w(1,ibnd), 1, v(1,ibnd,i), 1) 233 hm4para(1) = zdotc (ndim, w(1,ibnd), 1, v(1,ibnd,i), 1) 234#if defined(__MPI) 235 call mp_sum ( hm4para, intra_bgrp_comm ) 236#endif 237 hm(i,j) = hm4para(1) 238 ! w = w - hm_ij*v_i 239 call zaxpy (ndim, -hm(i,j), v(1,ibnd,i), 1, w(1,ibnd), 1) 240 ! 241 enddo 242 ! compute hm(j+1,j) 243! hm(j+1,j) = zdotc (ndim, w(1,ibnd), 1, w(1,ibnd), 1) 244 hm4para(1) = zdotc (ndim, w(1,ibnd), 1, w(1,ibnd), 1) 245#if defined(__MPI) 246 call mp_sum ( hm4para, intra_bgrp_comm ) 247#endif 248 hm(j+1,j) = hm4para(1) 249 ! compute v(j+1) 250 call dscal (2 * ndim, 1.d0/real(hm(j+1,j)), w (1, ibnd), 1) 251 call zcopy (ndim, w (1, ibnd), 1, v (1, ibnd, j+1), 1) 252 ! 253 enddo 254 ! 255 ! compute ym that minimize |beta*e_1 -hm*y| 256 ! 257 ! initilize vector e1 258 e1(1) = 1.d0 * bet(ibnd) 259 e1(2:m+1) = 0.d0 260 ! 261 ! transform hm to upper triangle matrix 262 do i = 1, m 263 ! 264 t = sqrt( abs(hm(i,i))**2 + abs(hm(i+1,i))**2 ) 265 c = hm(i,i) / t 266 s = hm(i+1,i) / t 267 ! 268 do j = i, m 269 ! 270 ei = hm(i,j) 271 hm(i,j) = hm(i,j) * c + hm(i+1,j) * s 272 hm(i+1,j) = - s * ei + c * hm(i+1,j) 273 enddo 274 ! 275 ei = e1(i) 276 e1(i) = e1(i)*c + e1(i+1)*s 277 e1(i+1) = - ei*s + e1(i+1)*c 278 ! 279 enddo 280 ! 281 res(ibnd) = e1(m+1) 282 ! 283 ! back subtitution to find ym (kept in e1) 284 e1(m+1) = (0.d0, 0.d0) 285 e1(m) = e1(m) / hm(m,m) 286 ! 287 do i = m-1, 1, -1 288 do j = m, i+1, -1 289 e1(i) = e1(i) - e1(j)*hm(i,j) 290 enddo 291 e1(i) = e1(i) / hm(i,i) 292 enddo 293 ! 294 ! compute the new dpsi 295 do i = 1, m 296 do j = 1, ndmx 297 dpsi(j, ibnd) = dpsi(j, ibnd) + e1(i)*v(j,ibnd,i) 298 enddo 299 enddo 300 ! 301 end if 302 ! 303 enddo ! of loop over bands 304 ! 305 enddo ! loop over iteration 306 ! 307100 continue 308 kter = kter_eff 309 ! 310 deallocate (bet, res) 311 deallocate (conv) 312 deallocate (r, v, w) 313 ! 314 call stop_clock ('gmres_solve') 315 ! 316 return 317 ! 318end subroutine gmressolve_all 319