TABLE OF CONTENTS

ABINIT/m_gwls_LanczosResolvents [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_LanczosResolvents

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2024 ABINIT group (JLJ, BR, MC)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

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

m_hamiltonian/build_preconditioned_Hamiltonian_Lanczos_basis [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  build_preconditioned_Hamiltonian_Lanczos_basis

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

243 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(dpc), intent(in) :: seed_vector(npw_g)
249 
250 integer :: l
251 integer :: ierr
252 integer :: mpi_communicator
253 
254 ! *************************************************************************
255 
256 
257 
258 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 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 &                            LR_alpha, LR_beta, Hamiltonian_Qk)
266 
267 
268 call diagonalize_lanczos_banded(LR_kmax,LR_nseeds,npw_g, LR_alpha,LR_beta,  & 
269 Hamiltonian_Qk, LR_Hamiltonian_eigenvalues,.false.)
270 
271 
272 
273 ! Update the Lanczos vectors with the preconditioning matrix
274 do l = 1, LR_kmax
275 Hamiltonian_Qk(:,l) =  precondition_C(:)*Hamiltonian_Qk(:,l)
276 end do
277 
278 ! compute the M matrix
279 
280 ! Compute Q^dagger . Q
281 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 call xmpi_sum(LR_M_matrix,mpi_communicator,ierr) ! sum on all processors working on FFT!
283 
284 do l = 1, LR_kmax
285 LR_M_matrix(l,:) = LR_M_matrix(l,:)*LR_Hamiltonian_eigenvalues(l)
286 end do
287 
288 end subroutine build_preconditioned_Hamiltonian_Lanczos_basis

m_hamiltonian/cleanup_LanczosResolvents [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  cleanup_LanczosResolvents

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

156 subroutine cleanup_LanczosResolvents
157 !----------------------------------------------------------------------------------------------------
158 !
159 ! This subroutine cleans up global arrays.
160 !
161 !
162 !----------------------------------------------------------------------------------------------------
163 implicit none
164 
165 ! *************************************************************************
166 
167 ABI_FREE(precondition_C)
168 ABI_FREE(precondition_one_on_C)
169 ABI_FREE(LR_alpha)
170 ABI_FREE(LR_beta )
171 ABI_FREE(LR_seeds)
172 ABI_FREE(Hamiltonian_Qk)
173 ABI_FREE(LR_Hamiltonian_eigenvalues)
174 ABI_FREE(LR_M_matrix)
175 
176 end subroutine cleanup_LanczosResolvents

m_hamiltonian/compute_resolvent_column_shift_lanczos [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  compute_resolvent_column_shift_lanczos

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

305 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(dpc), intent(in) :: list_z(nz)
319 
320 integer     , intent(in) :: nvec
321 complex(dpc), intent(in) :: list_left_vectors(npw_g,nvec)
322 
323 complex(dpc), intent(in) :: seed_vector(npw_g)
324 
325 complex(dpc), intent(out) :: matrix_elements_resolvent(nz,nvec)
326 
327 ! local variables
328 integer        :: iz, l
329 complex(dpc)   :: z
330 integer        :: ierr
331 integer        :: mpi_communicator
332 
333 complex(dpc)   :: right_vec(LR_kmax)
334 complex(dpc)   :: left_vecs(LR_kmax,nvec)
335 complex(dpc)   :: work_vec(LR_kmax)
336 complex(dpc)   :: shift_lanczos_matrix(LR_kmax,LR_kmax)
337 
338 complex(dpc)   :: work_array(npw_g)
339 
340 ! *************************************************************************
341 
342 
343 ! processes will communicate along FFT rows
344 mpi_communicator = mpi_enreg%comm_fft
345 
346 !----------------------------------------------------------------------------------------------------
347 ! First, build the Lanczos basis for the preconditioned Hamiltonian
348 !----------------------------------------------------------------------------------------------------
349 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 work_array(:) = precondition_one_on_C(:)**2*seed_vector 
358 
359 call ZGEMV('C',npw_g,LR_kmax,cmplx_1,Hamiltonian_Qk,npw_g,work_array,1,cmplx_0,right_vec,1)
360 call xmpi_sum(right_vec,mpi_communicator,ierr) ! sum on all processors working on FFT!
361 
362 !  Q^dagger | left_vectors > 
363 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 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 do iz = 1, nz
371 
372 z = list_z(iz)
373 
374 ! Generate the matrix to be inverted
375 shift_lanczos_matrix(:,:) = LR_M_matrix(:,:)
376 
377 do l = 1, LR_kmax
378 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 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 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 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 end subroutine compute_resolvent_column_shift_lanczos

m_hamiltonian/compute_resolvent_column_shift_lanczos_right_vectors [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  compute_resolvent_column_shift_lanczos_right_vectors

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

413 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(dpc), intent(in) :: seed_vector(npw_g)
430 complex(dpc), intent(out):: right_vec(LR_kmax)
431 
432 ! local variables
433 integer        :: mpi_communicator
434 integer        :: ierr
435 
436 complex(dpc)   :: work_array(npw_g)
437 
438 ! *************************************************************************
439 
440 
441 ! processes will communicate along FFT rows
442 mpi_communicator = mpi_enreg%comm_fft
443 
444 !----------------------------------------------------------------------------------------------------
445 ! First, build the Lanczos basis for the preconditioned Hamiltonian
446 !----------------------------------------------------------------------------------------------------
447 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 work_array(:) = precondition_one_on_C(:)**2*seed_vector(:)
455 
456 call ZGEMV('C', npw_g, LR_kmax, cmplx_1, Hamiltonian_Qk, npw_g, work_array, 1, cmplx_0, right_vec, 1)
457 call xmpi_sum(right_vec,mpi_communicator,ierr) ! sum on all processors working on FFT!
458 
459 
460 end subroutine compute_resolvent_column_shift_lanczos_right_vectors

m_hamiltonian/invert_general_matrix [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  invert_general_matrix

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

478 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(dpc), intent(inout) :: matrix(n,n)
486 
487 integer :: info
488 integer :: ipiv(n)
489 
490 integer        :: debug_unit
491 character(50)  :: debug_filename
492 
493 
494 
495 
496 complex(dpc) :: work(n)
497 
498 ! *************************************************************************
499 
500 call ZGETRF( n, n, matrix, n, ipiv, info )
501 if ( info /= 0) then        
502   debug_unit = get_unit()
503   write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
504 
505   open(debug_unit,file=trim(debug_filename),status='unknown')
506 
507   write(debug_unit,'(A)')      '*********************************************************************************************'
508   write(debug_unit,'(A,I4,A)') '*      ERROR: info = ',info,' in ZGETRF(1), gwls_LanczosResolvents'
509   write(debug_unit,'(A)')      '*********************************************************************************************'
510 
511   close(debug_unit)
512 
513 end if
514 
515 
516 call ZGETRI( n, matrix, n, ipiv, work, n, info )
517 
518 if ( info /= 0) then        
519   debug_unit = get_unit()
520   write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
521 
522   open(debug_unit,file=trim(debug_filename),status='unknown')
523 
524   write(debug_unit,'(A)')      '*********************************************************************************************'
525   write(debug_unit,'(A,I4,A)') '*      ERROR: info = ',info,' in ZGETRI(1), gwls_LanczosResolvents'
526   write(debug_unit,'(A)')      '*********************************************************************************************'
527 
528   close(debug_unit)
529 
530 end if
531 
532 
533 end subroutine invert_general_matrix

m_hamiltonian/matrix_function_preconditioned_Hamiltonian [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  matrix_function_preconditioned_Hamiltonian

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

192 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(dpc), intent(out) :: vector_out(Hsize)
204 complex(dpc), intent(in)  :: vector_in(Hsize)
205 
206 ! local variables
207 real(dp)  :: psikg(2,npw_g)
208 
209 complex(dpc)  :: tmp_vector(Hsize)
210 
211 ! *************************************************************************
212 
213 ! convert from one format to the other
214 
215 tmp_vector(:) = precondition_one_on_C(:)*vector_in(:)
216 
217 psikg(1,:) = dble (tmp_vector(:))
218 psikg(2,:) = dimag(tmp_vector(:))
219 
220 call Hpsik(psikg)
221 
222 tmp_vector(:) = cmplx_1*psikg(1,:)+cmplx_i*psikg(2,:)
223 
224 vector_out(:) = precondition_one_on_C(:)*tmp_vector(:)
225 
226 
227 end subroutine matrix_function_preconditioned_Hamiltonian

m_hamiltonian/setup_LanczosResolvents [ Functions ]

[ Top ] [ m_hamiltonian ] [ Functions ]

NAME

  setup_LanczosResolvents

FUNCTION

  .

INPUTS

OUTPUT

SOURCE

 92 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 LR_kmax   = kmax
113 LR_nseeds = 1      ! only one seed
114 
115 ! Prepare the algorithm
116 
117 ABI_MALLOC(precondition_C,        (npw_g))
118 ABI_MALLOC(precondition_one_on_C, (npw_g))
119 
120 
121 if ( prec ) then
122   call set_precondition()
123 
124   precondition_one_on_C(:) = sqrt(pcon(:))
125   !precondition_one_on_C(:) = pcon(:) ! yields very bad results!
126   precondition_C(:)        = cmplx_1/precondition_one_on_C(:)
127 
128 else 
129   precondition_C(:)        = cmplx_1
130   precondition_one_on_C(:) = cmplx_1
131 end if
132 
133 ABI_MALLOC(LR_alpha,       (LR_nseeds,LR_nseeds,LR_kmax))
134 ABI_MALLOC(LR_beta ,       (LR_nseeds,LR_nseeds,LR_kmax))
135 ABI_MALLOC(LR_seeds,       (npw_g,LR_nseeds))
136 ABI_MALLOC(Hamiltonian_Qk, (npw_g,LR_kmax))
137 ABI_MALLOC(LR_Hamiltonian_eigenvalues, (LR_kmax))
138 ABI_MALLOC(LR_M_matrix, (LR_kmax,LR_kmax))
139 
140 end subroutine setup_LanczosResolvents