Skip to content

Commit

Permalink
LS: Add sign calculation based on submatrix method
Browse files Browse the repository at this point in the history 
Use submatrix method to perform the sign iteration. This can be
activated by setting SIGN_METHOD to SUBMATRIX.

For now, Newton-Schulz iterations of order 2 and 3 are implemented for
the dense submatrices. This can easilly be extended by more iteration
schemes.

Co-authored-by: Robert Schade <robert.schade@uni-paderborn.de>
  • Loading branch information
@michaellass @robertschade @oschuett
2 people authored and oschuett committed Mar 12, 2020
1 parent c2469e9 commit fb312fe
Show file tree
Hide file tree
Showing 5 changed files with 176 additions and 20 deletions.
4 changes: 3 additions & 1 deletion src/dm_ls_scf_create.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ MODULE dm_ls_scf_create
ls_cluster_atomic, ls_cluster_molecular, ls_s_inversion_hotelling, ls_s_inversion_none, &
ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, &
ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_pexsi, ls_scf_sign, &
ls_scf_sign_ns, ls_scf_sign_proot, ls_scf_tc2, ls_scf_trs4
ls_scf_sign_ns, ls_scf_sign_proot, ls_scf_sign_submatrix, ls_scf_tc2, ls_scf_trs4
USE input_enumeration_types, ONLY: enum_i2c,&
enumeration_type
USE input_keyword_types, ONLY: keyword_get,&
Expand Down Expand Up @@ -390,6 +390,8 @@ SUBROUTINE ls_scf_init_read_write_input(input, ls_scf_env, unit_nr)
WRITE (unit_nr, '(T2,A,T61,A20)') "Sign method:", ADJUSTR("newton schulz")
CASE (ls_scf_sign_proot)
WRITE (unit_nr, '(T2,A,T61,A20)') "Sign method:", ADJUSTR("p-th root method")
CASE (ls_scf_sign_submatrix)
WRITE (unit_nr, '(T2,A,T61,A20)') "Sign method:", ADJUSTR("submatrix method")
CASE DEFAULT
CPABORT("Unkown sign method.")
END SELECT
Expand Down
17 changes: 9 additions & 8 deletions src/dm_ls_scf_methods.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,17 +25,14 @@ MODULE dm_ls_scf_methods
USE dm_ls_scf_types, ONLY: ls_cluster_atomic,&
ls_mstruct_type,&
ls_scf_env_type
USE input_constants, ONLY: ls_cluster_atomic,&
ls_s_preconditioner_atomic,&
ls_s_preconditioner_molecular,&
ls_s_preconditioner_none,&
ls_s_sqrt_ns,&
ls_s_sqrt_proot,&
ls_scf_sign_ns,&
ls_scf_sign_proot
USE input_constants, ONLY: &
ls_cluster_atomic, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, &
ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_sign_ns, &
ls_scf_sign_proot, ls_scf_sign_submatrix
USE iterate_matrix, ONLY: invert_Hotelling,&
matrix_sign_Newton_Schulz,&
matrix_sign_proot,&
matrix_sign_submatrix,&
matrix_sqrt_Newton_Schulz,&
matrix_sqrt_proot
USE kinds, ONLY: dp,&
Expand Down Expand Up @@ -518,6 +515,8 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o
CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
CASE (ls_scf_sign_proot)
CALL matrix_sign_proot(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
CASE (ls_scf_sign_submatrix)
CALL matrix_sign_submatrix(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
CASE DEFAULT
CPABORT("Unkown sign method.")
END SELECT
Expand All @@ -537,6 +536,8 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o
CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order)
CASE (ls_scf_sign_proot)
CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order)
CASE (ls_scf_sign_submatrix)
CALL matrix_sign_submatrix(matrix_sign, matrix_sinv_ks, threshold, sign_order)
CASE DEFAULT
CPABORT("Unkown sign method.")
END SELECT
Expand Down
2 changes: 1 addition & 1 deletion src/input_constants.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -885,7 +885,7 @@ MODULE input_constants
INTEGER, PARAMETER, PUBLIC :: ls_scf_sign = 17, ls_scf_trs4 = 18, &
ls_scf_tc2 = 19, ls_scf_pexsi = 20

INTEGER, PARAMETER, PUBLIC :: ls_scf_sign_ns = 1, ls_scf_sign_proot = 2
INTEGER, PARAMETER, PUBLIC :: ls_scf_sign_ns = 1, ls_scf_sign_proot = 2, ls_scf_sign_submatrix = 3

