Line data    Source code

       1             : !-------------------------------------------------------------------------------
       2             : ! Copyright (c) 2023 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_elsi
       8             : CONTAINS
       9           0 :    SUBROUTINE elsi_diag(solver,hmat, smat, ne, eig, ev)
      10             :       
      11             :       USE m_juDFT
      12             :       USE m_types_mpimat
      13             :       USE m_types_mat
      14             : #ifdef CPP_ELSI
      15             :       USE elsi
      16             :       USE mpi
      17             : #endif
      18             :       IMPLICIT NONE
      19             : 
      20             :       CLASS(t_mat), INTENT(INOUT)    :: hmat, smat
      21             :       CLASS(t_mat), ALLOCATABLE, INTENT(OUT)::ev
      22             :       REAL, INTENT(out)              :: eig(:)
      23             :       INTEGER, INTENT(INOUT)         :: ne
      24             :       INTEGER,INTENT(IN)             :: solver
      25             : 
      26             : #ifdef CPP_ELSI
      27             :       !...  Local variables
      28             :       !
      29             :       INTEGER           :: blk,nev,myid,np,ierr,i
      30             :       INTEGER,PARAMETER :: BLACS_DENSE=0
      31             :       TYPE(elsi_handle) :: eh
      32             :       CLASS(t_mat),ALLOCATABLE :: evec
      33             :       real,allocatable         :: eig_tmp(:)
      34             : 
      35             : 
      36             :       call timestart("ELSI")
      37             : 
      38             :       call hmat%u2l()
      39             :       call smat%u2l()
      40             : 
      41             :       SELECT TYPE (hmat)
      42             :       TYPE IS (t_mpimat)
      43             :          SELECT TYPE (smat)
      44             :          TYPE IS (t_mpimat)
      45             :             if (hmat%blacsdata%blacs_desc(5)==hmat%blacsdata%blacs_desc(6)) then
      46             :                 blk=hmat%blacsdata%blacs_desc(5)
      47             :             else
      48             :                 call judft_error("BUG: in ELSI the row/column blocksize must be equal")
      49             :             endif    
      50             :             call elsi_init(eh,solver,1,BLACS_DENSE,hmat%global_size1,1.0*ne,ne)
      51             :             call elsi_set_mpi(eh,hmat%blacsdata%mpi_com)
      52             :             call elsi_set_blacs(eh,hmat%blacsdata%blacs_desc(2),blk)
      53             :             ALLOCATE(t_mpimat::evec)
      54             :             CALL evec%init(hmat)
      55             :             allocate(eig_tmp(hmat%global_size1))
      56             :          CLASS DEFAULT
      57             :             call judft_error("BUG: Inconsistent matrixes in ELSI call")
      58             :          END SELECT
      59             :       TYPE IS (t_mat)
      60             :           SELECT TYPE(smat)
      61             :           TYPE IS (t_mat)
      62             :               call elsi_init(eh,solver,0,BLACS_DENSE,hmat%matsize1,1.*ne,ne)
      63             :               allocate(t_mat::evec)
      64             :               call evec%init(hmat)
      65             :               allocate(eig_tmp(hmat%matsize1))
      66             : 
      67             :           CLASS DEFAULT
      68             :             call judft_error("BUG: Inconsistent matrixes in ELSI call")
      69             :          END SELECT  
      70             :       END SELECT
      71             :       call timestart("elsi_ev")
      72             : #ifdef _OPENACC
      73             :       call elsi_set_elpa_gpu(eh,1)
      74             : #endif
      75             :       !Now perform diagonalization
      76             :       call elsi_set_output(eh,3)
      77             :       call elsi_set_output_unit(eh,7)
      78             :       call elsi_set_illcond_check(eh, 0)
      79             :       call elsi_reinit(eh)
      80             :       IF (hmat%l_real) THEN
      81             :         call elsi_ev_real(eh,hmat%data_r,smat%data_r,eig_tmp,evec%data_r)
      82             :       ELSE  
      83             :         call elsi_ev_complex(eh,hmat%data_c,smat%data_c,eig_tmp,evec%data_c)
      84             :       ENDIF  
      85             :       call timestop("elsi_ev")
      86             :       call timestart("redistribute")
      87             :       !Copy data into correct data structures
      88             :       eig=eig_tmp(:size(eig))
      89             :       SELECT TYPE(evec)
      90             :       TYPE IS (t_mat)
      91             :          allocate(t_mat::ev)
      92             :          call ev%init(evec%l_real,evec%matsize1,ne)
      93             :          if (evec%l_real) THEN
      94             :             ev%data_r=evec%data_r(:,:ne)
      95             :          else   
      96             :             ev%data_c=evec%data_c(:,:ne)
      97             :          endif
      98             :       TYPE IS (t_mpimat)
      99             :          ALLOCATE(t_mpimat::ev)
     100             :          CALL ev%init(evec%l_real,evec%global_size1,evec%global_size1,evec%blacsdata%mpi_com,.FALSE.)
     101             :          CALL ev%copy(evec,1,1)
     102             :          !determine ev assigned to this rank
     103             :          nev=ne
     104             :          ne=0
     105             :          CALL MPI_COMM_RANK(evec%blacsdata%mpi_com,myid,ierr)
     106             :          CALL MPI_COMM_SIZE(evec%blacsdata%mpi_com,np,ierr)
     107             :          DO i=myid+1,nev,np
     108             :             ne=ne+1
     109             :          !   eig(ne)=eigenvalues(i)
     110             :          ENDDO
     111             : 
     112             : 
     113             :       END select   
     114             :       call evec%free()
     115             :       call timestop("redistribute")
     116             :       call elsi_finalize(eh)
     117             :       call timestop("ELSI")
     118             :      
     119             : #endif
     120           0 :    END SUBROUTINE elsi_diag
     121             : END MODULE m_elsi