TABLE OF CONTENTS


ABINIT/m_gwls_GWlanczos [ Modules ]

[ Top ] [ Modules ]

NAME

 m_gwls_GWlanczos

FUNCTION

  .

COPYRIGHT

 Copyright (C) 2009-2018 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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 
28 module m_gwls_GWlanczos
29 !----------------------------------------------------------------------------------------------------
30 ! This module implements the Lanczos scheme to band diagonalize an implicit operator.
31 !----------------------------------------------------------------------------------------------------
32 !local modules
33 use m_gwls_utility
34 use m_gwls_TimingLog
35 use m_gwls_wf
36 use m_gwls_hamiltonian
37 use m_gwls_lineqsolver
38 use m_gwls_polarisability
39 use m_gwls_QR_factorization
40 
41 !abinit modules
42 use defs_basis
43 use defs_datatypes
44 use defs_abitypes
45 use defs_wvltypes
46 use m_abicore
47 use m_xmpi
48 use m_pawang
49 use m_errors
50 
51 use m_io_tools,         only : get_unit
52 use m_paw_dmft,         only : paw_dmft_type
53 use m_ebands,           only : ebands_init, ebands_free
54 
55 implicit none
56 save
57 private

m_gwls_GWlanczos/block_lanczos_algorithm [ Functions ]

[ Top ] [ m_gwls_GWlanczos ] [ Functions ]

NAME

  block_lanczos_algorithm

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputePoles,gwls_GenerateEpsilon,gwls_LanczosResolvents

CHILDREN

      zgemm,zhbev

SOURCE