INTEGER, PARAMETER, PUBLIC :: ls_scf_line_search_3point = 3, &
ls_scf_line_search_3point_2d = 6
Expand Down
9 changes: 5 additions & 4 deletions src/input_cp2k_ls.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ MODULE input_cp2k_ls
ls_s_inversion_sign_sqrt, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, &
ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_line_search_3point, &
ls_scf_line_search_3point_2d, ls_scf_pexsi, ls_scf_sign, ls_scf_sign_ns, &
ls_scf_sign_proot, ls_scf_tc2, ls_scf_trs4
ls_scf_sign_proot, ls_scf_sign_submatrix, ls_scf_tc2, ls_scf_trs4
USE input_keyword_types, ONLY: keyword_create,&
keyword_release,&
keyword_type
Expand Down Expand Up @@ -208,10 +208,11 @@ SUBROUTINE create_ls_scf_section(section)
usage="SIGN_METHOD NEWTONSCHULZ", &
default_i_val=ls_scf_sign_ns, &
citations=(/VandeVondele2012, Niklasson2003/), &
enum_c_vals=s2a("NEWTONSCHULZ", "PROOT"), &
enum_c_vals=s2a("NEWTONSCHULZ", "PROOT", "SUBMATRIX"), &
enum_desc=s2a("Newton-Schulz iteration.", &
"p-th order root iteration"), &
enum_i_vals=(/ls_scf_sign_ns, ls_scf_sign_proot/))
"p-th order root iteration", &
"Submatrix method"), &
enum_i_vals=(/ls_scf_sign_ns, ls_scf_sign_proot, ls_scf_sign_submatrix/))
CALL section_add_keyword(section, keyword)
CALL keyword_release(keyword)

Expand Down
164 changes: 158 additions & 6 deletions src/iterate_matrix.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,18 +17,20 @@ MODULE iterate_matrix
cp_logger_get_default_unit_nr,&
cp_logger_type
USE dbcsr_api, ONLY: &
dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_filter, &
dbcsr_frobenius_norm, dbcsr_gershgorin_norm, dbcsr_get_diag, dbcsr_get_info, &
dbcsr_get_matrix_type, dbcsr_get_occupation, dbcsr_multiply, dbcsr_norm, &
dbcsr_norm_maxabsnorm, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, &
dbcsr_set_diag, dbcsr_trace, dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, &
dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_filter, dbcsr_frobenius_norm, &
dbcsr_gershgorin_norm, dbcsr_get_diag, dbcsr_get_info, dbcsr_get_matrix_type, &
dbcsr_get_occupation, dbcsr_multiply, dbcsr_norm, dbcsr_norm_maxabsnorm, dbcsr_p_type, &
dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_set_diag, dbcsr_trace, dbcsr_transposed, &
dbcsr_type, dbcsr_type_no_symmetry
USE kinds, ONLY: dp,&
int_8
USE machine, ONLY: m_flush,&
m_walltime
USE mathconstants, ONLY: ifac
USE mathlib, ONLY: abnormal_value
USE message_passing, ONLY: mp_sum
USE submatrix_dissection, ONLY: submatrix_dissection_type
#include "./base/base_uses.f90"

IMPLICIT NONE
Expand All @@ -42,7 +44,7 @@ MODULE iterate_matrix
END INTERFACE

PUBLIC :: invert_Hotelling, matrix_sign_Newton_Schulz, matrix_sqrt_Newton_Schulz, &
matrix_sqrt_proot, matrix_sign_proot, &
matrix_sqrt_proot, matrix_sign_proot, matrix_sign_submatrix, &
purify_mcweeny, invert_Taylor, determinant

