Commit ca6b706f authored by Victor Yu's avatar Victor Yu

Refine dense DM extrapolation interface

parent 9a45b873
......@@ -302,13 +302,11 @@ void c_elsi_orthonormalize_ev_complex(elsi_handle handle_c,
double _Complex *evec);
void c_elsi_extrapolate_dm_real(elsi_handle handle_c,
double *ovlp0,
double *ovlp1,
double *ovlp,
double *dm);
void c_elsi_extrapolate_dm_complex(elsi_handle handle_c,
double _Complex *ovlp0,
double _Complex *ovlp1,
double _Complex *ovlp,
double _Complex *dm);
void c_elsi_construct_dm_real(elsi_handle handle_c,
......
......@@ -1545,19 +1545,17 @@ subroutine elsi_orthonormalize_ev_complex_c_wrapper(h_c,ovlp_c,evec_c)&
end subroutine
subroutine elsi_extrapolate_dm_real_c_wrapper(h_c,ovlp0_c,ovlp1_c,dm_c)&
subroutine elsi_extrapolate_dm_real_c_wrapper(h_c,ovlp_c,dm_c)&
bind(C,name="c_elsi_extrapolate_dm_real")
implicit none
type(c_ptr), value, intent(in) :: h_c
type(c_ptr), value, intent(in) :: ovlp0_c
type(c_ptr), value, intent(in) :: ovlp1_c
type(c_ptr), value, intent(in) :: ovlp_c
type(c_ptr), value, intent(in) :: dm_c
type(elsi_handle), pointer :: h_f
real(kind=c_double), pointer :: ovlp0_f(:,:)
real(kind=c_double), pointer :: ovlp1_f(:,:)
real(kind=c_double), pointer :: ovlp_f(:,:)
real(kind=c_double), pointer :: dm_f(:,:)
integer(kind=c_int) :: lrow
......@@ -1568,27 +1566,24 @@ subroutine elsi_extrapolate_dm_real_c_wrapper(h_c,ovlp0_c,ovlp1_c,dm_c)&
lrow = h_f%bh%n_lrow
lcol = h_f%bh%n_lcol
call c_f_pointer(ovlp0_c,ovlp0_f,shape=[lrow,lcol])
call c_f_pointer(ovlp1_c,ovlp1_f,shape=[lrow,lcol])
call c_f_pointer(ovlp_c,ovlp_f,shape=[lrow,lcol])
call c_f_pointer(dm_c,dm_f,shape=[lrow,lcol])
call elsi_extrapolate_dm_real(h_f,ovlp0_f,ovlp1_f,dm_f)
call elsi_extrapolate_dm_real(h_f,ovlp_f,dm_f)
end subroutine
subroutine elsi_extrapolate_dm_complex_c_wrapper(h_c,ovlp0_c,ovlp1_c,dm_c)&
subroutine elsi_extrapolate_dm_complex_c_wrapper(h_c,ovlp_c,dm_c)&
bind(C,name="c_elsi_extrapolate_dm_complex")
implicit none
type(c_ptr), value, intent(in) :: h_c
type(c_ptr), value, intent(in) :: ovlp0_c
type(c_ptr), value, intent(in) :: ovlp1_c
type(c_ptr), value, intent(in) :: ovlp_c
type(c_ptr), value, intent(in) :: dm_c
type(elsi_handle), pointer :: h_f
complex(kind=c_double), pointer :: ovlp0_f(:,:)
complex(kind=c_double), pointer :: ovlp1_f(:,:)
complex(kind=c_double), pointer :: ovlp_f(:,:)
complex(kind=c_double), pointer :: dm_f(:,:)
integer(kind=c_int) :: lrow
......@@ -1599,11 +1594,10 @@ subroutine elsi_extrapolate_dm_complex_c_wrapper(h_c,ovlp0_c,ovlp1_c,dm_c)&
lrow = h_f%bh%n_lrow
lcol = h_f%bh%n_lcol
call c_f_pointer(ovlp0_c,ovlp0_f,shape=[lrow,lcol])
call c_f_pointer(ovlp1_c,ovlp1_f,shape=[lrow,lcol])
call c_f_pointer(ovlp_c,ovlp_f,shape=[lrow,lcol])
call c_f_pointer(dm_c,dm_f,shape=[lrow,lcol])
call elsi_extrapolate_dm_complex(h_f,ovlp0_f,ovlp1_f,dm_f)
call elsi_extrapolate_dm_complex(h_f,ovlp_f,dm_f)
end subroutine
......
......@@ -9,8 +9,8 @@
!!
module ELSI_SETUP
use ELSI_CONSTANTS, only: PEXSI_SOLVER,SINGLE_PROC,MULTI_PROC,PEXSI_CSC,&
SIESTA_CSC,UNSET
use ELSI_CONSTANTS, only: ELPA_SOLVER,PEXSI_SOLVER,SINGLE_PROC,MULTI_PROC,&
PEXSI_CSC,SIESTA_CSC,UNSET
use ELSI_DATATYPE, only: elsi_handle,elsi_param_t,elsi_basic_t
use ELSI_ELPA, only: elsi_cleanup_elpa
use ELSI_IO, only: elsi_say,elsi_final_print,fjson_close_file,&
......@@ -400,12 +400,14 @@ subroutine elsi_reinit(eh)
call elsi_deallocate(eh%bh,eh%col_ptr_sp2,"col_ptr_sp2")
end if
if(allocated(eh%ovlp_real_copy)) then
call elsi_deallocate(eh%bh,eh%ovlp_real_copy,"ovlp_real_copy")
end if
if(.not. eh%ph%solver == ELPA_SOLVER) then
if(allocated(eh%ovlp_real_copy)) then
call elsi_deallocate(eh%bh,eh%ovlp_real_copy,"ovlp_real_copy")
end if
if(allocated(eh%ovlp_cmplx_copy)) then
call elsi_deallocate(eh%bh,eh%ovlp_cmplx_copy,"ovlp_cmplx_copy")
if(allocated(eh%ovlp_cmplx_copy)) then
call elsi_deallocate(eh%bh,eh%ovlp_cmplx_copy,"ovlp_cmplx_copy")
end if
end if
end if
......
This diff is collapsed.
......@@ -14,8 +14,7 @@ module ELSI_TOOLS
use ELSI_DATATYPE, only: elsi_handle
use ELSI_ELPA, only: elsi_update_dm_elpa
use ELSI_MALLOC, only: elsi_allocate
use ELSI_NTPOLY, only: elsi_init_ntpoly,elsi_update_dm_ntpoly,Matrix_ps,&
CopyMatrix,DestructMatrix
use ELSI_NTPOLY, only: elsi_init_ntpoly,elsi_update_dm_ntpoly,CopyMatrix
use ELSI_OCC, only: elsi_mu_and_occ,elsi_entropy
use ELSI_PRECISION, only: i4,r8
use ELSI_REDIST, only: elsi_blacs_to_ntpoly_hs,elsi_ntpoly_to_blacs_dm
......@@ -80,172 +79,89 @@ end subroutine
!>
!! This routine extrapolates density matrix for a new overlap.
!!
subroutine elsi_extrapolate_dm_real(eh,ovlp0,ovlp1,dm)
subroutine elsi_extrapolate_dm_real(eh,ovlp,dm)
implicit none
type(elsi_handle), intent(inout) :: eh !< Handle
real(kind=r8), intent(inout) :: ovlp0(eh%bh%n_lrow,eh%bh%n_lcol) !< Old overlap
real(kind=r8), intent(inout) :: ovlp1(eh%bh%n_lrow,eh%bh%n_lcol) !< New overlap
real(kind=r8), intent(inout) :: ovlp(eh%bh%n_lrow,eh%bh%n_lcol) !< New overlap
real(kind=r8), intent(inout) :: dm(eh%bh%n_lrow,eh%bh%n_lcol) !< Density matrix
integer(kind=i4) :: solver_save
real(kind=r8) :: dummy(1,1)
logical :: unit_ovlp_save
type(Matrix_ps) :: nt_ovlp0
character(len=*), parameter :: caller = "elsi_extrapolate_dm_real"
call elsi_check_init(eh%bh,eh%handle_init,caller)
solver_save = eh%ph%solver
if(eh%ph%solver == SIPS_SOLVER .and. eh%ph%sips_n_elpa > 0) then
solver_save = SIPS_SOLVER
eh%ph%solver = ELPA_SOLVER
! Overlap will be destroyed by Cholesky
if(.not. allocated(eh%ovlp_real_copy)) then
call elsi_allocate(eh%bh,eh%ovlp_real_copy,eh%bh%n_lrow,eh%bh%n_lcol,&
"ovlp_real_copy",caller)
eh%ovlp_real_copy = ovlp1
end if
end if
if(eh%ph%solver == OMM_SOLVER .and. eh%ph%omm_n_elpa > 0) then
solver_save = OMM_SOLVER
eh%ph%solver = ELPA_SOLVER
if(eh%ph%omm_flavor == 0) then
! Overlap will be destroyed by Cholesky
if(.not. allocated(eh%ovlp_real_copy)) then
call elsi_allocate(eh%bh,eh%ovlp_real_copy,eh%bh%n_lrow,&
eh%bh%n_lcol,"ovlp_real_copy",caller)
eh%ovlp_real_copy = ovlp1
end if
end if
end if
select case(eh%ph%solver)
case(ELPA_SOLVER)
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,ovlp0,ovlp1,dm)
case default
if(eh%ph%solver == NTPOLY_SOLVER) then
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp0,dummy,nt_ovlp0,&
eh%ph%nt_dm)
eh%ph%unit_ovlp = unit_ovlp_save
else
call elsi_init_ntpoly(eh%ph,eh%bh)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,ovlp0,eh%ph%nt_dm,nt_ovlp0)
end if
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,&
eh%ovlp_real_copy,ovlp,dm)
eh%ovlp_real_copy = ovlp
case default
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,ovlp,eh%ph%nt_ham,eh%ph%nt_dm)
eh%ph%first_blacs_to_ntpoly = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp1,dummy,eh%ph%nt_ovlp,&
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp,dm,eh%ph%nt_ovlp,&
eh%ph%nt_dm)
eh%ph%unit_ovlp = unit_ovlp_save
call elsi_update_dm_ntpoly(eh%ph,eh%bh,nt_ovlp0,eh%ph%nt_ovlp,&
eh%ph%nt_dm,eh%ph%nt_ham)
if(eh%ph%solver == NTPOLY_SOLVER) then
call DestructMatrix(nt_ovlp0)
end if
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_ham,dm)
call elsi_update_dm_ntpoly(eh%ph,eh%bh,eh%ph%nt_ovlp_copy,eh%ph%nt_ovlp,&
eh%ph%nt_ham,eh%ph%nt_dm)
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_dm,dm)
call CopyMatrix(eh%ph%nt_ovlp,eh%ph%nt_ovlp_copy)
end select
eh%ph%solver = solver_save
end subroutine
!>
!! This routine extrapolates density matrix for a new overlap.
!!
subroutine elsi_extrapolate_dm_complex(eh,ovlp0,ovlp1,dm)
subroutine elsi_extrapolate_dm_complex(eh,ovlp,dm)
implicit none
type(elsi_handle), intent(inout) :: eh !< Handle
complex(kind=r8), intent(inout) :: ovlp0(eh%bh%n_lrow,eh%bh%n_lcol) !< Old overlap
complex(kind=r8), intent(inout) :: ovlp1(eh%bh%n_lrow,eh%bh%n_lcol) !< New overlap
complex(kind=r8), intent(inout) :: ovlp(eh%bh%n_lrow,eh%bh%n_lcol) !< New overlap
complex(kind=r8), intent(inout) :: dm(eh%bh%n_lrow,eh%bh%n_lcol) !< Density matrix
integer(kind=i4) :: solver_save
complex(kind=r8) :: dummy(1,1)
logical :: unit_ovlp_save
type(Matrix_ps) :: nt_ovlp0
character(len=*), parameter :: caller = "elsi_extrapolate_dm_complex"
call elsi_check_init(eh%bh,eh%handle_init,caller)
solver_save = eh%ph%solver
if(eh%ph%solver == OMM_SOLVER .and. eh%ph%omm_n_elpa > 0) then
solver_save = OMM_SOLVER
eh%ph%solver = ELPA_SOLVER
if(eh%ph%omm_flavor == 0) then
! Overlap will be destroyed by Cholesky
if(.not. allocated(eh%ovlp_cmplx_copy)) then
call elsi_allocate(eh%bh,eh%ovlp_cmplx_copy,eh%bh%n_lrow,&
eh%bh%n_lcol,"ovlp_cmplx_copy",caller)
eh%ovlp_cmplx_copy = ovlp1
end if
end if
end if
select case(eh%ph%solver)
case(ELPA_SOLVER)
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,ovlp0,ovlp1,dm)
case default
if(eh%ph%solver == NTPOLY_SOLVER) then
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp0,dummy,nt_ovlp0,&
eh%ph%nt_dm)
eh%ph%unit_ovlp = unit_ovlp_save
else
call elsi_init_ntpoly(eh%ph,eh%bh)
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,ovlp0,eh%ph%nt_dm,nt_ovlp0)
end if
call elsi_update_dm_elpa(eh%ph,eh%bh,eh%row_map,eh%col_map,&
eh%ovlp_cmplx_copy,ovlp,dm)
eh%ovlp_cmplx_copy = ovlp
case default
unit_ovlp_save = eh%ph%unit_ovlp
eh%ph%unit_ovlp = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,dm,ovlp,eh%ph%nt_ham,eh%ph%nt_dm)
eh%ph%first_blacs_to_ntpoly = .true.
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp1,dummy,eh%ph%nt_ovlp,&
call elsi_blacs_to_ntpoly_hs(eh%ph,eh%bh,ovlp,dm,eh%ph%nt_ovlp,&
eh%ph%nt_dm)
eh%ph%unit_ovlp = unit_ovlp_save
call elsi_update_dm_ntpoly(eh%ph,eh%bh,nt_ovlp0,eh%ph%nt_ovlp,&
eh%ph%nt_dm,eh%ph%nt_ham)
if(eh%ph%solver == NTPOLY_SOLVER) then
call DestructMatrix(nt_ovlp0)
end if
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_ham,dm)
call elsi_update_dm_ntpoly(eh%ph,eh%bh,eh%ph%nt_ovlp_copy,eh%ph%nt_ovlp,&
eh%ph%nt_ham,eh%ph%nt_dm)
call elsi_ntpoly_to_blacs_dm(eh%bh,eh%ph%nt_dm,dm)
call CopyMatrix(eh%ph%nt_ovlp,eh%ph%nt_ovlp_copy)
end select
eh%ph%solver = solver_save
end subroutine
!>
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment