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