CONTAINS
Expand Down Expand Up @@ -1072,6 +1074,156 @@ SUBROUTINE matrix_sign_proot(matrix_sign, matrix, threshold, sign_order)
END SUBROUTINE matrix_sign_proot
! **************************************************************************************************
!> \brief compute the sign of a dense matrix using Newton-Schulz iterations
!> \param matrix_sign ...
!> \param matrix ...
!> \param threshold ...
!> \param sign_order ...
!> \author Michael Lass
! **************************************************************************************************
SUBROUTINE dense_matrix_sign_Newton_Schulz(matrix_sign, matrix, threshold, sign_order)
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix_sign
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
REAL(KIND=dp), INTENT(IN) :: threshold
INTEGER, INTENT(IN), OPTIONAL :: sign_order
CHARACTER(LEN=*), PARAMETER :: routineN = 'dense_matrix_sign_Newton_Schulz', &
routineP = moduleN//':'//routineN
INTEGER :: handle, i, j, sz
REAL(KIND=dp) :: frob_matrix, frob_matrix_base, &
gersh_matrix, prefactor, scaling_factor
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp1, tmp2
REAL(KIND=dp), DIMENSION(1) :: work
REAL(KIND=dp), EXTERNAL :: dlange
CALL timeset(routineN, handle)
! scale the matrix to get into the convergence range
sz = SIZE(matrix, 1)
frob_matrix = dlange('F', sz, sz, matrix, sz, work) !dbcsr_frobenius_norm(matrix_sign)
gersh_matrix = dlange('1', sz, sz, matrix, sz, work) !dbcsr_gershgorin_norm(matrix_sign)
scaling_factor = 1/MIN(frob_matrix, gersh_matrix)
!print*,"norms",sz,frob_matrix,gersh_matrix,scaling_factor
matrix_sign = matrix*scaling_factor
ALLOCATE (tmp1(sz, sz))
ALLOCATE (tmp2(sz, sz))
DO i = 1, 100
CALL dgemm('N', 'N', sz, sz, sz, -1.0_dp, matrix_sign, sz, matrix_sign, sz, 0.0_dp, tmp1, sz)
! check convergence (frob norm of what should be the identity matrix minus identity matrix)
frob_matrix_base = dlange('F', sz, sz, tmp1, sz, work)
DO j = 1, sz
tmp1(j, j) = tmp1(j, j) + 1.0_dp
ENDDO
frob_matrix = dlange('F', sz, sz, tmp1, sz, work)
IF (sign_order .EQ. 2) THEN
prefactor = 0.5_dp
! update the above to 3*I-X*X
DO j = 1, sz
tmp1(j, j) = tmp1(j, j) + 2.0_dp
ENDDO
ELSE IF (sign_order .EQ. 3) THEN
tmp2(:, :) = tmp1
tmp1 = tmp1*4.0_dp/3.0_dp
DO j = 1, sz
tmp1(j, j) = tmp1(j, j) + 8.0_dp/3.0_dp
ENDDO
CALL dgemm('N', 'N', sz, sz, sz, 1.0_dp, tmp2, sz, tmp2, sz, 1.0_dp, tmp1, sz)
prefactor = 3.0_dp/8.0_dp
ELSE
CPABORT("requested order is not implemented.")
ENDIF
CALL dgemm('N', 'N', sz, sz, sz, prefactor, matrix_sign, sz, tmp1, sz, 0.0_dp, tmp2, sz)
matrix_sign = tmp2
! frob_matrix/frob_matrix_base < SQRT(threshold)
IF (frob_matrix*frob_matrix < (threshold*frob_matrix_base*frob_matrix_base)) EXIT
PRINT *, "Submatrix Sign iter", i, frob_matrix/frob_matrix_base
ENDDO
DEALLOCATE (tmp1)
DEALLOCATE (tmp2)
CALL timestop(handle)
END SUBROUTINE dense_matrix_sign_Newton_Schulz
! **************************************************************************************************
!> \brief Submatrix method
!> \param matrix_sign ...
!> \param matrix ...
!> \param threshold ...
!> \param sign_order ...
!> \par History
!> 2019.03 created [Robert Schade]
!> 2019.06 impl. submatrix method [Michael Lass]
!> \author Robert Schade, Michael Lass
! **************************************************************************************************
SUBROUTINE matrix_sign_submatrix(matrix_sign, matrix, threshold, sign_order)
TYPE(dbcsr_type), INTENT(INOUT) :: matrix_sign, matrix
REAL(KIND=dp), INTENT(IN) :: threshold
INTEGER, INTENT(IN), OPTIONAL :: sign_order
CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_sign_submatrix', &
routineP = moduleN//':'//routineN
INTEGER :: group, handle, i, myrank, nblkcols, &
order, sm_size, unit_nr
INTEGER, ALLOCATABLE, DIMENSION(:) :: my_sms
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sm, sm_sign
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_distribution_type) :: dist
TYPE(submatrix_dissection_type) :: dissection
CALL timeset(routineN, handle)
logger => cp_get_default_logger()
IF (logger%para_env%mepos == logger%para_env%source) THEN
unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
ELSE
unit_nr = -1
ENDIF
IF (PRESENT(sign_order)) THEN
order = sign_order
ELSE
order = 2
ENDIF
CALL dbcsr_get_info(matrix=matrix, nblkcols_total=nblkcols, distribution=dist, group=group)
CALL dbcsr_distribution_get(dist=dist, mynode=myrank)
CALL dissection%init(matrix)
CALL dissection%get_sm_ids_for_rank(myrank, my_sms)
!$OMP PARALLEL PRIVATE(sm, sm_sign, sm_size) SHARED(dissection, myrank, my_sms)
!$OMP DO SCHEDULE(GUIDED)
DO i = 1, SIZE(my_sms)
PRINT *, "Rank", myrank, "processing submatrix", my_sms(i)
CALL dissection%generate_submatrix(my_sms(i), sm)
sm_size = SIZE(sm, 1)
ALLOCATE (sm_sign(sm_size, sm_size))
CALL dense_matrix_sign_Newton_Schulz(sm_sign, sm, threshold, sign_order)
CALL dissection%copy_resultcol(my_sms(i), sm_sign)
DEALLOCATE (sm, sm_sign)
ENDDO
!$OMP END DO
!$OMP END PARALLEL
CALL dissection%communicate_results(matrix_sign)
CALL dissection%final
CALL timestop(handle)
END SUBROUTINE matrix_sign_submatrix
! **************************************************************************************************
!> \brief compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations
!> the order of the algorithm should be 2..5, 3 or 5 is recommended
Expand Down

0 comments on commit fb312fe

Please sign in to comment.