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 : !!***
|
