Branch data Line data Source code
1 : : !!****m* ABINIT/m_gwls_lineqsolver
2 : : !! NAME
3 : : !! m_gwls_lineqsolver
4 : : !!
5 : : !! FUNCTION
6 : : !! .
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2009-2026 ABINIT group (JLJ, BR, MC)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : :
23 : :
24 : : !---------------------------------------------------------------------
25 : : ! Module to solve A.x = b efficiently, where A will involve
26 : : ! the Hamiltonian.
27 : : !---------------------------------------------------------------------
28 : :
29 : :
30 : : module m_gwls_lineqsolver
31 : : !----------------------------------------------------------------------------------------------------
32 : : ! This module contains routines to solve A x = b interatively to solve the Sternheimer equation
33 : : ! in various contexts...
34 : : !----------------------------------------------------------------------------------------------------
35 : :
36 : :
37 : : ! local modules
38 : : use m_gwls_utility
39 : : use m_gwls_wf
40 : : use m_gwls_hamiltonian
41 : :
42 : : ! abinit modules
43 : : use defs_basis
44 : : use m_abicore
45 : : use m_xmpi
46 : : use m_bandfft_kpt
47 : : use m_cgtools
48 : :
49 : : use m_time, only : timab
50 : : use m_io_tools, only : get_unit
51 : :
52 : : implicit none
53 : : save
54 : : private
55 : : !!***
56 : :
57 : : logical :: activate_inf_shift_poles = .false.
58 : : real(dp) :: inf_shift_poles = 1.0d-4
59 : : !!***
60 : :
61 : : public :: sqmr, qmr, activate_inf_shift_poles, inf_shift_poles
62 : : !!***
63 : :
64 : : contains
65 : :
66 : : !!****f* m_hamiltonian/sqmr
67 : : !! NAME
68 : : !! sqmr
69 : : !!
70 : : !! FUNCTION
71 : : !! .
72 : : !!
73 : : !! INPUTS
74 : : !!
75 : : !! OUTPUT
76 : : !!
77 : : !! SOURCE
78 : :
79 : 1840 : subroutine sqmr(b,x,lambda,project_on_what,omega,omega_imaginary,kill_Pc_x)
80 : : !--------------------------------------------------------------------------------
81 : : ! This subroutine solves the linear algebra problem
82 : : !
83 : : ! A x = b
84 : : !
85 : : ! Where:
86 : : ! INPUT
87 : : ! -----
88 : : ! real(dp) b right-hand-side of the equation to be solved
89 : : ! real(dp) omega *OPTIONAL* frequency used in building A
90 : : ! logical omega_imaginary *OPTIONAL* is the frequency imaginary?
91 : : ! real(dp) lambda value to be subtracted from the Hamiltonian
92 : : ! integer project_on_what flag which determines the projection scheme.
93 : : !
94 : : ! OUTPUT
95 : : ! -----
96 : : ! real(dp) x solution
97 : : !
98 : : ! Note that blocksize corresponds to the number of band processors; it is a global
99 : : ! variable defined in gwls_hamiltonian. The name is inspired from lobpcgwf.F90.
100 : : !
101 : : ! with:
102 : : ! omega omega_imaginary Operator
103 : : ! ------------------------------------------------------
104 : : ! absent - A = (H - lambda)
105 : : ! present - A = (H - lambda)^2 - omega^2
106 : : ! present present, true A = (H - lambda)^2 + omega^2
107 : : !
108 : : ! project_on_what action
109 : : ! ------------------------------------------------------
110 : : ! 0 no projection
111 : : ! 1 projection on conduction states
112 : : ! 2 projection out of subspace degenerate with lambda
113 : : ! 3 projection on states beyond all the states explicitly stored
114 : : !
115 : : ! NOTE: It is the developper's responsibility to apply (H-ev) on the input
116 : : ! if the frequency is not zero.
117 : : !--------------------------------------------------------------------------------
118 : :
119 : : !External variables
120 : : real(dp), intent(in) :: b(2,npw_g)
121 : : real(dp), intent(in) :: lambda
122 : : real(dp), intent(out) :: x(2,npw_g)
123 : : integer, intent(in) :: project_on_what
124 : : real(dp), intent(in), optional :: omega
125 : : logical, optional :: omega_imaginary, kill_Pc_x
126 : :
127 : : !Local variables
128 : : real(dp) :: norm, tmp(2), residual
129 : 1840 : real(dp), allocatable :: g(:), theta(:), rho(:), sigma(:), c(:)
130 : 1840 : real(dp), allocatable :: t(:,:), delta(:,:), r(:,:), d(:,:), w(:,:), wmb(:,:)
131 : : integer :: ii,ipw, k, l
132 : : real(dp):: signe
133 : : real(dp):: norm_Axb
134 : :
135 : :
136 : : real(dp):: norm_b, tol14
137 : :
138 : : integer :: min_index
139 : : logical :: singular
140 : : logical :: precondition_on
141 : : logical :: has_omega
142 : :
143 : : integer :: pow
144 : :
145 : : logical :: imaginary
146 : :
147 : : integer,save :: counter = 0
148 : : integer :: io_unit
149 : : character(128) :: filename
150 : : logical :: file_exists
151 : : logical :: head_node
152 : :
153 : : integer :: ierr
154 : :
155 : : integer :: mpi_communicator, mpi_rank, mpi_group
156 : :
157 : :
158 : :
159 : : ! timing
160 : : real(dp) :: tsec(2)
161 : : integer :: GWLS_TIMAB, OPTION_TIMAB
162 : :
163 : : ! *************************************************************************
164 : :
165 : : ! The processors communicate over FFT!
166 : 1840 : mpi_communicator = mpi_enreg%comm_fft
167 : :
168 : : ! what is the rank of this processor, within its group?
169 : 1840 : mpi_rank = mpi_enreg%me_fft
170 : :
171 : : ! Which group does this processor belong to, given the communicator?
172 : 1840 : mpi_group = mpi_enreg%me_band
173 : :
174 : : ! Do we have omega?
175 : 1840 : has_omega=present(omega)
176 : :
177 : : ! Test if the input has finite norm
178 : 1840 : tol14 = 1.0D-14
179 : 1840 : tmp = cg_zdotc(npw_g,b,b)
180 : 1840 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
181 : 1840 : norm_b = tmp(1)
182 : :
183 [ + + ]: 1840 : if (norm_b < tol14) then
184 : : ! Because of band parallelism, it is possible that sqmr gets a zero norm argument.
185 : : ! A | x> = 0 implies |x > = 0.
186 [ + + + + ]: 667248 : x(:,:) = zero
187 : : return
188 : : end if
189 : :
190 : 736 : GWLS_TIMAB = 1523
191 : 736 : OPTION_TIMAB = 1
192 : 736 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
193 : :
194 : : ! only the head node should write to the log file
195 : 736 : head_node = ( mpi_rank == 0 )
196 : :
197 : : !Memory allocation for local variables
198 [ - + + - : 2208 : ABI_MALLOC(g, (nline))
+ - + - -
+ - + ]
199 [ - + + - : 1472 : ABI_MALLOC(theta,(nline))
- + - + ]
200 [ - + + - : 1472 : ABI_MALLOC(rho, (nline))
- + - + ]
201 [ - + + - : 1472 : ABI_MALLOC(sigma,(nline))
- + - + ]
202 [ - + + - : 1472 : ABI_MALLOC(c, (nline))
- + - + ]
203 : :
204 [ - + + - : 2208 : ABI_MALLOC(t, (2,npw_g))
+ - + - -
+ - + ]
205 [ - + + - : 1472 : ABI_MALLOC(delta,(2,npw_g))
- + - + ]
206 [ - + + - : 1472 : ABI_MALLOC(r, (2,npw_g))
- + - + ]
207 [ - + + - : 1472 : ABI_MALLOC(d, (2,npw_g))
- + - + ]
208 [ - + + - : 1472 : ABI_MALLOC(w, (2,npw_g))
- + - + ]
209 [ - + + - : 1472 : ABI_MALLOC(wmb, (2,npw_g))
- + - + ]
210 : :
211 : :
212 : :
213 : : !Some vectors won't be filled (first iteration missing) so it's useful to initialise them.
214 [ + + ]: 74336 : g = zero
215 [ + + ]: 74336 : theta = zero
216 [ + + ]: 74336 : rho = zero
217 [ + + ]: 74336 : sigma = zero
218 [ + + ]: 74336 : c = zero
219 [ + + + + ]: 444832 : x = zero
220 [ + + + + ]: 444832 : delta = zero
221 [ + + + + ]: 444832 : r = zero
222 [ + + + + ]: 444832 : d = zero
223 [ + + + + ]: 444832 : w = zero
224 : :
225 : :
226 : : ! Determine if the frequency is imaginary
227 [ - + ]: 736 : if (has_omega .and. present(omega_imaginary)) then
228 : 0 : imaginary = omega_imaginary
229 : : else
230 : : imaginary = .false.
231 : : end if
232 : :
233 : :
234 : : ! Define the sign in front of (H-eig(v))**2.
235 : : ! If omega_imaginary is not given, we assume that omega is real (and sign=-1).
236 [ - + ]: 736 : if (has_omega) then
237 [ # # ]: 0 : if ( imaginary ) then
238 : : signe = one
239 : : else
240 : 0 : signe =-one
241 : : end if
242 : : end if
243 : :
244 : :
245 : : !Check for singularity problems
246 : : if (has_omega) then
247 [ # # # # : 0 : norm = minval(abs((eig(1:nbandv)-lambda)**2 + signe*(omega)**2))
# # # # ]
248 [ # # # # : 0 : min_index = minloc(abs((eig(1:nbandv)-lambda)**2 + signe*(omega)**2),1)
# # # # #
# ]
249 : : else
250 [ + - - + : 4416 : norm = minval(abs(eig(1:nbandv)-lambda))
- - + + ]
251 [ + - - + : 3680 : min_index = minloc(abs(eig(1:nbandv)-lambda),1)
- - + + +
+ ]
252 : : end if
253 : 736 : singular = norm < 1.0d-12
254 : :
255 : : !--------------------------------------------------------------------------------
256 : : ! If the linear operator has a kernel, then the intermediate vectors obtained in
257 : : ! SQMR must be projected out of this subspace several time at each iterations,
258 : : ! otherwise SQMR is unstable.
259 : : !
260 : : ! This is true even if the seed vector has been initially projected out of this
261 : : ! subspace, since the preconditionning will re-introduce a non-zero component in
262 : : ! the subspace of the kernel of the linear operator.
263 : : ! ===> Use project_on_what==2 in such cases.
264 : : !
265 : : ! Here, test if the operator is singular and if we are NOT projecting out of
266 : : ! the kernel.
267 : : ! ===> If true, stop the code.
268 : : !
269 : : !--------------------------------------------------------------------------------
270 : :
271 : : ! Quit if the operator has an uncontrolled kernel, a sign that the routine is being
272 : : ! misused by a developper...
273 [ + + + - : 736 : if (singular .and. ( (project_on_what==1 .and. (min_index > nbandv)) .or. project_on_what==0 )) then
+ - - + ]
274 : 0 : write(std_out,*) "ERROR - SQMR: Quasi-singuar problem treated, min. eigenvalue of A is ", norm," < 1d-12."
275 : 0 : write(std_out,*) " Yet, there is no projection out of the kernel of A. "
276 : :
277 [ # # # # ]: 0 : if (project_on_what==1 .and. (min_index > nbandv)) then
278 : 0 : write(std_out,*) " "
279 : 0 : write(std_out,*) " There is a projection on the conduction states, but A is singular in this "
280 : 0 : write(std_out,*) " subspace (the kernel contains state i=",min_index," > ",nbandv,"=# of valence states)."
281 : : end if
282 : :
283 : 0 : write(std_out,*) " "
284 : 0 : write(std_out,*) " In this situation, SQMR will be unstable. Use project_on_what==2 as an "
285 : 0 : write(std_out,*) " input argument of SQMR."
286 : 0 : write(std_out,*) " "
287 : 0 : write(std_out,*) " Decision taken to exit..."
288 : 0 : stop
289 : : end if
290 : :
291 : : !--------------------------------------------------------------------------------
292 : : ! Open a log file for the output of SQMR; only write if head of group!
293 : : !--------------------------------------------------------------------------------
294 [ + + ]: 736 : if (head_node) then
295 : :
296 : 576 : io_unit = get_unit()
297 : :
298 : 576 : write(filename,'(A,I0.4,A)') "SQMR_GROUP=",mpi_group,".log"
299 : :
300 : 576 : inquire(file=filename,exist=file_exists)
301 : :
302 [ + + ]: 576 : if (file_exists) then
303 : 573 : open( io_unit,file=filename,position='append',status=files_status_old)
304 : : else
305 : 3 : open( io_unit,file=filename,status=files_status_new)
306 : 3 : write(io_unit,10) "#======================================================================================="
307 : 3 : write(io_unit,10) "# "
308 : 3 : write(io_unit,10) "# This file contains information regarding the application of the SQMR scheme, "
309 : 3 : write(io_unit,10) "# for this MPI group. "
310 : 3 : write(io_unit,10) "#======================================================================================="
311 : 3 : flush(io_unit)
312 : : end if
313 : :
314 : 576 : counter = counter + 1
315 : 576 : write(io_unit,10) "# "
316 : 576 : write(io_unit,11) "# Call # ", counter
317 : 576 : write(io_unit,12) "# lambda = ",lambda," Ha "
318 [ - + ]: 576 : if (has_omega) then
319 : 0 : write(io_unit,12) "# omega = ",omega," Ha "
320 [ # # ]: 0 : if (imaginary) then
321 : 0 : write(io_unit,10) "# omega is imaginary "
322 : : else
323 : 0 : write(io_unit,10) "# omega is real "
324 : : end if
325 : : else
326 : 576 : write(io_unit,10) "# omega is absent "
327 : : end if
328 : :
329 : 576 : write(io_unit,13) "# project_on_what = ",project_on_what," "
330 : 576 : write(io_unit,13) "# "
331 [ - + ]: 576 : if (has_omega ) then
332 [ # # ]: 0 : if (imaginary) then
333 : 0 : write(io_unit,10) "# SOLVE ((H-lambda)^2 + omega^2) x = b"
334 : : else
335 : 0 : write(io_unit,10) "# SOLVE ((H-lambda)^2 - omega^2) x = b"
336 : : end if
337 : : else
338 : 576 : write(io_unit,10) "# SOLVE (H-lambda) x = b"
339 : : end if
340 : :
341 : 576 : flush(io_unit)
342 : : end if ! head_node
343 : : !--------------------------------------------------------------------------------
344 : : ! Precondition to accelerate convergence
345 : : !--------------------------------------------------------------------------------
346 : 736 : precondition_on = .true.
347 [ - + ]: 736 : if(imaginary) then
348 [ # # ]: 0 : if(omega > 10.0_dp) then
349 : : precondition_on = .false.
350 : : end if
351 : : end if
352 : :
353 : : ! DEBUG
354 : 736 : pow = project_on_what
355 : : !precondition_on = .false.
356 : :
357 : : !Prepare to precondition
358 : : if (precondition_on) then
359 [ - + ]: 736 : if ( imaginary ) then
360 : 0 : call set_precondition(lambda,omega)
361 : : else
362 : 736 : call set_precondition()
363 : : end if
364 : : else
365 : 0 : call unset_precondition()
366 : : end if
367 : :
368 : : !--------------------------------------------------------------------------------
369 : : !Initialisation
370 : : !--------------------------------------------------------------------------------
371 : :
372 : 736 : k = 1
373 : 736 : l = 1
374 : :
375 [ + + ]: 148768 : do ipw=1,npw_g
376 [ + + ]: 444832 : do ii=1,2
377 : 444096 : r(ii,ipw) = b(ii,ipw)
378 : : end do
379 : : end do
380 : :
381 [ + + ]: 736 : if (head_node) then
382 : 576 : write(io_unit,10) "# "
383 : 576 : write(io_unit,10) "# iteration approximate residual"
384 : 576 : write(io_unit,10) "#----------------------------------------"
385 : 576 : flush(io_unit)
386 : : end if
387 : :
388 : : do ! outer loop
389 : 736 : call precondition(d,r)
390 : :
391 : : ! g(k) = norm_k(r)
392 : 736 : tmp = cg_zdotc(npw_g,r,r)
393 : 736 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
394 : 736 : g(k) = dsqrt(tmp(1))
395 : :
396 : :
397 : : !tmp = scprod_k(r,d)
398 : 736 : tmp = cg_zdotc(npw_g,r,d)
399 : 736 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
400 : 736 : rho(k) = tmp(1)
401 : :
402 [ + + ]: 736 : if (head_node) then
403 : 576 : write(io_unit,16) k, g(k)**2
404 : 576 : flush(io_unit)
405 : : end if
406 : :
407 : : do ! inner loop
408 : 16285 : k=k+1
409 : 16285 : l=l+1
410 : :
411 : : ! Apply the A operator
412 [ - + ]: 16285 : if (has_omega) then
413 : 0 : call Hpsik(w,d,lambda)
414 : 0 : call Hpsik(w,cte=lambda)
415 [ # # ]: 0 : do ipw=1,npw_g
416 [ # # ]: 0 : do ii=1,2
417 : 0 : w(ii,ipw) = w(ii,ipw) + d(ii,ipw)*signe*omega**2
418 : : end do
419 : : end do
420 : : else
421 : 16285 : call Hpsik(w,d,lambda)
422 : : end if
423 : :
424 : : ! Apply projections, if requested
425 : : !if(dtset%gwcalctyp /= 2) then !This is a test to obtain the time taken by the orthos.
426 [ + - ]: 16285 : if(pow == 1) call pc_k_valence_kernel(w)
427 : : !if(pow == 2) call pc_k(w,eig_e=lambda)
428 : : !if(pow == 3) call pc_k(w,n=nband,above=.true.)
429 : : !end if
430 : :
431 : : ! Apply SQMR scheme
432 : : !tmp = scprod_k(d,w)
433 : 16285 : tmp = cg_zdotc(npw_g,d,w)
434 : 16285 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
435 : :
436 : 16285 : sigma(k-1) = tmp(1)
437 [ + + ]: 3274274 : do ipw=1,npw_g
438 [ + + ]: 9790252 : do ii=1,2
439 : 9773967 : r(ii,ipw) = r(ii,ipw)-(rho(k-1)/sigma(k-1))*w(ii,ipw)
440 : : end do
441 : : end do
442 : :
443 : : ! The following two lines must have a bug! We cannot distribute the norm this way!
444 : : ! theta(k) = norm_k(r)/g(k-1)
445 : : ! call xmpi_sum(theta(k), mpi_communicator,ierr) ! sum on all processors working on FFT!
446 : 16285 : tmp = cg_zdotc(npw_g,r,r)
447 : 16285 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
448 : 16285 : theta(k) = dsqrt(tmp(1))/g(k-1)
449 : :
450 : 16285 : c(k) = one/dsqrt(one+theta(k)**2)
451 : 16285 : g(k) = g(k-1)*theta(k)*c(k)
452 [ + + ]: 3274274 : do ipw=1,npw_g
453 [ + + ]: 9790252 : do ii=1,2
454 : 6515978 : delta(ii,ipw) = delta(ii,ipw)*(c(k)*theta(k-1))**2+d(ii,ipw)*rho(k-1)/sigma(k-1)*c(k)**2
455 : 9773967 : x(ii,ipw) = x(ii,ipw)+delta(ii,ipw)
456 : : end do
457 : : end do
458 : :
459 [ + + ]: 16285 : if (head_node) then
460 : 12677 : write(io_unit,16) k, g(k)**2
461 : 12677 : flush(io_unit)
462 : : end if
463 : :
464 : : ! Test exit condition
465 [ + + - + ]: 16285 : if(g(k)**2<tolwfr .or. k>= nline) exit
466 : : !if(k>=nline) exit
467 : :
468 : : ! Safety test every 100 iterations, check that estimated residual is of the right order of magnitude.
469 : : ! If not, restart SQMR.
470 [ - + ]: 15549 : if(mod(l,100)==0) then
471 [ # # ]: 0 : if(has_omega) then
472 : 0 : call Hpsik(w,x,lambda)
473 : 0 : call Hpsik(w,cte=lambda)
474 [ # # ]: 0 : do ipw=1,npw_g
475 [ # # ]: 0 : do ii=1,2
476 : 0 : w(ii,ipw) = w(ii,ipw) + x(ii,ipw)*signe*omega**2
477 : : end do
478 : : end do
479 : : else
480 : 0 : call Hpsik(w,x,lambda)
481 : : end if
482 : :
483 [ # # ]: 0 : if(pow == 1) call pc_k_valence_kernel(w)
484 : : !if(pow == 2) call pc_k(w,eig_e=lambda)
485 : : !if(pow == 3) call pc_k(w,n=nband,above=.true.)
486 : :
487 : : !if(norm_k(w-b)**2 > 10*g(k)**2) exit
488 [ # # ]: 0 : do ipw=1,npw_g
489 [ # # ]: 0 : do ii=1,2
490 : 0 : wmb(ii,ipw) = w(ii,ipw) - b(ii,ipw)
491 : : end do
492 : : end do
493 : 0 : tmp = cg_zdotc(npw_g,wmb,wmb)
494 : 0 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
495 [ # # ]: 0 : if(tmp(1) > 10*g(k)**2) exit
496 : :
497 : : end if
498 : :
499 : : ! Get ready for next cycle
500 : 15549 : call precondition(w,r)
501 : : !if(dtset%gwcalctyp /= 2) then
502 [ + - ]: 15549 : if(pow == 1) call pc_k_valence_kernel(w)
503 : : !if(pow == 2) call pc_k(w,eig_e=lambda)
504 : : !if(pow == 3) call pc_k(w,n=nband,above=.true.)
505 : : !end if
506 : :
507 : : !tmp = scprod_k(r,w)
508 : 15549 : tmp = cg_zdotc(npw_g,r,w)
509 : 15549 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
510 : 15549 : rho(k) = tmp(1)
511 : :
512 [ + + ]: 3126242 : do ipw=1,npw_g
513 [ + + ]: 9345420 : do ii=1,2
514 : 9329871 : d(ii,ipw) = w(ii,ipw)+d(ii,ipw)*rho(k)/rho(k-1)
515 : : end do
516 : : end do
517 : :
518 : : end do ! end inner loop
519 : :
520 : : ! Exit condition
521 [ - + - - ]: 736 : if(g(k)**2<tolwfr .or. k>=nline) exit
522 : : !if(k>=nline) exit
523 : :
524 [ # # ]: 0 : if (head_node) write(io_unit,10) " ---- RESTART of SQMR -----"
525 : :
526 : : !norm_Axb = norm_k(w-b)**2
527 [ # # ]: 0 : do ipw=1,npw_g
528 [ # # ]: 0 : do ii=1,2
529 : 0 : wmb(ii,ipw) = w(ii,ipw) - b(ii,ipw)
530 : : end do
531 : : end do
532 : 0 : tmp = cg_zdotc(npw_g,wmb,wmb)
533 : 0 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
534 : 0 : norm_Axb = tmp(1)
535 : :
536 [ # # ]: 0 : if (head_node) then
537 : 0 : write(io_unit,*) "|Ax-b|^2 :",norm_Axb
538 : 0 : write(io_unit,*) "g(k)^2 :",g(k)**2
539 : 0 : flush(io_unit)
540 : : end if
541 : :
542 : 0 : k = k+1
543 : 0 : l = 1
544 : :
545 : : ! Apply the operator
546 [ # # ]: 0 : if(has_omega) then
547 : 0 : call Hpsik(r,x,lambda)
548 : 0 : call Hpsik(r,cte=lambda)
549 [ # # ]: 0 : do ipw=1,npw_g
550 [ # # ]: 0 : do ii=1,2
551 : 0 : r(ii,ipw) = r(ii,ipw) + x(ii,ipw)*signe*omega**2
552 : : end do
553 : : end do
554 : : else
555 : 0 : call Hpsik(r,x,lambda)
556 : : end if
557 : :
558 [ # # ]: 0 : if(pow == 1) call pc_k_valence_kernel(r)
559 : : !if(pow == 2) call pc_k(r,eig_e=lambda)
560 : : !if(pow == 3) call pc_k(r,n=nband,above=.true.)
561 : :
562 [ - - ]: 736 : do ipw=1,npw_g
563 [ # # ]: 0 : do ii=1,2
564 : 0 : r(ii,ipw) = b(ii,ipw) - r(ii,ipw)
565 : : end do
566 : : end do
567 : :
568 : : end do ! outer loop
569 : :
570 : :
571 : 736 : ktot = ktot+k
572 [ - + - - ]: 736 : if(k >= nline .and. head_node ) then
573 : 0 : write(io_unit,10) " **** Iterations were not enough to converge! ****"
574 : : end if
575 : :
576 [ - + ]: 736 : if( present(kill_Pc_x) ) then
577 [ # # # # ]: 0 : if (.not. kill_Pc_x .and. pow == 1) call pc_k_valence_kernel(x)
578 : : end if
579 : :
580 [ + - ]: 736 : if( .not. present(kill_Pc_x) .and. pow == 1 ) call pc_k_valence_kernel(x)
581 : :
582 : :
583 : :
584 : :
585 : : !if(pow == 2) call pc_k(x,eig_e=lambda)
586 : : !if(pow == 3) call pc_k(x,n=nband,above=.true.)
587 : :
588 [ - + ]: 736 : if(has_omega) then
589 : 0 : call Hpsik(w,x,lambda)
590 : 0 : call Hpsik(w,cte=lambda)
591 [ # # ]: 0 : do ipw=1,npw_g
592 [ # # ]: 0 : do ii=1,2
593 : 0 : w(ii,ipw) = w(ii,ipw) + x(ii,ipw)*signe*omega**2
594 : : end do
595 : : end do
596 : : else
597 : 736 : call Hpsik(w,x,lambda)
598 : : end if
599 [ + - ]: 736 : if(pow == 1) call pc_k_valence_kernel(w)
600 : : !if(pow == 2) call pc_k(w,eig_e=lambda)
601 : : !if(pow == 3) call pc_k(w,n=nband,above=.true.)
602 : :
603 : : !residual = norm_k(w-b)**2
604 [ + + ]: 148768 : do ipw=1,npw_g
605 [ + + ]: 444832 : do ii=1,2
606 : 444096 : wmb(ii,ipw) = w(ii,ipw) - b(ii,ipw)
607 : : end do
608 : : end do
609 : 736 : tmp = cg_zdotc(npw_g,wmb,wmb)
610 : 736 : call xmpi_sum(tmp, mpi_communicator,ierr) ! sum on all processors working on FFT!
611 : 736 : residual = tmp(1)
612 : :
613 [ + + ]: 736 : if (head_node) then
614 : 576 : write(io_unit,15) "iterations :", k
615 : 576 : write(io_unit,14) "tolwfr :", tolwfr
616 : 576 : write(io_unit,14) "residuals (estimated) :", g(k)**2
617 : 576 : write(io_unit,14) "residuals : |Ax-b|^2 :", residual
618 : 576 : close(io_unit)
619 : : end if
620 : :
621 : :
622 : :
623 : 736 : ABI_FREE(g)
624 : 736 : ABI_FREE(theta)
625 : 736 : ABI_FREE(rho)
626 : 736 : ABI_FREE(sigma)
627 : 736 : ABI_FREE(c)
628 [ - + ]: 736 : ABI_FREE(t)
629 [ - + ]: 736 : ABI_FREE(delta)
630 [ - + ]: 736 : ABI_FREE(r)
631 [ - + ]: 736 : ABI_FREE(d)
632 [ - + ]: 736 : ABI_FREE(w)
633 [ - + ]: 736 : ABI_FREE(wmb)
634 : :
635 : :
636 : 736 : OPTION_TIMAB = 2
637 : 736 : call timab(GWLS_TIMAB,OPTION_TIMAB,tsec)
638 : :
639 : :
640 : :
641 : : 10 format(A)
642 : : 11 format(A,I8)
643 : : 12 format(A,E24.16,A)
644 : : 13 format(A,I2,A)
645 : : 14 format(20X,A,E24.16)
646 : : 15 format(20X,A,I8)
647 : : 16 format(5X,I5,15X,E12.3)
648 : :
649 : : end subroutine sqmr
650 : : !!***
651 : :
652 : : !!****f* m_hamiltonian/qmr
653 : : !! NAME
654 : : !! qmr
655 : : !!
656 : : !! FUNCTION
657 : : !! .
658 : : !!
659 : : !! INPUTS
660 : : !!
661 : : !! OUTPUT
662 : : !!
663 : : !! SOURCE
664 : :
665 : 320 : subroutine qmr(b,x,lambda) !,project_on_what)
666 : : !--------------------------------------------------------------------------------
667 : : ! This subroutine solves the linear algebra problem
668 : : !
669 : : ! A x = b
670 : : !
671 : : ! where A := (H - lambda) can be non-hermitian.
672 : : ! Thus, complex values of lambda are allowed and
673 : : ! non-hermitian operators could be handled instead of H.
674 : : !
675 : : ! Arguments :
676 : : ! INPUT
677 : : ! -----
678 : : ! real(dp) b(2,npw_k) right-hand-side of the equation to be solved
679 : : ! real(dp) lambda(2) value to be subtracted from the Hamiltonian (complex).
680 : : ! integer project_on_what flag which determines the projection scheme.
681 : : !
682 : : ! OUTPUT
683 : : ! -----
684 : : ! real(dp) x(2,npw_k) solution
685 : : !
686 : : ! project_on_what action
687 : : ! ------------------------------------------------------
688 : : ! 0 no projection
689 : : ! 1 projection on conduction states
690 : : ! 2 projection out of subspace degenerate with lambda
691 : : ! 3 projection on states beyond all the states explicitly stored
692 : : !--------------------------------------------------------------------------------
693 : :
694 : : !External variables
695 : : real(dp), intent(in) :: b(2,npw_k)
696 : : real(dp), intent(in) :: lambda(2)
697 : : real(dp), intent(out) :: x(2,npw_k)
698 : : !integer, intent(in) :: project_on_what !Unused yet, no projections done.
699 : :
700 : : !Local variables
701 : 320 : complex(dp), allocatable :: xc(:), r(:), v(:), w(:), z(:), p(:), q(:), y(:), t(:), d(:), s(:)
702 : 320 : complex(dp), allocatable :: beta(:), eta(:), delta(:), epsilonn(:)
703 : : complex(dp) :: lambdac
704 : 320 : real(dp), allocatable :: rho(:), zeta(:), gama(:), theta(:), resid(:)
705 : : integer :: i
706 : : integer :: ierr
707 : :
708 : : integer :: mpi_communicator
709 : : !logical :: precondition_on
710 : :
711 : : ! *************************************************************************
712 : :
713 [ + + + + : 247040 : if(sum(b**2) < tol12) then
+ + ]
714 [ + + + + ]: 148224 : x=zero
715 : : else
716 : :
717 : : !Allocation
718 [ - + + - : 384 : ABI_MALLOC(xc,(npw_k))
+ - + - -
+ - + ]
719 [ - + + - : 256 : ABI_MALLOC(r ,(npw_k))
- + - + ]
720 [ - + + - : 256 : ABI_MALLOC(v ,(npw_k))
- + - + ]
721 [ - + + - : 256 : ABI_MALLOC(w ,(npw_k))
- + - + ]
722 [ - + + - : 256 : ABI_MALLOC(z ,(npw_k))
- + - + ]
723 [ - + + - : 256 : ABI_MALLOC(p ,(npw_k))
- + - + ]
724 [ - + + - : 256 : ABI_MALLOC(q ,(npw_k))
- + - + ]
725 [ - + + - : 256 : ABI_MALLOC(y ,(npw_k))
- + - + ]
726 [ - + + - : 256 : ABI_MALLOC(t ,(npw_k))
- + - + ]
727 [ - + + - : 256 : ABI_MALLOC(d ,(npw_k))
- + - + ]
728 [ - + + - : 256 : ABI_MALLOC(s ,(npw_k))
- + - + ]
729 : :
730 [ - + + - : 384 : ABI_MALLOC(beta ,(nline))
+ - + - -
+ - + ]
731 [ + - - + : 384 : ABI_MALLOC(rho ,(nline+1))
+ - + - +
- - + -
+ ]
732 [ - + + - : 256 : ABI_MALLOC(zeta ,(nline+1))
- + - + ]
733 [ - + + - : 256 : ABI_MALLOC(gama ,(nline+1))
- + - + ]
734 [ - + + - : 384 : ABI_MALLOC(eta ,(nline+1))
+ - + - -
+ - + ]
735 [ - + + - : 256 : ABI_MALLOC(theta ,(nline+1))
- + - + ]
736 [ - + + - : 256 : ABI_MALLOC(delta ,(nline))
- + - + ]
737 [ - + + - : 256 : ABI_MALLOC(epsilonn,(nline))
- + - + ]
738 [ - + + - : 256 : ABI_MALLOC(resid ,(nline+1))
- + - + ]
739 : :
740 : : !Initialization
741 [ + + + + ]: 98816 : x = zero
742 [ + + ]: 33024 : xc = zero
743 [ + + ]: 33024 : r = zero
744 [ + + ]: 33024 : v = zero
745 [ + + ]: 33024 : w = zero
746 [ + + ]: 33024 : z = zero
747 [ + + ]: 33024 : p = zero
748 [ + + ]: 33024 : q = zero
749 [ + + ]: 33024 : y = zero
750 [ + + ]: 33024 : t = zero
751 [ + + ]: 33024 : d = zero
752 [ + + ]: 33024 : s = zero
753 : :
754 [ + + ]: 12928 : beta = zero
755 [ + + ]: 13056 : rho = zero
756 [ + + ]: 13056 : zeta = zero
757 [ + + ]: 13056 : gama = zero
758 [ + + ]: 13056 : eta = zero
759 [ + + ]: 13056 : theta = zero
760 [ + + ]: 12928 : delta = zero
761 [ + + ]: 12928 : epsilonn = zero
762 [ + + ]: 13056 : resid = zero
763 : :
764 : 128 : call unset_precondition()
765 : :
766 : :
767 : : !mpi_communicator = mpi_enreg%comm_fft
768 : 128 : mpi_communicator = mpi_enreg%comm_bandfft
769 : :
770 : 128 : lambdac = dcmplx(lambda(1),lambda(2))
771 : :
772 : 128 : i = 1
773 [ - + + + ]: 33152 : r = dcmplx(b(1,:),b(2,:))
774 [ - + + + ]: 33152 : v = r
775 : :
776 : 128 : rho(i) = norm_kc(v)
777 : 128 : call xmpi_sum(rho(i),mpi_communicator ,ierr) ! sum on all processors working on FFT!
778 : :
779 [ - + + + ]: 33152 : w = r
780 : 128 : call precondition_cplx(z,w)
781 : 128 : zeta(i) = norm_kc(z)
782 : 128 : call xmpi_sum(zeta(i),mpi_communicator,ierr) ! sum on all processors working on FFT!
783 : :
784 : 128 : gama(i) = one
785 : 128 : eta(i) = -one
786 : : !theta(i) = zero
787 : : !p = zero
788 : : !q = zero
789 : :
790 [ - + ]: 3612 : do i=1,nline
791 [ + + ]: 931896 : v = v/rho(i)
792 [ + + ]: 931896 : w = w/zeta(i)
793 [ + + ]: 931896 : z = z/zeta(i)
794 : 3612 : delta(i) = scprod_kc(z,v)
795 : 3612 : call xmpi_sum(delta(i),mpi_communicator,ierr) ! sum on all processors working on FFT!
796 : :
797 : 3612 : call precondition_cplx(y,v)
798 [ + + ]: 3612 : if(i/=1) then
799 [ - + + + ]: 902356 : p = y - (zeta(i)*delta(i)/epsilonn(i-1))*p
800 [ - + + + ]: 902356 : q = z - ( rho(i)*delta(i)/epsilonn(i-1))*q
801 : : else
802 [ - + + + ]: 33152 : p = y
803 [ - + + + ]: 33152 : q = z
804 : : end if
805 : 3612 : call Hpsikc(t,p,lambdac)
806 : 3612 : epsilonn(i) = scprod_kc(q,t)
807 : 3612 : call xmpi_sum(epsilonn(i),mpi_communicator,ierr) ! sum on all processors working on FFT!
808 : :
809 : 3612 : beta(i) = epsilonn(i)/delta(i)
810 [ - + + + ]: 935508 : v = t - beta(i)*v
811 : 3612 : rho(i+1) = norm_kc(v)
812 : 3612 : call xmpi_sum(rho(i+1),mpi_communicator,ierr) ! sum on all processors working on FFT!
813 : :
814 [ - + + + ]: 935508 : call Hpsikc(z,q,conjg(lambdac)) ; w = z - beta(i)*w
815 : 3612 : call precondition_cplx(z,w)
816 : 3612 : zeta(i+1) = norm_kc(z)
817 : 3612 : call xmpi_sum(zeta(i+1), mpi_communicator,ierr) ! sum on all processors working on FFT!
818 : :
819 : 3612 : theta(i+1) = rho(i+1)/(gama(i)*abs(beta(i)))
820 : 3612 : gama(i+1) = 1./sqrt(1+theta(i+1)**2)
821 : 3612 : eta(i+1) = -eta(i)*rho(i)*gama(i+1)**2/(beta(i)*gama(i)**2)
822 [ - + + + ]: 935508 : d = eta(i+1)*p + ((theta(i)*gama(i+1))**2)*d
823 [ - + + + ]: 935508 : s = eta(i+1)*t + ((theta(i)*gama(i+1))**2)*s
824 [ - + + + ]: 935508 : xc = xc + d
825 [ - + + + ]: 935508 : r = r - s
826 : 3612 : resid(i) = norm_kc(r)**2
827 : 3612 : call xmpi_sum(resid(i), mpi_communicator,ierr) ! sum on all processors working on FFT!
828 : :
829 : : !write(std_out,*) "QMR residual**2 = ",resid(i),"; i = ",i-1
830 [ + + ]: 18060 : if(resid(i) < tolwfr) exit
831 : : end do
832 : :
833 [ - + ]: 128 : if(i>=nline) then
834 : 0 : write(std_out,*) " **** Iterations were not enough to converge! ****"
835 : : end if
836 : :
837 : 128 : call Hpsikc(r,xc,lambdac)
838 [ - + + + ]: 33152 : v = dcmplx(b(1,:),b(2,:))
839 [ + - + + ]: 33024 : resid(nline+1) = norm_kc(r - v)**2
840 : 128 : call xmpi_sum(resid(nline+1), mpi_communicator,ierr) ! sum on all processors working on FFT!
841 : :
842 : 128 : write(std_out,*) "QMR residual**2 (at end) = ",resid(nline+1),"; # iterations = ",i-1
843 : :
844 [ + + ]: 33024 : x(1,:) = dble(xc)
845 [ + + ]: 33024 : x(2,:) = dimag(xc)
846 : :
847 : : !Deallocate
848 [ - + ]: 128 : ABI_FREE(xc)
849 [ - + ]: 128 : ABI_FREE(r)
850 [ - + ]: 128 : ABI_FREE(v)
851 [ - + ]: 128 : ABI_FREE(w)
852 [ - + ]: 128 : ABI_FREE(z)
853 [ - + ]: 128 : ABI_FREE(p)
854 [ - + ]: 128 : ABI_FREE(q)
855 [ - + ]: 128 : ABI_FREE(y)
856 [ - + ]: 128 : ABI_FREE(t)
857 [ - + ]: 128 : ABI_FREE(d)
858 [ - + ]: 128 : ABI_FREE(s)
859 : :
860 [ - + ]: 128 : ABI_FREE(beta )
861 : 128 : ABI_FREE(rho )
862 : 128 : ABI_FREE(zeta )
863 : 128 : ABI_FREE(gama )
864 : 128 : ABI_FREE(eta )
865 [ - + ]: 128 : ABI_FREE(theta )
866 [ - + ]: 128 : ABI_FREE(delta )
867 [ - + ]: 128 : ABI_FREE(epsilonn)
868 : 512 : ABI_FREE(resid )
869 : :
870 : : end if
871 : :
872 : 320 : end subroutine qmr
873 : : !!***
874 : :
875 : : end module m_gwls_lineqsolver
876 : : !!***
|
