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