Branch data Line data Source code
1 : : !!****m* ABINIT/m_gwls_GWlanczos
2 : : !! NAME
3 : : !! m_gwls_GWlanczos
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 : : module m_gwls_GWlanczos
23 : : !----------------------------------------------------------------------------------------------------
24 : : ! This module implements the Lanczos scheme to band diagonalize an implicit operator.
25 : : !----------------------------------------------------------------------------------------------------
26 : : !local modules
27 : : use m_gwls_utility
28 : : use m_gwls_TimingLog
29 : : use m_gwls_wf
30 : : use m_gwls_hamiltonian
31 : : use m_gwls_lineqsolver
32 : : use m_gwls_polarisability
33 : : use m_gwls_QR_factorization
34 : :
35 : : !abinit modules
36 : : use defs_basis
37 : : use defs_wvltypes
38 : : use m_abicore
39 : : use m_xmpi
40 : : use m_pawang
41 : : use m_errors
42 : :
43 : : use m_io_tools, only : get_unit
44 : :
45 : : implicit none
46 : : save
47 : : private
48 : : !!***
49 : :
50 : : !!***
51 : : public :: block_lanczos_algorithm
52 : :
53 : : public :: diagonalize_lanczos_banded
54 : : public :: get_seeds
55 : : !!***
56 : : contains
57 : :
58 : :
59 : : !!****f* m_gwls_GWlanczos/get_seeds
60 : : !! NAME
61 : : !! get_seeds
62 : : !!
63 : : !! FUNCTION
64 : : !! .
65 : : !!
66 : : !! INPUTS
67 : : !!
68 : : !! OUTPUT
69 : : !!
70 : : !! SOURCE
71 : :
72 : 29 : subroutine get_seeds(first_seed, nseeds, seeds)
73 : : !----------------------------------------------------------------------------------------------------
74 : : !
75 : : ! This subroutine compute the seeds using the eigenstates of the Hamiltonian
76 : : !
77 : : !----------------------------------------------------------------------------------------------------
78 : : integer, intent(in) :: first_seed, nseeds
79 : : complex(dp), intent(out) :: seeds(npw_k,nseeds)
80 : :
81 : 29 : real(dp) , allocatable :: psik_out(:,:)
82 : 29 : real(dp) , allocatable :: psikb_e(:,:)
83 : 29 : real(dp) , allocatable :: psig_e(:,:)
84 : 29 : real(dp) , allocatable :: psikb_s(:,:)
85 : 29 : real(dp) , allocatable :: psig_s(:,:)
86 : :
87 : : ! local variables
88 : : integer :: n
89 : : integer :: i, j, nsblk
90 : : ! *************************************************************************
91 : :
92 : : ! Generate the seeds for the Lanczos algorithm
93 [ - + + - : 87 : ABI_MALLOC(psik_out,(2,npw_k))
+ - + - -
+ - + ]
94 [ - + + - : 87 : ABI_MALLOC(psikb_e,(2,npw_kb))
+ - + - -
+ - + ]
95 [ - + + - : 87 : ABI_MALLOC(psig_e,(2,npw_g))
+ - + - -
+ - + ]
96 [ - + + - : 58 : ABI_MALLOC(psikb_s,(2,npw_kb))
- + - + ]
97 [ - + + - : 58 : ABI_MALLOC(psig_s,(2,npw_g))
- + - + ]
98 : :
99 [ + + ]: 29 : nsblk = ceiling(1.0*nseeds/blocksize)
100 : :
101 [ + + ]: 68 : do i=1,nsblk
102 [ + + ]: 98 : do j=1,blocksize
103 [ + + + + ]: 22457 : psikb_e(:,(j-1)*npw_k+1:j*npw_k) = cg(:,(e-1)*npw_k+1:e*npw_k)
104 : : end do
105 : :
106 [ + + + + ]: 22398 : psig_e = zero
107 : 39 : call wf_block_distribute(psikb_e, psig_e,1) ! LA -> FFT
108 : :
109 [ + + ]: 98 : do j=1,blocksize
110 : 59 : n = (i-1)*blocksize + j-1 + first_seed
111 [ + + ]: 98 : if ((i-1)*blocksize + j <= nseeds) then
112 [ + + + + ]: 19318 : psikb_s(:,(j-1)*npw_k+1:j*npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k)
113 : : else
114 [ + + + + ]: 3100 : psikb_s(:,(j-1)*npw_k+1:j*npw_k) = zero
115 : : end if
116 : : end do
117 : :
118 [ + + + + ]: 22398 : psig_s = zero
119 : 39 : call wf_block_distribute(psikb_s, psig_s,1) ! LA -> FFT
120 : :
121 : : ! Fourier transform valence wavefunction, to real space
122 : 39 : call g_to_r(psir1,psig_s)
123 : :
124 [ + + + + : 267501 : psir1(2,:,:,:) = -psir1(2,:,:,:)
+ + ]
125 : :
126 : 39 : call gr_to_g(psig_s,psir1,psig_e)
127 : :
128 : : ! return to LA configuration, in order to apply Coulomb potential
129 : 39 : call wf_block_distribute(psikb_s, psig_s,2) ! FFT -> LA
130 : :
131 [ + + ]: 127 : do j=1, blocksize
132 : 59 : n = (i-1)*blocksize + j
133 [ + + ]: 98 : if(n<=nseeds) then
134 [ + - - + : 19361 : psik_out = psikb_s(:,(j-1)*npw_k+1:j*npw_k)
- - + + +
+ ]
135 : 43 : call sqrt_vc_k(psik_out)
136 [ + + ]: 6468 : seeds(:,n) = cmplx_1*psik_out(1,:) + cmplx_i*psik_out(2,:)
137 : : end if
138 : : end do
139 : : end do
140 : :
141 [ - + ]: 29 : ABI_FREE(psik_out)
142 [ - + ]: 29 : ABI_FREE(psikb_e)
143 [ - + ]: 29 : ABI_FREE(psig_e)
144 [ - + ]: 29 : ABI_FREE(psikb_s)
145 [ - + ]: 29 : ABI_FREE(psig_s)
146 : :
147 : 29 : end subroutine get_seeds
148 : : !!***
149 : :
150 : : !!****f* m_gwls_GWlanczos/block_lanczos_algorithm
151 : : !! NAME
152 : : !! block_lanczos_algorithm
153 : : !!
154 : : !! FUNCTION
155 : : !! .
156 : : !!
157 : : !! INPUTS
158 : : !!
159 : : !! OUTPUT
160 : : !!
161 : : !! SOURCE
162 : :
163 [ - + - - : 173 : subroutine block_lanczos_algorithm(mpi_communicator,matrix_function,kmax,nseeds,Hsize,seeds,alpha,beta,Lbasis,X0,beta0,Qk)
- - ]
164 : : !----------------------------------------------------------------------------------------------------
165 : : !
166 : : ! This subroutine implements the Block Lanczos algorithm for an arbitrary, user-supplied function which
167 : : ! returns the action of the implicit matrix, which is assumed to be Hermitian.
168 : : !
169 : : !
170 : : !----------------------------------------------------------------------------------------------------
171 : :
172 : : !-----------------------------------------
173 : : ! interface with implicit matrix function
174 : : !-----------------------------------------
175 : : !********************************************************
176 : : !*** NOTE: ***
177 : : !*** ***
178 : : !*** There *appears* to be a bug in makemake; ***
179 : : !*** when the name of the function in the interface ***
180 : : !*** below contains the word "implicit", makemake ***
181 : : !*** yields an error. I suppose makemake simply ***
182 : : !*** parses for the word "implicit" without regards ***
183 : : !*** for the fact that it may simply be part of a ***
184 : : !*** naive developper's chosen name for the routine. ***
185 : : !*** Correspondingly, I change the name to ***
186 : : !*** "matrix_function" to avoid problems. ***
187 : : !*** ***
188 : : !*** Bruno Rousseau ***
189 : : !*** 08/13/2012 ***
190 : : !********************************************************
191 : : interface
192 : : subroutine matrix_function(v_out,v_in,l)
193 : :
194 : : use defs_basis
195 : :
196 : : integer, intent(in) :: l
197 : : complex(dp), intent(out) :: v_out(l)
198 : : complex(dp), intent(in) :: v_in(l)
199 : :
200 : : end subroutine matrix_function
201 : : end interface
202 : :
203 : : !------------------------------
204 : : ! input/output variables
205 : : !------------------------------
206 : :
207 : : integer, intent(in) :: mpi_communicator
208 : : integer, intent(in) :: kmax ! number of Lanczos blocks
209 : : integer, intent(in) :: nseeds ! size of each blocks
210 : : integer, intent(in) :: Hsize ! size of the Hilbert space in which the matrix lives
211 : :
212 : : complex(dp), intent(inout):: seeds(Hsize,nseeds) ! seed vectors for the algorithm
213 : : ! overwritten by X_{k+1} on output
214 : :
215 : : !logical, intent(in) :: ortho ! should the Lanczos vector be orthogonalized?
216 : :
217 : : complex(dp), intent(out) :: alpha(nseeds,nseeds,kmax) ! the alpha array from the Lanczos algorithm
218 : : complex(dp), intent(out) :: beta(nseeds,nseeds,kmax) ! the beta array from the Lanczos algorithm
219 : : complex(dp), intent(out) :: Lbasis(Hsize,nseeds*kmax) ! array containing the Lanczos basis
220 : :
221 : :
222 : : complex(dp), intent(in),optional :: X0(Hsize,nseeds)
223 : : complex(dp), intent(in),optional :: beta0(nseeds,nseeds)
224 : : complex(dp), intent(in),optional :: Qk(:,:) ! array containing vectors to which
225 : :
226 : : ! the basis must be orthonormalized
227 : :
228 : :
229 : :
230 : : !------------------------------
231 : : ! local variables
232 : : !------------------------------
233 : : integer :: k, seed1
234 : : integer :: dum(2), lk
235 : :
236 : 173 : complex(dp), allocatable :: xk(:,:), xkm1(:,:), rk(:,:)
237 : :
238 : : integer :: ntime, itime
239 : : real(dp) :: total_time1, total_time2
240 : : real(dp) :: time1, time2
241 : : integer :: ierr
242 : 173 : real(dp),allocatable :: list_time(:)
243 : :
244 : : ! *************************************************************************
245 : :
246 : 173 : call cpu_time(total_time1)
247 : :
248 : :
249 : 173 : ntime = 7
250 [ + - ]: 173 : ABI_MALLOC(list_time,(ntime))
251 [ + + ]: 1384 : list_time(:) = zero
252 : :
253 [ - + ]: 173 : if(present(Qk)) then
254 [ # # ]: 0 : dum = shape(Qk)
255 : 0 : lk = dum(2)
256 : : end if
257 : :
258 [ - + + - : 692 : ABI_MALLOC( xk, (Hsize,nseeds))
- + + - +
- + - - +
- + ]
259 [ - + + - : 519 : ABI_MALLOC( xkm1,(Hsize,nseeds))
- + + - -
+ - + ]
260 [ - + + - : 519 : ABI_MALLOC( rk ,(Hsize,nseeds))
- + + - -
+ + - ]
261 : :
262 : :
263 [ + + + + : 6165 : alpha = cmplx_0
+ + ]
264 [ + + + + : 6165 : beta = cmplx_0
+ + ]
265 : :
266 : : !------------------------------------------------
267 : : ! Orthonormalize the seeds
268 : : !------------------------------------------------
269 : : ! initialize the xk array with the seeds
270 [ + + + + ]: 35569 : xk(:,:) = seeds(:,:)
271 : :
272 : :
273 : : ! orthonormalize the block using the QR algorithm
274 : : ! xk is overwritten by Q, the array of orthonormal vectors
275 : :
276 : 173 : call extract_QR(mpi_communicator, Hsize,nseeds,xk)
277 : :
278 : : !------------------------------------------------
279 : : ! Loop on all blocks
280 : : !------------------------------------------------
281 : :
282 : :
283 [ + + ]: 2077 : do k = 1, kmax
284 : :
285 : 1904 : itime = 0
286 : :
287 : : ! tabulate basis, computed at previous step
288 [ + + + + ]: 384224 : Lbasis(:,nseeds*(k-1)+1:nseeds*k) = xk(:,:)
289 : :
290 : :
291 : 1904 : itime = itime+1
292 : 1904 : call cpu_time(time1)
293 : :
294 : : ! Initialize the residual array
295 [ + + ]: 3864 : do seed1 = 1, nseeds
296 : : ! If we are constructing the $\hat \epsilon(i\omega = 0)$ matrix (and the Lanczos basis at the same time),
297 : : ! note the index in which the Sternheimer solutions will be stored (for use in the projected Sternheimer section).
298 [ + + ]: 1960 : if(write_solution) index_solution = (k-1)*nseeds + seed1
299 : :
300 : 3864 : call matrix_function(rk(:,seed1),xk(:,seed1),Hsize)
301 : : end do
302 : :
303 : 1904 : call cpu_time(time2)
304 : 1904 : list_time(itime) = list_time(itime) + time2-time1
305 : :
306 : 1904 : itime = itime+1
307 : 1904 : call cpu_time(time1)
308 : : ! compute the alpha array, alpha = X^d.A.X
309 : :
310 : : call ZGEMM( 'C', & ! take Hermitian conjugate of first array
311 : : 'N', & ! leave second array as is
312 : : nseeds, & ! the number of rows of the matrix op( A )
313 : : nseeds, & ! the number of columns of the matrix op( B )
314 : : Hsize, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
315 : : cmplx_1, & ! alpha constant
316 : : xk, & ! matrix A
317 : : Hsize, & ! LDA
318 : : rk, & ! matrix B
319 : : Hsize, & ! LDB
320 : : cmplx_0, & ! beta constant
321 : : alpha(:,:,k), & ! matrix C
322 : 1904 : nseeds) ! LDC
323 : 1904 : call xmpi_sum(alpha(:,:,k),mpi_communicator,ierr) ! sum on all processors
324 : :
325 : :
326 : :
327 : 1904 : call cpu_time(time2)
328 : 1904 : list_time(itime) = list_time(itime) + time2-time1
329 : :
330 : 1904 : itime = itime+1
331 : 1904 : call cpu_time(time1)
332 : : ! update the residual array, rk = rk-X.alpha
333 : : call ZGEMM( 'N', & ! leave first array as is
334 : : 'N', & ! leave second array as is
335 : : Hsize, & ! the number of rows of the matrix op( A )
336 : : nseeds, & ! the number of columns of the matrix op( B )
337 : : nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
338 : : -cmplx_1, & ! alpha constant
339 : : xk, & ! matrix A
340 : : Hsize, & ! LDA
341 : : alpha(:,:,k), & ! matrix B
342 : : nseeds, & ! LDB
343 : : cmplx_1, & ! beta constant
344 : : rk, & ! matrix C
345 : 1904 : Hsize) ! LDC
346 : :
347 : 1904 : call cpu_time(time2)
348 : 1904 : list_time(itime) = list_time(itime) + time2-time1
349 : :
350 [ - + - - ]: 1904 : if (k .eq. 1 .and. present(X0) .and. present(beta0)) then
351 : : ! if k == 1, and X0,beta0 are present,
352 : : ! update the residual array, r1 = r1-X_{0}.beta^d_{0}
353 : : call ZGEMM( 'N', & ! leave first array as is
354 : : 'C', & ! Hermitian conjugate the second array
355 : : Hsize, & ! the number of rows of the matrix op( A )
356 : : nseeds, & ! the number of columns of the matrix op( B )
357 : : nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
358 : : -cmplx_1, & ! alpha constant
359 : : X0, & ! matrix A
360 : : Hsize, & ! LDA
361 : : beta0(:,:), & ! matrix B
362 : : nseeds, & ! LDB
363 : : cmplx_1, & ! beta constant
364 : : rk, & ! matrix C
365 : 0 : Hsize) ! LDC
366 : : end if
367 : :
368 : :
369 : 1904 : itime = itime+1
370 : 1904 : call cpu_time(time1)
371 [ + + ]: 1904 : if (k .gt. 1) then
372 : :
373 : : ! if k > 1, update the residual array, rk = rk-X_{k-1}.beta^d_{k-1}
374 : : call ZGEMM( 'N', & ! leave first array as is
375 : : 'C', & ! Hermitian conjugate the second array
376 : : Hsize, & ! the number of rows of the matrix op( A )
377 : : nseeds, & ! the number of columns of the matrix op( B )
378 : : nseeds, & ! the number of columns of the matrix op( A ) == rows of matrix op( B )
379 : : -cmplx_1, & ! alpha constant
380 : : xkm1, & ! matrix A
381 : : Hsize, & ! LDA
382 : : beta(:,:,k-1), & ! matrix B
383 : : nseeds, & ! LDB
384 : : cmplx_1, & ! beta constant
385 : : rk, & ! matrix C
386 : 1731 : Hsize) ! LDC
387 : :
388 : : end if
389 : 1904 : call cpu_time(time2)
390 : 1904 : list_time(itime) = list_time(itime) + time2-time1
391 : :
392 : : ! store xk for next iteration
393 [ + + + + ]: 384224 : xkm1(:,:) = xk(:,:)
394 : :
395 : :
396 : : ! Orthonormalize THE RESIDUAL to all previously calculated directions
397 : :
398 : 1904 : itime = itime+1
399 : 1904 : call cpu_time(time1)
400 : : !if ( ortho .and. (dtset%gwcalctyp/=1) ) then !This is a test to obtain the CPU time taken by the orthogonalizations.
401 : :
402 [ - + ]: 1904 : if(present(Qk)) then
403 : : ! Orthonormalize to all previously calculated directions, if
404 : : ! this is a restarted Lanczos step
405 [ # # # # : 0 : call orthogonalize(mpi_communicator, Hsize,lk,nseeds,Qk,rk)
# # # # #
# ]
406 : : end if
407 : :
408 : 1904 : call orthogonalize(mpi_communicator, Hsize,k*nseeds,nseeds,Lbasis(:,1:k*nseeds),rk)
409 : :
410 : : !end if
411 : 1904 : call cpu_time(time2)
412 : 1904 : list_time(itime) = list_time(itime) + time2-time1
413 : :
414 : :
415 : 1904 : itime = itime+1
416 : 1904 : call cpu_time(time1)
417 : :
418 : : ! perform QR decomposition to extract X_{k+1} and beta_{k}
419 : 1904 : call extract_QR(mpi_communicator, Hsize,nseeds,rk,beta(:,:,k))
420 : :
421 : 1904 : call cpu_time(time2)
422 : 1904 : list_time(itime) = list_time(itime) + time2-time1
423 : : ! copy the Q matrix (written on rk) in xk, which becomes X_{k+1}
424 [ + + + + ]: 384397 : xk(:,:) = rk(:,:)
425 : :
426 : :
427 : : end do !end loop on k
428 : :
429 : : ! overwrite the seeds with the last vector block.
430 [ + + + + ]: 35569 : seeds(:,:) = xk(:,:)
431 : :
432 [ - + ]: 173 : ABI_FREE( xk )
433 [ - + ]: 173 : ABI_FREE( xkm1)
434 [ - + ]: 173 : ABI_FREE( rk )
435 : 173 : call cpu_time(total_time2)
436 : :
437 : 173 : list_time(7) = total_time2-total_time1
438 : :
439 : 173 : call write_block_lanczos_timing_log(list_time,ntime)
440 : :
441 : 173 : ABI_FREE(list_time)
442 : :
443 : 173 : end subroutine block_lanczos_algorithm
444 : : !!***
445 : :
446 : : !!****f* m_gwls_GWlanczos/diagonalize_lanczos_banded
447 : : !! NAME
448 : : !! diagonalize_lanczos_banded
449 : : !!
450 : : !! FUNCTION
451 : : !! .
452 : : !!
453 : : !! INPUTS
454 : : !!
455 : : !! OUTPUT
456 : : !!
457 : : !! SOURCE
458 : :
459 : 177 : subroutine diagonalize_lanczos_banded(kmax,nseeds,Hsize,alpha,beta,Lbasis,eigenvalues,debug)
460 : : !-----------------------------------------------------------------------------------
461 : : ! Given the result of the Lanczos algorithm, this subroutine diagonalize the banded
462 : : ! matrix as well as updates the basis.
463 : : !-----------------------------------------------------------------------------------
464 : : integer, intent(in) :: kmax ! number of Lanczos blocks
465 : : integer, intent(in) :: nseeds ! size of each blocks
466 : : integer, intent(in) :: Hsize ! size of the Hilbert space in which the matrix lives
467 : : logical, intent(in) :: debug
468 : :
469 : : complex(dp), intent(in) :: alpha(nseeds,nseeds,kmax) ! the alpha array from the Lanczos algorithm
470 : : complex(dp), intent(in) :: beta (nseeds,nseeds,kmax) ! the beta array from the Lanczos algorithm
471 : :
472 : : complex(dp), intent(inout) :: Lbasis(Hsize,nseeds*kmax) ! array containing the Lanczos basis
473 : :
474 : :
475 : : real(dp), intent(out) :: eigenvalues(nseeds*kmax)
476 : :
477 : :
478 : : ! local variables
479 : :
480 : : integer :: kd ! number of superdiagonal above the diagonal in banded storage
481 : : integer :: ldab ! dimension of banded storage matrix
482 : :
483 : 177 : complex(dp), allocatable :: band_storage_matrix(:,:)
484 : 177 : complex(dp), allocatable :: saved_band_storage_matrix(:,:)
485 : :
486 : 177 : complex(dp), allocatable :: eigenvectors(:,:)
487 : :
488 : 177 : complex(dp), allocatable :: Lbasis_tmp(:,:)
489 : :
490 : : integer :: i, j
491 : : integer :: k
492 : : integer :: s1, s2
493 : : integer :: info
494 : :
495 : :
496 : 177 : complex(dp), allocatable :: work(:)
497 : 177 : real(dp), allocatable :: rwork(:)
498 : :
499 : : integer :: io_unit
500 : : character(128) :: filename
501 : : logical :: file_exists
502 : :
503 : : integer :: debug_unit
504 : : character(50) :: debug_filename
505 : :
506 : : ! *************************************************************************
507 : :
508 : :
509 : :
510 : : ! number of superdiagonals
511 : 177 : kd = nseeds
512 : 177 : ldab = kd + 1
513 : :
514 [ - + + - : 708 : ABI_MALLOC( band_storage_matrix, (ldab,nseeds*kmax))
- + + - +
- + - - +
- + ]
515 [ - + + - : 531 : ABI_MALLOC(saved_band_storage_matrix, (ldab,nseeds*kmax))
- + + - -
+ - + ]
516 : : !---------------------------------------------------------
517 : : ! Store banded matrix in banded format
518 : : !---------------------------------------------------------
519 : : ! for UPLO = 'L', AB(1+i-j,j) = A(i,j) for j<=i<=min(n,j+kd).
520 : :
521 [ + + + + ]: 6297 : band_storage_matrix(:,:) = cmplx_0
522 : :
523 : : !-----------------------------------------
524 : : ! Store the alpha and beta matrices
525 : : !-----------------------------------------
526 : :
527 : : ! loop on all blocks
528 [ - + ]: 1928 : do k=1,kmax
529 : :
530 : : ! alpha blocks
531 [ + + ]: 3912 : do s2 = 1, nseeds
532 : 1984 : j = (k-1)*nseeds+s2
533 : :
534 [ + + ]: 5980 : do s1 = s2, nseeds
535 : 2068 : i = j+s1-s2
536 : 4052 : band_storage_matrix(1+i-j,j) = alpha(s1,s2,k)
537 : : end do
538 : : end do
539 : :
540 : : ! exit when k = kmax, as this beta block does not contribute.
541 [ + + ]: 1928 : if (k .eq. kmax) exit
542 : :
543 : : ! beta blocks
544 [ + + ]: 3721 : do s2 = 1, nseeds
545 : 1793 : j = (k-1)*nseeds+s2
546 : :
547 [ + + ]: 5400 : do s1 = 1, s2
548 : 1856 : i = j+s1-s2+nseeds
549 : 3649 : band_storage_matrix(1+i-j,j) = beta(s1,s2,k)
550 : : end do
551 : :
552 : : end do
553 : : end do
554 : :
555 [ + + + + ]: 6297 : saved_band_storage_matrix(:,:) = band_storage_matrix(:,:)
556 : :
557 : : !-----------------------------------------
558 : : ! Diagonalize the banded matrix
559 : : !-----------------------------------------
560 : :
561 [ - + + - : 708 : ABI_MALLOC(eigenvectors, (nseeds*kmax,nseeds*kmax))
- + + - +
- + - - +
- + ]
562 : :
563 [ - + + - : 531 : ABI_MALLOC(work,(nseeds*kmax))
+ - + - -
+ - + ]
564 [ + - - + : 531 : ABI_MALLOC(rwork,(3*nseeds*kmax-2))
+ - + - +
- - + -
+ ]
565 : :
566 : : call ZHBEV( 'V', & ! compute eigenvalues and eigenvectors
567 : : 'L', & ! lower triangular part of matrix is stored in banded_matrix
568 : : nseeds*kmax, & ! dimension of matrix
569 : : kd, & ! number of superdiagonals in banded matrix
570 : : band_storage_matrix, & ! matrix in banded storage
571 : : ldab, & ! leading dimension of banded_matrix
572 : : eigenvalues, & ! eigenvalues of matrix
573 : : eigenvectors, & ! eigenvectors of matrix
574 : : nseeds*kmax, & ! dimension of eigenvector matrix
575 : 177 : work, rwork, info ) ! work arrays and info
576 : :
577 : :
578 [ - + ]: 177 : if ( info /= 0) then
579 : 0 : debug_unit = get_unit()
580 : 0 : write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
581 : :
582 : 0 : open(debug_unit,file=trim(debug_filename),status='unknown')
583 : :
584 : 0 : write(debug_unit,'(A)') '*********************************************************************************************'
585 : 0 : write(debug_unit,'(A,I4,A)') '* ERROR: info = ',info,' in ZHBEV (1), gwls_GWlanczos'
586 : 0 : write(debug_unit,'(A)') '*********************************************************************************************'
587 : :
588 : 0 : close(debug_unit)
589 : :
590 : : end if
591 : :
592 : :
593 : :
594 : : !----------------------------------------------------------------------------------
595 : : ! update the Lanczos basis to reflect diagonalization of T matrix
596 : : !
597 : : ! Note that by definition
598 : : !
599 : : ! Q^H . A . Q = T ==> A = Q . T . Q^H
600 : : !
601 : : ! where Q (Lbasis) contains the Lanczos basis.
602 : : !
603 : : ! Diagonalizing T, ie T = U . LAMBDA . U^H leads to
604 : : !
605 : : ! A = [ Q.U] . LAMBDA . [Q.U]^H
606 : : !
607 : : ! The updated basis is thus Q.U == Lbasis . eigenvectors
608 : : !----------------------------------------------------------------------------------
609 : :
610 : : ! NEVER use matmul!!! It sends temporary arrays to the stack, which can be much smaller
611 : : ! than needed; this leads to mysterious segfaults!
612 : :
613 : : ! Lbasis = matmul(Lbasis,eigenvectors)
614 : :
615 : : ! use temporary array, which is PROPERLY ALLOCATED, to perform matrix multiplication
616 [ - + + - : 708 : ABI_MALLOC(Lbasis_tmp, (Hsize,nseeds*kmax))
- + + - +
- + - - +
- + ]
617 : :
618 : : ! Compute C = A * B, where A = Lbasis, B = eigenvectors, and C = Lbasis_tmp
619 : : call ZGEMM( 'N', & ! leave array A as is
620 : : 'N', & ! leave array B as is
621 : : Hsize, & ! number of rows of A
622 : : nseeds*kmax, & ! number of columns of B
623 : : nseeds*kmax, & ! number of columns of A /rows of B
624 : : cmplx_1, & ! constant alpha
625 : : Lbasis, & ! matrix A
626 : : Hsize, & ! LDA
627 : : eigenvectors, & ! matrix B
628 : : nseeds*kmax, & ! LDB
629 : : cmplx_0, & ! constant beta
630 : : Lbasis_tmp, & ! matrix C
631 : 177 : Hsize) ! LDC
632 : :
633 : : ! overwrite initial array
634 [ + + + + ]: 388689 : Lbasis(:,:) = Lbasis_tmp(:,:)
635 : :
636 : :
637 [ - + ]: 177 : ABI_FREE(Lbasis_tmp)
638 : :
639 [ - + - - ]: 177 : if ( debug .and. mpi_enreg%me == 0) then
640 : : !----------------------------------------------------------------------------------
641 : : ! For the purpose of debugging, print relevant results to file to check
642 : : ! data consistency. This may not be necessary once the code has been shown to
643 : : ! work properly.
644 : : !----------------------------------------------------------------------------------
645 : :
646 : 0 : io_unit = get_unit()
647 : 0 : i = 0
648 : 0 : file_exists = .true.
649 [ # # ]: 0 : do while (file_exists)
650 : 0 : i = i+1
651 : 0 : write(filename,'(A,I0.4,A)') "diagonalize_banded_matrix_",i,".log"
652 : 0 : inquire(file=filename,exist=file_exists)
653 : : end do
654 : :
655 : :
656 : 0 : open(io_unit,file=filename,status=files_status_new)
657 : 0 : write(io_unit,10) "#======================================================================================="
658 : 0 : write(io_unit,10) "# "
659 : 0 : write(io_unit,10) "# This file contains information pertaining to the diagonalization of a banded "
660 : 0 : write(io_unit,10) "# matrix, expressed in terms of the Lanczos alpha and beta block arrays. "
661 : 0 : write(io_unit,10) "# "
662 : 0 : write(io_unit,10) "#======================================================================================="
663 : 0 : write(io_unit,10) "# "
664 : 0 : write(io_unit,12) "# diagonalization info : ",info," "
665 : 0 : write(io_unit,10) "# "
666 : 0 : write(io_unit,10) "# l lambda_l "
667 : 0 : write(io_unit,10) "#======================================================================================="
668 : :
669 [ # # ]: 0 : do i = 1, nseeds*kmax
670 : 0 : write(io_unit,13) i, eigenvalues(i)
671 : : end do
672 : :
673 : 0 : write(io_unit,10) " "
674 : 0 : write(io_unit,10) "# "
675 : 0 : write(io_unit,10) "# alpha and beta blocks "
676 : 0 : write(io_unit,10) "# "
677 : 0 : write(io_unit,10) "#======================================================================================="
678 : :
679 : : ! loop on all blocks
680 [ # # ]: 0 : do k=1,kmax
681 : 0 : write(io_unit,10) "# "
682 : 0 : write(io_unit,12) "# block k = ",k," "
683 : 0 : write(io_unit,10) "# "
684 : 0 : write(io_unit,15) "# alpha: ||alpha^H-alpha|| = ", &
685 [ # # # # ]: 0 : sqrt(sum(abs(alpha(:,:,k)-transpose(conjg(alpha(:,:,k))))**2))
686 [ # # ]: 0 : do s1 = 1, nseeds
687 : 0 : write(io_unit,14) alpha(s1,:,k)
688 : : end do
689 : :
690 : 0 : write(io_unit,10) "# "
691 : 0 : write(io_unit,10) "# beta "
692 [ # # ]: 0 : do s1 = 1, nseeds
693 : 0 : write(io_unit,14) beta(s1,:,k)
694 : : end do
695 : :
696 : :
697 : : end do
698 : :
699 : :
700 : :
701 : :
702 : 0 : write(io_unit,10) "# "
703 : 0 : write(io_unit,10) "# band storage matrix: "
704 : 0 : write(io_unit,10) "#======================================================================================="
705 : :
706 [ # # ]: 0 : do i = 1, ldab
707 : 0 : write(io_unit,14) saved_band_storage_matrix(i,:)
708 : : end do
709 : :
710 : 0 : close(io_unit)
711 : : end if
712 : :
713 : :
714 : : ! clean up memory
715 [ - + ]: 177 : ABI_FREE(eigenvectors)
716 [ - + ]: 177 : ABI_FREE( work)
717 [ - + ]: 177 : ABI_FREE(rwork)
718 [ - + ]: 177 : ABI_FREE(band_storage_matrix)
719 [ - + ]: 177 : ABI_FREE(saved_band_storage_matrix)
720 : :
721 : :
722 : : 10 format(A)
723 : : 12 format(A,I5,A)
724 : : 13 format(I5,5X,F24.12)
725 : : 14 format(4X,1000(F12.8,SP,F12.8,1X,'i',2X))
726 : : 15 format(A,ES24.8)
727 : :
728 : 177 : end subroutine diagonalize_lanczos_banded
729 : : !!***
730 : :
731 : : end module m_gwls_GWlanczos
732 : : !!***
|
