Branch data Line data Source code
1 : : !!****m* ABINIT/m_gwls_LanczosResolvents
2 : : !! NAME
3 : : !! m_gwls_LanczosResolvents
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 : : module m_gwls_LanczosResolvents
24 : : !----------------------------------------------------------------------------------------------------
25 : : ! This module contains routines to compute the matrix elements of the resolent, namely
26 : : !
27 : : ! < phi_l | 1/ (H-z) | psi_l' >
28 : : !
29 : : ! where H is the Hamiltonian, and z should be such that we are far from its poles.
30 : : !
31 : : !----------------------------------------------------------------------------------------------------
32 : : ! local modules
33 : :
34 : : use m_gwls_utility
35 : : use m_gwls_wf
36 : : use m_gwls_TimingLog
37 : : use m_gwls_hamiltonian
38 : : use m_gwls_lineqsolver
39 : : use m_gwls_GWlanczos
40 : :
41 : : use defs_basis
42 : : use m_abicore
43 : : use m_xmpi
44 : : use m_io_tools, only : get_unit
45 : :
46 : : implicit none
47 : : save
48 : : private
49 : : !!***
50 : :
51 : :
52 : : integer :: LR_kmax, LR_nseeds
53 : :
54 : : complex(dp), allocatable :: precondition_C(:) ! this operator is diagonal!
55 : : complex(dp), allocatable :: precondition_one_on_C(:) ! this operator is diagonal!
56 : :
57 : :
58 : : complex(dp), allocatable, public :: Hamiltonian_Qk(:,:) ! Lanczos basis of the Hamiltonian
59 : :
60 : : complex(dp), allocatable :: LR_alpha(:,:,:)
61 : : complex(dp), allocatable :: LR_beta (:,:,:)
62 : : complex(dp), allocatable :: LR_seeds(:,:)
63 : : real(dp), allocatable :: LR_Hamiltonian_eigenvalues(:)
64 : :
65 : : complex(dp), allocatable, public :: LR_M_matrix(:,:)
66 : : !!***
67 : :
68 : : public :: setup_LanczosResolvents
69 : : public :: cleanup_LanczosResolvents
70 : : public :: compute_resolvent_column_shift_lanczos
71 : : public :: compute_resolvent_column_shift_lanczos_right_vectors
72 : : public :: build_preconditioned_Hamiltonian_Lanczos_basis
73 : : public :: invert_general_matrix
74 : : !!***
75 : :
76 : : contains
77 : :
78 : :
79 : : !!****f* m_hamiltonian/setup_LanczosResolvents
80 : : !! NAME
81 : : !! setup_LanczosResolvents
82 : : !!
83 : : !! FUNCTION
84 : : !! .
85 : : !!
86 : : !! INPUTS
87 : : !!
88 : : !! OUTPUT
89 : : !!
90 : : !! SOURCE
91 : :
92 : 26 : subroutine setup_LanczosResolvents(kmax, prec)
93 : : !----------------------------------------------------------------------------------------------------
94 : : !
95 : : ! This subroutine prepares global work arrays in order to use the Lanczos scheme to
96 : : ! compute matrix elements of inverse operators.
97 : : !
98 : : !
99 : : !----------------------------------------------------------------------------------------------------
100 : : implicit none
101 : :
102 : : !------------------------------
103 : : ! input/output variables
104 : : !------------------------------
105 : :
106 : : integer, intent(in) :: kmax ! number of Lanczos steps
107 : : logical, intent(in) :: prec ! should there be preconditioning?
108 : :
109 : : ! *************************************************************************
110 : :
111 : :
112 : 26 : LR_kmax = kmax
113 : 26 : LR_nseeds = 1 ! only one seed
114 : :
115 : : ! Prepare the algorithm
116 : :
117 [ - + + - : 78 : ABI_MALLOC(precondition_C, (npw_g))
+ - + - -
+ - + -
+ ]
118 [ - + + - : 52 : ABI_MALLOC(precondition_one_on_C, (npw_g))
- + - + -
+ ]
119 : :
120 : :
121 [ - + ]: 26 : if ( prec ) then
122 : 0 : call set_precondition()
123 : :
124 [ # # ]: 0 : precondition_one_on_C(:) = sqrt(pcon(:))
125 : : !precondition_one_on_C(:) = pcon(:) ! yields very bad results!
126 [ # # ]: 0 : precondition_C(:) = cmplx_1/precondition_one_on_C(:)
127 : :
128 : : else
129 [ + + ]: 4652 : precondition_C(:) = cmplx_1
130 [ + + ]: 4652 : precondition_one_on_C(:) = cmplx_1
131 : : end if
132 : :
133 [ - + + - : 130 : ABI_MALLOC(LR_alpha, (LR_nseeds,LR_nseeds,LR_kmax))
- + + - -
+ + - + -
+ - - + -
+ - + ]
134 [ - + + - : 104 : ABI_MALLOC(LR_beta , (LR_nseeds,LR_nseeds,LR_kmax))
- + + - -
+ + - - +
- + - + ]
135 [ - + + - : 104 : ABI_MALLOC(LR_seeds, (npw_g,LR_nseeds))
- + + - +
- + - - +
- + - + ]
136 [ - + + - : 104 : ABI_MALLOC(Hamiltonian_Qk, (npw_g,LR_kmax))
- + + - +
- + - - +
- + - + ]
137 [ - + + - : 78 : ABI_MALLOC(LR_Hamiltonian_eigenvalues, (LR_kmax))
+ - + - -
+ - + -
+ ]
138 [ - + + - : 104 : ABI_MALLOC(LR_M_matrix, (LR_kmax,LR_kmax))
- + + - +
- + - - +
- + - + ]
139 : :
140 : 26 : end subroutine setup_LanczosResolvents
141 : : !!***
142 : :
143 : : !!****f* m_hamiltonian/cleanup_LanczosResolvents
144 : : !! NAME
145 : : !! cleanup_LanczosResolvents
146 : : !!
147 : : !! FUNCTION
148 : : !! .
149 : : !!
150 : : !! INPUTS
151 : : !!
152 : : !! OUTPUT
153 : : !!
154 : : !! SOURCE
155 : :
156 : 26 : subroutine cleanup_LanczosResolvents
157 : : !----------------------------------------------------------------------------------------------------
158 : : !
159 : : ! This subroutine cleans up global arrays.
160 : : !
161 : : !
162 : : !----------------------------------------------------------------------------------------------------
163 : : implicit none
164 : :
165 : : ! *************************************************************************
166 : :
167 [ - + ]: 26 : ABI_FREE(precondition_C)
168 [ - + ]: 26 : ABI_FREE(precondition_one_on_C)
169 [ - + ]: 26 : ABI_FREE(LR_alpha)
170 [ - + ]: 26 : ABI_FREE(LR_beta )
171 [ - + ]: 26 : ABI_FREE(LR_seeds)
172 [ - + ]: 26 : ABI_FREE(Hamiltonian_Qk)
173 [ - + ]: 26 : ABI_FREE(LR_Hamiltonian_eigenvalues)
174 [ - + ]: 26 : ABI_FREE(LR_M_matrix)
175 : :
176 : 26 : end subroutine cleanup_LanczosResolvents
177 : : !!***
178 : :
179 : : !!****f* m_hamiltonian/matrix_function_preconditioned_Hamiltonian
180 : : !! NAME
181 : : !! matrix_function_preconditioned_Hamiltonian
182 : : !!
183 : : !! FUNCTION
184 : : !! .
185 : : !!
186 : : !! INPUTS
187 : : !!
188 : : !! OUTPUT
189 : : !!
190 : : !! SOURCE
191 : :
192 : 1728 : subroutine matrix_function_preconditioned_Hamiltonian(vector_out,vector_in,Hsize)
193 : : !----------------------------------------------------------------------------------------------------
194 : : ! This subroutine is a simple wrapper around the preconditioned Hamiltonian operator which acts as
195 : : !
196 : : ! C^{-1} . H . C^{-1}
197 : : !
198 : : ! where C^{-2} ~ 1 / T, where T is the kinetic energy.
199 : : !
200 : : !----------------------------------------------------------------------------------------------------
201 : : implicit none
202 : : integer, intent(in) :: Hsize
203 : : complex(dp), intent(out) :: vector_out(Hsize)
204 : : complex(dp), intent(in) :: vector_in(Hsize)
205 : :
206 : : ! local variables
207 : 3456 : real(dp) :: psikg(2,npw_g)
208 : :
209 : 3456 : complex(dp) :: tmp_vector(Hsize)
210 : :
211 : : ! *************************************************************************
212 : :
213 : : ! convert from one format to the other
214 : :
215 [ + + ]: 347136 : tmp_vector(:) = precondition_one_on_C(:)*vector_in(:)
216 : :
217 [ + + ]: 347136 : psikg(1,:) = dble (tmp_vector(:))
218 [ + + ]: 347136 : psikg(2,:) = dimag(tmp_vector(:))
219 : :
220 : 1728 : call Hpsik(psikg)
221 : :
222 [ + + ]: 347136 : tmp_vector(:) = cmplx_1*psikg(1,:)+cmplx_i*psikg(2,:)
223 : :
224 [ + + ]: 347136 : vector_out(:) = precondition_one_on_C(:)*tmp_vector(:)
225 : :
226 : :
227 : 1728 : end subroutine matrix_function_preconditioned_Hamiltonian
228 : : !!***
229 : :
230 : : !!****f* m_hamiltonian/build_preconditioned_Hamiltonian_Lanczos_basis
231 : : !! NAME
232 : : !! build_preconditioned_Hamiltonian_Lanczos_basis
233 : : !!
234 : : !! FUNCTION
235 : : !! .
236 : : !!
237 : : !! INPUTS
238 : : !!
239 : : !! OUTPUT
240 : : !!
241 : : !! SOURCE
242 : :
243 : 288 : subroutine build_preconditioned_Hamiltonian_Lanczos_basis(seed_vector)
244 : : !----------------------------------------------------------------------------------------------------
245 : : ! This function Computes the Lanczos basis of the preconditioned Hamiltonian.
246 : : !----------------------------------------------------------------------------------------------------
247 : : implicit none
248 : : complex(dp), intent(in) :: seed_vector(npw_g)
249 : :
250 : : integer :: l
251 : : integer :: ierr
252 : : integer :: mpi_communicator
253 : :
254 : : ! *************************************************************************
255 : :
256 : :
257 : :
258 : 144 : mpi_communicator = mpi_enreg%comm_fft
259 : :
260 : : ! Compute the Lanczos basis as well as the eigenvalues of the T matrix
261 : : ! Seed must also be preconditioned!
262 [ + + ]: 28928 : LR_seeds(:,1) = precondition_one_on_C(:)*seed_vector(:)
263 : :
264 : : call block_lanczos_algorithm(mpi_communicator, matrix_function_preconditioned_Hamiltonian, LR_kmax, LR_nseeds, npw_g, LR_seeds, &
265 : 144 : & LR_alpha, LR_beta, Hamiltonian_Qk)
266 : :
267 : :
268 : : call diagonalize_lanczos_banded(LR_kmax,LR_nseeds,npw_g, LR_alpha,LR_beta, &
269 : 144 : Hamiltonian_Qk, LR_Hamiltonian_eigenvalues,.false.)
270 : :
271 : :
272 : :
273 : : ! Update the Lanczos vectors with the preconditioning matrix
274 [ + + ]: 1872 : do l = 1, LR_kmax
275 [ + + ]: 347280 : Hamiltonian_Qk(:,l) = precondition_C(:)*Hamiltonian_Qk(:,l)
276 : : end do
277 : :
278 : : ! compute the M matrix
279 : :
280 : : ! Compute Q^dagger . Q
281 : 144 : call ZGEMM('C','N',LR_kmax,LR_kmax,npw_g,cmplx_1,Hamiltonian_Qk,npw_g,Hamiltonian_Qk,npw_g,cmplx_0,LR_M_matrix,LR_kmax)
282 : 144 : call xmpi_sum(LR_M_matrix,mpi_communicator,ierr) ! sum on all processors working on FFT!
283 : :
284 [ + + ]: 1872 : do l = 1, LR_kmax
285 [ + + ]: 24912 : LR_M_matrix(l,:) = LR_M_matrix(l,:)*LR_Hamiltonian_eigenvalues(l)
286 : : end do
287 : :
288 : 144 : end subroutine build_preconditioned_Hamiltonian_Lanczos_basis
289 : : !!***
290 : :
291 : :
292 : : !!****f* m_hamiltonian/compute_resolvent_column_shift_lanczos
293 : : !! NAME
294 : : !! compute_resolvent_column_shift_lanczos
295 : : !!
296 : : !! FUNCTION
297 : : !! .
298 : : !!
299 : : !! INPUTS
300 : : !!
301 : : !! OUTPUT
302 : : !!
303 : : !! SOURCE
304 : :
305 : 72 : subroutine compute_resolvent_column_shift_lanczos(nz, list_z, nvec, list_left_vectors, seed_vector, matrix_elements_resolvent)
306 : : !----------------------------------------------------------------------------------------------------
307 : : ! This function Computes one column of the resolvent using the shift lanczos algorithm,
308 : : ! namely
309 : : !
310 : : !
311 : : ! I_l(z) = < left_vector_l | xi(z) >
312 : : !
313 : : ! where |xi(z) > is obtained from the shift lanczos scheme.
314 : : !----------------------------------------------------------------------------------------------------
315 : : implicit none
316 : :
317 : : integer , intent(in) :: nz
318 : : complex(dp), intent(in) :: list_z(nz)
319 : :
320 : : integer , intent(in) :: nvec
321 : : complex(dp), intent(in) :: list_left_vectors(npw_g,nvec)
322 : :
323 : : complex(dp), intent(in) :: seed_vector(npw_g)
324 : :
325 : : complex(dp), intent(out) :: matrix_elements_resolvent(nz,nvec)
326 : :
327 : : ! local variables
328 : : integer :: iz, l
329 : : complex(dp) :: z
330 : : integer :: ierr
331 : : integer :: mpi_communicator
332 : :
333 : 144 : complex(dp) :: right_vec(LR_kmax)
334 : 144 : complex(dp) :: left_vecs(LR_kmax,nvec)
335 : 144 : complex(dp) :: work_vec(LR_kmax)
336 : 144 : complex(dp) :: shift_lanczos_matrix(LR_kmax,LR_kmax)
337 : :
338 : 144 : complex(dp) :: work_array(npw_g)
339 : :
340 : : ! *************************************************************************
341 : :
342 : :
343 : : ! processes will communicate along FFT rows
344 : 72 : mpi_communicator = mpi_enreg%comm_fft
345 : :
346 : : !----------------------------------------------------------------------------------------------------
347 : : ! First, build the Lanczos basis for the preconditioned Hamiltonian
348 : : !----------------------------------------------------------------------------------------------------
349 : 72 : call build_preconditioned_Hamiltonian_Lanczos_basis(seed_vector)
350 : :
351 : :
352 : : !----------------------------------------------------------------------------------------------------
353 : : ! Next, generate arrays which do not depend on z, the external shift.
354 : : !----------------------------------------------------------------------------------------------------
355 : :
356 : : ! Q^dagger . C^{-2} | seed >
357 [ + + ]: 14464 : work_array(:) = precondition_one_on_C(:)**2*seed_vector
358 : :
359 : 72 : call ZGEMV('C',npw_g,LR_kmax,cmplx_1,Hamiltonian_Qk,npw_g,work_array,1,cmplx_0,right_vec,1)
360 : 72 : call xmpi_sum(right_vec,mpi_communicator,ierr) ! sum on all processors working on FFT!
361 : :
362 : : ! Q^dagger | left_vectors >
363 : 72 : call ZGEMM('C','N',LR_kmax,nvec,npw_g,cmplx_1,Hamiltonian_Qk,npw_g,list_left_vectors,npw_g,cmplx_0,left_vecs,LR_kmax)
364 : 72 : call xmpi_sum(left_vecs,mpi_communicator,ierr) ! sum on all processors working on FFT!
365 : :
366 : : !----------------------------------------------------------------------------------------------------
367 : : ! Use shift Lanczos to compute all matrix elements!
368 : : !----------------------------------------------------------------------------------------------------
369 : :
370 [ + + ]: 144 : do iz = 1, nz
371 : :
372 : 72 : z = list_z(iz)
373 : :
374 : : ! Generate the matrix to be inverted
375 [ + + + + ]: 5256 : shift_lanczos_matrix(:,:) = LR_M_matrix(:,:)
376 : :
377 [ + + ]: 648 : do l = 1, LR_kmax
378 : 648 : shift_lanczos_matrix(l,l) = shift_lanczos_matrix(l,l)-z
379 : : end do
380 : :
381 : : ! since z could be complex, the matrix is not necessarily hermitian. Invert using general Lapack scheme
382 : 72 : call invert_general_matrix(LR_kmax,shift_lanczos_matrix)
383 : : ! the matrix now contains the inverse!
384 : :
385 : :
386 : : ! | work_vec > = M^{-1} . | right_vec >
387 : 72 : call ZGEMV('N',LR_kmax,LR_kmax,cmplx_1,shift_lanczos_matrix,LR_kmax,right_vec,1,cmplx_0,work_vec,1)
388 : :
389 : :
390 : : ! matrix_elements = < right_vecs | work_vec >
391 [ - + - - : 144 : call ZGEMV('C',LR_kmax,nvec,cmplx_1,left_vecs,LR_kmax,work_vec,1,cmplx_0,matrix_elements_resolvent(iz,:),1)
- - - + -
- ]
392 : :
393 : : end do
394 : :
395 : :
396 : 72 : end subroutine compute_resolvent_column_shift_lanczos
397 : : !!***
398 : :
399 : :
400 : : !!****f* m_hamiltonian/compute_resolvent_column_shift_lanczos_right_vectors
401 : : !! NAME
402 : : !! compute_resolvent_column_shift_lanczos_right_vectors
403 : : !!
404 : : !! FUNCTION
405 : : !! .
406 : : !!
407 : : !! INPUTS
408 : : !!
409 : : !! OUTPUT
410 : : !!
411 : : !! SOURCE
412 : :
413 : 72 : subroutine compute_resolvent_column_shift_lanczos_right_vectors(seed_vector, right_vec)
414 : : !----------------------------------------------------------------------------------------------------
415 : : ! The goal of this routine is to participate in the computation of a column of the resolvent
416 : : ! using the shift lanczos algorithm. This routine should be called FIRST.
417 : : !
418 : : ! This routine:
419 : : !
420 : : ! 1) builds the preconditioned Hamiltonian Lanczos basis, using the seed_vector, and
421 : : ! stores result in the Module variables.
422 : : !
423 : : ! 2) prepares and returns the Lanczos basis-projected right vectors, ready for further processing.
424 : : !
425 : : !
426 : : !----------------------------------------------------------------------------------------------------
427 : : implicit none
428 : :
429 : : complex(dp), intent(in) :: seed_vector(npw_g)
430 : : complex(dp), intent(out):: right_vec(LR_kmax)
431 : :
432 : : ! local variables
433 : : integer :: mpi_communicator
434 : : integer :: ierr
435 : :
436 : 144 : complex(dp) :: work_array(npw_g)
437 : :
438 : : ! *************************************************************************
439 : :
440 : :
441 : : ! processes will communicate along FFT rows
442 : 72 : mpi_communicator = mpi_enreg%comm_fft
443 : :
444 : : !----------------------------------------------------------------------------------------------------
445 : : ! First, build the Lanczos basis for the preconditioned Hamiltonian
446 : : !----------------------------------------------------------------------------------------------------
447 : 72 : call build_preconditioned_Hamiltonian_Lanczos_basis(seed_vector)
448 : :
449 : : !----------------------------------------------------------------------------------------------------
450 : : ! Next, generate arrays which do not depend on z, the external shift.
451 : : !----------------------------------------------------------------------------------------------------
452 : :
453 : : ! Q^dagger . C^{-2} | seed >
454 [ + + ]: 14464 : work_array(:) = precondition_one_on_C(:)**2*seed_vector(:)
455 : :
456 : 72 : call ZGEMV('C', npw_g, LR_kmax, cmplx_1, Hamiltonian_Qk, npw_g, work_array, 1, cmplx_0, right_vec, 1)
457 : 72 : call xmpi_sum(right_vec,mpi_communicator,ierr) ! sum on all processors working on FFT!
458 : :
459 : :
460 : 72 : end subroutine compute_resolvent_column_shift_lanczos_right_vectors
461 : : !!***
462 : :
463 : :
464 : :
465 : : !!****f* m_hamiltonian/invert_general_matrix
466 : : !! NAME
467 : : !! invert_general_matrix
468 : : !!
469 : : !! FUNCTION
470 : : !! .
471 : : !!
472 : : !! INPUTS
473 : : !!
474 : : !! OUTPUT
475 : : !!
476 : : !! SOURCE
477 : :
478 : 2152 : subroutine invert_general_matrix(n,matrix)
479 : : !----------------------------------------------------------------------------------------------------
480 : : ! Simple wrapper around lapack routines to invert a general matrix
481 : : !----------------------------------------------------------------------------------------------------
482 : : implicit none
483 : :
484 : : integer, intent(in) :: n
485 : : complex(dp), intent(inout) :: matrix(n,n)
486 : :
487 : : integer :: info
488 : 4304 : integer :: ipiv(n)
489 : :
490 : : integer :: debug_unit
491 : : character(50) :: debug_filename
492 : 4304 : complex(dp) :: work(n)
493 : :
494 : : ! *************************************************************************
495 : :
496 : 2152 : call ZGETRF( n, n, matrix, n, ipiv, info )
497 [ - + ]: 2152 : if ( info /= 0) then
498 : 0 : debug_unit = get_unit()
499 : 0 : write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
500 : :
501 : 0 : open(debug_unit,file=trim(debug_filename),status='unknown')
502 : :
503 : 0 : write(debug_unit,'(A)') '*********************************************************************************************'
504 : 0 : write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZGETRF(1), gwls_LanczosResolvents'
505 : 0 : write(debug_unit,'(A)') '*********************************************************************************************'
506 : :
507 : 0 : close(debug_unit)
508 : :
509 : : end if
510 : :
511 : :
512 : 2152 : call ZGETRI( n, matrix, n, ipiv, work, n, info )
513 : :
514 [ - + ]: 2152 : if ( info /= 0) then
515 : 0 : debug_unit = get_unit()
516 : 0 : write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
517 : :
518 : 0 : open(debug_unit,file=trim(debug_filename),status='unknown')
519 : :
520 : 0 : write(debug_unit,'(A)') '*********************************************************************************************'
521 : 0 : write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZGETRI(1), gwls_LanczosResolvents'
522 : 0 : write(debug_unit,'(A)') '*********************************************************************************************'
523 : :
524 : 0 : close(debug_unit)
525 : :
526 : : end if
527 : :
528 : :
529 : 2152 : end subroutine invert_general_matrix
530 : : !!***
531 : :
532 : :
533 : : end module m_gwls_LanczosResolvents
534 : : !!***
|