197 subroutine block_lanczos_algorithm(mpi_communicator,matrix_function,kmax,nseeds,Hsize,seeds,alpha,beta,Lbasis,X0,beta0,Qk)
198 !----------------------------------------------------------------------------------------------------
199 !
200 ! This subroutine implements the Block Lanczos algorithm for an arbitrary, user-supplied function which
201 ! returns the action of the implicit matrix, which is assumed to be Hermitian.
202 !
203 !
204 !----------------------------------------------------------------------------------------------------
205 
206 !This section has been created automatically by the script Abilint (TD).
207 !Do not modify the following lines by hand.
208 #undef ABI_FUNC
209 #define ABI_FUNC 'block_lanczos_algorithm'
210 !End of the abilint section
211 
212 implicit none
213 
214 !-----------------------------------------
215 ! interface with implicit matrix function
216 !-----------------------------------------
217 !********************************************************
218 !***        NOTE:                                           ***
219 !***                                                        ***
220 !***        There *appears* to be a bug in makemake;        ***
221 !***        when the name of the function in the interface  ***
222 !***        below contains the word "implicit", makemake    ***
223 !***        yields an error. I suppose makemake simply      ***
224 !***        parses for the word "implicit" without regards  ***
225 !***        for the fact that it may simply be part of a    ***
226 !***        naive developper's chosen name for the routine. ***
227 !***        Correspondingly, I change the name to           ***
228 !***        "matrix_function" to avoid problems.            ***
229 !***                                                        ***
230 !***                                     Bruno Rousseau     ***
231 !***                                     08/13/2012         ***
232 !********************************************************
233 interface
234   subroutine matrix_function(v_out,v_in,l)
235 
236   use defs_basis
237 
238   integer,     intent(in)   :: l
239   complex(dpc), intent(out) :: v_out(l)
240   complex(dpc), intent(in)  :: v_in(l)
241 
242   end subroutine matrix_function
243 end interface
244 
245 !------------------------------
246 ! input/output variables
247 !------------------------------
248 
249 integer, intent(in) :: mpi_communicator
250 integer, intent(in) :: kmax        ! number of Lanczos blocks
251 integer, intent(in) :: nseeds      ! size of each blocks
252 integer, intent(in) :: Hsize       ! size of the Hilbert space in which the matrix lives
253 
254 complex(dpc), intent(inout):: seeds(Hsize,nseeds) ! seed vectors for the algorithm
255 ! overwritten by X_{k+1} on output
256 
257 !logical,      intent(in) :: ortho           ! should the Lanczos vector be orthogonalized?
258 
259 complex(dpc), intent(out) :: alpha(nseeds,nseeds,kmax)  ! the alpha array from the Lanczos algorithm
260 complex(dpc), intent(out) :: beta(nseeds,nseeds,kmax)   ! the  beta array from the Lanczos algorithm
261 complex(dpc), intent(out) :: Lbasis(Hsize,nseeds*kmax)  ! array containing the Lanczos basis
262 
263 
264 complex(dpc), intent(in),optional :: X0(Hsize,nseeds)
265 complex(dpc), intent(in),optional :: beta0(nseeds,nseeds)
266 complex(dpc), intent(in),optional :: Qk(:,:)  ! array containing vectors to which
267 
268 ! the basis must be orthonormalized
269 
270 
271 
272 !------------------------------
273 ! local variables
274 !------------------------------
275 integer     :: k, seed1
276 integer     :: dum(2), lk
277 
278 complex(dpc), allocatable :: xk(:,:), xkm1(:,:), rk(:,:)
279 
280 integer     :: ntime, itime
281 real(dp)    :: total_time1, total_time2
282 real(dp)    :: time1, time2
283 integer     :: ierr
284 real(dp),allocatable :: list_time(:)
285 
286 ! *************************************************************************
287 
288 call cpu_time(total_time1)
289 
290 
291 ntime = 7
292 ABI_ALLOCATE(list_time,(ntime))
293 list_time(:) = zero
294 
295 if(present(Qk)) then
296   dum   = shape(Qk)
297   lk    = dum(2)
298 end if
299 
300 ABI_ALLOCATE( xk,  (Hsize,nseeds))
301 ABI_ALLOCATE( xkm1,(Hsize,nseeds))
302 ABI_ALLOCATE( rk  ,(Hsize,nseeds))
303 
304 
305 alpha = cmplx_0
306 beta  = cmplx_0
307 
308 !------------------------------------------------
309 ! Orthonormalize the seeds
310 !------------------------------------------------
311 ! initialize the xk array with the seeds
312 xk(:,:) = seeds(:,:)
313 
314 
315 ! orthonormalize the block using the QR algorithm
316 ! xk is overwritten by Q, the array of orthonormal vectors
317 
318 call extract_QR(mpi_communicator, Hsize,nseeds,xk)
319 
320 !------------------------------------------------
321 ! Loop on all blocks
322 !------------------------------------------------
323 
324 
325 do k = 1, kmax
326 
327 itime = 0
328 
329 ! tabulate basis, computed at previous step
330 Lbasis(:,nseeds*(k-1)+1:nseeds*k) = xk(:,:)
331 
332 
333 itime = itime+1
334 call cpu_time(time1)
335 
336 ! Initialize the residual array
337 do seed1 = 1, nseeds
338 ! If we are constructing the $\hat \epsilon(i\omega = 0)$ matrix (and the Lanczos basis at the same time),
339 ! note the index in which the Sternheimer solutions will be stored (for use in the projected Sternheimer section).
340 if(write_solution) index_solution = (k-1)*nseeds + seed1
341 
342 call matrix_function(rk(:,seed1),xk(:,seed1),Hsize)
343 end do
344 
345 call cpu_time(time2)
346 list_time(itime) = list_time(itime) + time2-time1
347 
348 itime = itime+1
349 call cpu_time(time1)
350 ! compute the alpha array, alpha = X^d.A.X
351 
352 call ZGEMM(              'C',   & ! take Hermitian conjugate of first array
353 'N',   & ! leave second array as is
354 nseeds,   & ! the number  of rows of the  matrix op( A )
355 nseeds,   & ! the number  of columns of the  matrix op( B )
356 Hsize,   & ! the number  of columns of the  matrix op( A ) == rows of matrix op( B )
357 cmplx_1,   & ! alpha constant
358 xk,   & ! matrix A
359 Hsize,   & ! LDA
360 rk,   & ! matrix B
361 Hsize,   & ! LDB
362 cmplx_0,   & ! beta constant
363 alpha(:,:,k),   & ! matrix C
364 nseeds)     ! LDC
365 call xmpi_sum(alpha(:,:,k),mpi_communicator,ierr) ! sum on all processors
366 
367 
368 
369 call cpu_time(time2)
370 list_time(itime) = list_time(itime) + time2-time1
371 
372 itime = itime+1
373 call cpu_time(time1)
374 ! update the residual array, rk = rk-X.alpha
375 call ZGEMM(              'N',   & ! leave first array as is
376 'N',   & ! leave second array as is
377 Hsize,   & ! the number  of rows of the  matrix op( A )
378 nseeds,   & ! the number  of columns of the  matrix op( B )
379 nseeds,   & ! the number  of columns of the  matrix op( A ) == rows of matrix op( B )
380 -cmplx_1,   & ! alpha constant
381 xk,   & ! matrix A
382 Hsize,   & ! LDA
383 alpha(:,:,k),   & ! matrix B
384 nseeds,   & ! LDB
385 cmplx_1,   & ! beta constant
386 rk,   & ! matrix C
387 Hsize)     ! LDC
388 
389 call cpu_time(time2)
390 list_time(itime) = list_time(itime) + time2-time1
391 
392 if (k .eq. 1 .and. present(X0) .and. present(beta0)) then
393   ! if k == 1, and X0,beta0 are present,
394   !  update the residual array, r1 = r1-X_{0}.beta^d_{0}
395   call ZGEMM(                'N',   & ! leave first array as is
396   'C',   & ! Hermitian conjugate the second array
397   Hsize,   & ! the number  of rows of the  matrix op( A )
398   nseeds,   & ! the number  of columns of the  matrix op( B )
399   nseeds,   & ! the number  of columns of the  matrix op( A ) == rows of matrix op( B )
400   -cmplx_1,   & ! alpha constant
401   X0,   & ! matrix A
402   Hsize,   & ! LDA
403   beta0(:,:),   & ! matrix B
404   nseeds,   & ! LDB
405   cmplx_1,   & ! beta constant
406   rk,   & ! matrix C
407   Hsize)     ! LDC
408 end if
409 
410 
411 itime = itime+1
412 call cpu_time(time1)
413 if (k .gt. 1) then
414 
415   ! if k > 1, update the residual array, rk = rk-X_{k-1}.beta^d_{k-1}
416   call ZGEMM(                'N',   & ! leave first array as is
417   'C',   & ! Hermitian conjugate the second array
418   Hsize,   & ! the number  of rows of the  matrix op( A )
419   nseeds,   & ! the number  of columns of the  matrix op( B )
420   nseeds,   & ! the number  of columns of the  matrix op( A ) == rows of matrix op( B )
421   -cmplx_1,   & ! alpha constant
422   xkm1,   & ! matrix A
423   Hsize,   & ! LDA
424   beta(:,:,k-1),   & ! matrix B
425   nseeds,   & ! LDB
426   cmplx_1,   & ! beta constant
427   rk,   & ! matrix C
428   Hsize)     ! LDC
429 
430 end if
431 call cpu_time(time2)
432 list_time(itime) = list_time(itime) + time2-time1
433 
434 ! store xk for next iteration
435 xkm1(:,:) = xk(:,:)
436 
437 
438 ! Orthonormalize THE RESIDUAL to all previously calculated directions 
439 
440 itime = itime+1
441 call cpu_time(time1)
442 !if ( ortho .and. (dtset%gwcalctyp/=1) ) then !This is a test to obtain the CPU time taken by the orthogonalizations.  
443 
444 if(present(Qk)) then
445   ! Orthonormalize to all previously calculated directions, if
446   ! this is a restarted Lanczos step
447   call orthogonalize(mpi_communicator, Hsize,lk,nseeds,Qk,rk)
448 end if
449 
450 call orthogonalize(mpi_communicator, Hsize,k*nseeds,nseeds,Lbasis(:,1:k*nseeds),rk)
451 
452 !end if
453 call cpu_time(time2)
454 list_time(itime) = list_time(itime) + time2-time1
455 
456 
457 itime = itime+1
458 call cpu_time(time1)
459 
460 ! perform QR decomposition to extract X_{k+1} and beta_{k}
461 call extract_QR(mpi_communicator, Hsize,nseeds,rk,beta(:,:,k))
462 
463 call cpu_time(time2)
464 list_time(itime) = list_time(itime) + time2-time1
465 ! copy the Q matrix (written on rk) in xk, which becomes X_{k+1}
466 xk(:,:) = rk(:,:)
467 
468 
469 end do !end loop on k
470 
471 ! overwrite the seeds with the last vector block.
472 seeds(:,:) = xk(:,:)
473 
474 ABI_DEALLOCATE( xk  )
475 ABI_DEALLOCATE( xkm1)
476 ABI_DEALLOCATE( rk  )
477 call cpu_time(total_time2)
478 
479 list_time(7) = total_time2-total_time1
480 
481 call write_block_lanczos_timing_log(list_time,ntime)
482 
483 ABI_DEALLOCATE(list_time)
484 
485 end subroutine block_lanczos_algorithm

m_gwls_GWlanczos/diagonalize_lanczos_banded [ Functions ]

[ Top ] [ m_gwls_GWlanczos ] [ Functions ]

NAME

  diagonalize_lanczos_banded

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputePoles,gwls_GenerateEpsilon,gwls_LanczosResolvents

CHILDREN

      zgemm,zhbev

SOURCE

507 subroutine diagonalize_lanczos_banded(kmax,nseeds,Hsize,alpha,beta,Lbasis,eigenvalues,debug)
508 !-----------------------------------------------------------------------------------
509 ! Given the result of the Lanczos algorithm, this subroutine diagonalize the banded
510 ! matrix as well as updates the basis.
511 !-----------------------------------------------------------------------------------
512 
513 !This section has been created automatically by the script Abilint (TD).
514 !Do not modify the following lines by hand.
515 #undef ABI_FUNC
516 #define ABI_FUNC 'diagonalize_lanczos_banded'
517 !End of the abilint section
518 
519 implicit none
520 
521 integer, intent(in)  :: kmax        ! number of Lanczos blocks
522 integer, intent(in)  :: nseeds      ! size of each blocks
523 integer, intent(in)  :: Hsize       ! size of the Hilbert space in which the matrix lives
524 logical, intent(in)  :: debug
525 
526 complex(dpc), intent(in) :: alpha(nseeds,nseeds,kmax)  ! the alpha array from the Lanczos algorithm
527 complex(dpc), intent(in) :: beta (nseeds,nseeds,kmax)  ! the  beta array from the Lanczos algorithm
528 
529 complex(dpc), intent(inout) :: Lbasis(Hsize,nseeds*kmax)  ! array containing the Lanczos basis
530 
531 
532 real(dp), intent(out) :: eigenvalues(nseeds*kmax)
533 
534 
535 ! local variables
536 
537 integer :: kd   ! number of superdiagonal above the diagonal in banded storage
538 integer :: ldab ! dimension of banded storage matrix
539 
540 complex(dpc), allocatable :: band_storage_matrix(:,:)
541 complex(dpc), allocatable :: saved_band_storage_matrix(:,:)
542 
543 complex(dpc), allocatable :: eigenvectors(:,:)
544 
545 complex(dpc), allocatable :: Lbasis_tmp(:,:)  
546 
547 integer :: i, j
548 integer :: k
549 integer :: s1, s2
550 integer :: info
551 
552 
553 complex(dpc), allocatable :: work(:)
554 real(dp),     allocatable :: rwork(:)
555 
556 integer        :: io_unit
557 character(128) :: filename
558 logical        :: file_exists
559 
560 integer        :: debug_unit
561 character(50)  :: debug_filename
562 
563 ! *************************************************************************
564 
565 
566 
567 ! number of superdiagonals
568 kd   = nseeds
569 ldab = kd + 1
570 
571 ABI_ALLOCATE(      band_storage_matrix, (ldab,nseeds*kmax))
572 ABI_ALLOCATE(saved_band_storage_matrix, (ldab,nseeds*kmax))
573 !---------------------------------------------------------
574 ! Store banded matrix in banded format
575 !---------------------------------------------------------
576 ! for UPLO = 'L', AB(1+i-j,j)    = A(i,j) for j<=i<=min(n,j+kd).
577 
578 band_storage_matrix(:,:) = cmplx_0
579 
580 !-----------------------------------------
581 ! Store the alpha and beta matrices
582 !-----------------------------------------
583 
584 ! loop on all blocks
585 do k=1,kmax
586 
587 ! alpha blocks
588 do s2 = 1, nseeds
589 j = (k-1)*nseeds+s2
590 
591 do s1 = s2, nseeds
592 i = j+s1-s2
593 band_storage_matrix(1+i-j,j) = alpha(s1,s2,k)
594 end do
595 end do
596 
597 ! exit when k = kmax, as this beta block does not contribute.
598 if  (k .eq. kmax) exit
599 
600 ! beta blocks
601 do s2 = 1, nseeds
602 j = (k-1)*nseeds+s2
603 
604 do s1 = 1, s2
605 i = j+s1-s2+nseeds
606 band_storage_matrix(1+i-j,j) = beta(s1,s2,k)
607 end do
608 
609 end do
610 end do
611 
612 saved_band_storage_matrix(:,:) = band_storage_matrix(:,:)
613 
614 !-----------------------------------------
615 ! Diagonalize the banded matrix
616 !-----------------------------------------
617 
618 ABI_ALLOCATE(eigenvectors, (nseeds*kmax,nseeds*kmax))
619 
620 ABI_ALLOCATE(work,(nseeds*kmax))
621 ABI_ALLOCATE(rwork,(3*nseeds*kmax-2))
622 
623 call ZHBEV(                     'V',      & ! compute eigenvalues and eigenvectors
624 'L',      & ! lower triangular part of matrix is stored in banded_matrix
625 nseeds*kmax,      & ! dimension of matrix
626 kd,      & ! number of superdiagonals in banded matrix
627 band_storage_matrix,      & ! matrix in banded storage
628 ldab,      & ! leading dimension of banded_matrix
629 eigenvalues,      & ! eigenvalues of matrix
630 eigenvectors,      & ! eigenvectors of matrix
631 nseeds*kmax,      & ! dimension of eigenvector matrix
632 work, rwork, info )  ! work arrays and info
633 
634 
635 if ( info /= 0) then        
636   debug_unit = get_unit()
637   write(debug_filename,'(A,I4.4,A)') 'LAPACK_DEBUG_PROC=',mpi_enreg%me,'.log'
638 
639   open(debug_unit,file=trim(debug_filename),status='unknown')
640 
641   write(debug_unit,'(A)')      '*********************************************************************************************'
642   write(debug_unit,'(A,I4,A)') '*      ERROR: info = ',info,' in ZHBEV (1), gwls_GWlanczos'
643   write(debug_unit,'(A)')      '*********************************************************************************************'
644 
645   close(debug_unit)
646 
647 end if
648 
649 
650 
651 !----------------------------------------------------------------------------------
652 ! update the Lanczos basis to reflect diagonalization of T matrix
653 !
654 !        Note that by definition
655 !
656 !                        Q^H . A . Q = T  ==>   A = Q . T . Q^H
657 !
658 !        where Q (Lbasis) contains the Lanczos basis.
659 !
660 !        Diagonalizing T, ie  T = U . LAMBDA . U^H leads to
661 !
662 !                         A  = [ Q.U] . LAMBDA . [Q.U]^H
663 !
664 !     The updated basis is thus Q.U == Lbasis . eigenvectors
665 !----------------------------------------------------------------------------------
666 
667 ! NEVER use matmul!!! It sends temporary arrays to the stack, which can be much smaller
668 ! than needed; this leads to mysterious segfaults!
669 
670 ! Lbasis = matmul(Lbasis,eigenvectors)
671 
672 ! use temporary array, which is PROPERLY ALLOCATED, to perform matrix multiplication
673 ABI_ALLOCATE(Lbasis_tmp, (Hsize,nseeds*kmax))  
674 
675 ! Compute C = A * B, where A = Lbasis, B = eigenvectors, and C = Lbasis_tmp
676 call ZGEMM(     'N',     & ! leave array A as is
677 'N',     & ! leave array B as is
678 Hsize,     & ! number of rows of A
679 nseeds*kmax,     & ! number of columns of B
680 nseeds*kmax,     & ! number of columns of A /rows of B
681 cmplx_1,     & ! constant alpha
682 Lbasis,     & ! matrix A
683 Hsize,     & ! LDA
684 eigenvectors,     & ! matrix B
685 nseeds*kmax,     & ! LDB
686 cmplx_0,     & ! constant beta
687 Lbasis_tmp,     & ! matrix C
688 Hsize)       ! LDC
689 
690 ! overwrite initial array
691 Lbasis(:,:) = Lbasis_tmp(:,:)
692 
693 
694 ABI_DEALLOCATE(Lbasis_tmp)
695 
696 if ( debug .and. mpi_enreg%me == 0)  then
697   !----------------------------------------------------------------------------------
698   ! For the purpose of debugging, print relevant results to file to check
699   ! data consistency. This may not be necessary once the code has been shown to
700   ! work properly.
701   !----------------------------------------------------------------------------------
702 
703   io_unit  = get_unit()
704   i = 0
705   file_exists = .true.
706   do while (file_exists)
707   i = i+1
708   write(filename,'(A,I0.4,A)') "diagonalize_banded_matrix_",i,".log"
709   inquire(file=filename,exist=file_exists)
710   end do
711 
712 
713   open(io_unit,file=filename,status=files_status_new)
714   write(io_unit,10) "#======================================================================================="
715   write(io_unit,10) "#                                                                                       "
716   write(io_unit,10) "#   This file contains information pertaining to the diagonalization of a banded        "
717   write(io_unit,10) "#   matrix, expressed in terms of the Lanczos alpha and beta block arrays.              "
718   write(io_unit,10) "#                                                                                       "
719   write(io_unit,10) "#======================================================================================="
720   write(io_unit,10) "#                                                                                       "
721   write(io_unit,12) "#   diagonalization info : ",info,"                                                     "
722   write(io_unit,10) "#                                                                                       "
723   write(io_unit,10) "#   l                 lambda_l                                                          "
724   write(io_unit,10) "#======================================================================================="
725 
726   do i = 1, nseeds*kmax
727   write(io_unit,13) i, eigenvalues(i)
728   end do
729 
730   write(io_unit,10) "                                                                                        "
731   write(io_unit,10) "#                                                                                       "
732   write(io_unit,10) "#    alpha and beta blocks                                                              "
733   write(io_unit,10) "#                                                                                       "
734   write(io_unit,10) "#======================================================================================="
735 
736   ! loop on all blocks
737   do k=1,kmax
738   write(io_unit,10) "#                        "
739   write(io_unit,12) "#   block k = ",k,"      "
740   write(io_unit,10) "#                        "
741   write(io_unit,15) "#                  alpha:   ||alpha^H-alpha|| = ",  &
742   sqrt(sum(abs(alpha(:,:,k)-transpose(conjg(alpha(:,:,k))))**2))
743   do s1 = 1, nseeds
744   write(io_unit,14) alpha(s1,:,k)
745   end do
746 
747   write(io_unit,10) "#                        "
748   write(io_unit,10) "#                  beta  "
749   do s1 = 1, nseeds
750   write(io_unit,14) beta(s1,:,k)
751   end do
752 
753 
754   end do
755 
756 
757 
758 
759   write(io_unit,10) "#                                                                                       "
760   write(io_unit,10) "#   band storage matrix:                                                                "
761   write(io_unit,10) "#======================================================================================="
762 
763   do i = 1, ldab
764   write(io_unit,14) saved_band_storage_matrix(i,:)
765   end do
766 
767   close(io_unit)
768 end if
769 
770 
771 ! clean up memory
772 ABI_DEALLOCATE(eigenvectors)
773 ABI_DEALLOCATE( work)
774 ABI_DEALLOCATE(rwork)
775 ABI_DEALLOCATE(band_storage_matrix)
776 ABI_DEALLOCATE(saved_band_storage_matrix)
777 
778 
779 10 format(A)
780 12 format(A,I5,A)
781 13 format(I5,5X,F24.12)
782 14 format(4X,1000(F12.8,SP,F12.8,1X,'i',2X))
783 15 format(A,ES24.8)
784 
785 end subroutine diagonalize_lanczos_banded

m_gwls_GWlanczos/get_seeds [ Functions ]

[ Top ] [ m_gwls_GWlanczos ] [ Functions ]

NAME

  get_seeds

FUNCTION

  .

INPUTS

OUTPUT

PARENTS

      gwls_ComputePoles,gwls_GenerateEpsilon

CHILDREN

      zgemm,zhbev

SOURCE

 88 subroutine get_seeds(first_seed, nseeds, seeds)
 89 !----------------------------------------------------------------------------------------------------
 90 !
 91 ! This subroutine compute the seeds using the eigenstates of the Hamiltonian
 92 !
 93 !----------------------------------------------------------------------------------------------------
 94 
 95 !This section has been created automatically by the script Abilint (TD).
 96 !Do not modify the following lines by hand.
 97 #undef ABI_FUNC
 98 #define ABI_FUNC 'get_seeds'
 99 !End of the abilint section
100 
101 implicit none
102 
103 integer,      intent(in)  :: first_seed, nseeds
104 complex(dpc), intent(out) :: seeds(npw_k,nseeds)
105 
106 real(dp)    , allocatable :: psik_out(:,:)
107 real(dp)    , allocatable :: psikb_e(:,:)
108 real(dp)    , allocatable :: psig_e(:,:)
109 real(dp)    , allocatable :: psikb_s(:,:)
110 real(dp)    , allocatable :: psig_s(:,:)
111 
112 
113 
114 ! local variables
115 integer  :: n
116 integer  :: i, j, nsblk
117 
118 ! *************************************************************************
119 
120 ! Generate the seeds for the Lanczos algorithm
121 ABI_ALLOCATE(psik_out,(2,npw_k))
122 ABI_ALLOCATE(psikb_e,(2,npw_kb))
123 ABI_ALLOCATE(psig_e,(2,npw_g))
124 ABI_ALLOCATE(psikb_s,(2,npw_kb))
125 ABI_ALLOCATE(psig_s,(2,npw_g))
126 
127 nsblk = ceiling(1.0*nseeds/blocksize) 
128 
129 do i=1,nsblk
130 do j=1,blocksize
131 psikb_e(:,(j-1)*npw_k+1:j*npw_k) = cg(:,(e-1)*npw_k+1:e*npw_k)
132 end do
133 
134 psig_e = zero
135 call wf_block_distribute(psikb_e,  psig_e,1) ! LA -> FFT 
136 
137 do j=1,blocksize
138 n = (i-1)*blocksize + j-1 + first_seed
139 if ((i-1)*blocksize + j <= nseeds) then
140   psikb_s(:,(j-1)*npw_k+1:j*npw_k) = cg(:,(n-1)*npw_k+1:n*npw_k)
141 else
142   psikb_s(:,(j-1)*npw_k+1:j*npw_k) = zero
143 end if
144 end do
145 
146 psig_s = zero
147 call wf_block_distribute(psikb_s,  psig_s,1) ! LA -> FFT 
148 
149 ! Fourier transform valence wavefunction, to real space
150 call g_to_r(psir1,psig_s)
151 
152 psir1(2,:,:,:) = -psir1(2,:,:,:)
153 
154 call gr_to_g(psig_s,psir1,psig_e)
155 
156 ! return to LA configuration, in order to apply Coulomb potential
157 call wf_block_distribute(psikb_s,  psig_s,2) ! FFT -> LA
158 
159 do j=1, blocksize
160 n = (i-1)*blocksize + j
161 if(n<=nseeds) then
162   psik_out = psikb_s(:,(j-1)*npw_k+1:j*npw_k)
163   call sqrt_vc_k(psik_out)
164   seeds(:,n) = cmplx_1*psik_out(1,:) + cmplx_i*psik_out(2,:) 
165 end if
166 end do
167 end do
168 
169 ABI_DEALLOCATE(psik_out)
170 ABI_DEALLOCATE(psikb_e)
171 ABI_DEALLOCATE(psig_e)
172 ABI_DEALLOCATE(psikb_s)
173 ABI_DEALLOCATE(psig_s)
174 
175 end subroutine get_seeds