LS: Add direct solver for sign method on submatrices

Since the submatrices are dense and we require the entire spectrum to
calculate the sign function, calculating the sign function using a
direct eigendecomposition is competitive to iterative methods.

This method can be used by setting SUBMATRIX_SIGN_METHOD to DIRECT. The
setting defaults to NEWTONSCHULZ. To perform the eigendecomposition, we
need a symmetric matrix. Therefore this method can only be used when
SIGN_SYMMETRIC is set as well.

Co-authored-by: Robert Schade <robert.schade@uni-paderborn.de>

src/dm_ls_scf.F

@@ -639,7 +639,7 @@ SUBROUTINE ls_scf_main(qs_env, ls_scf_env)
                      CALL density_matrix_sign(ls_scf_env%matrix_p(ispin), ls_scf_env%mu_spin(ispin), ls_scf_env%fixed_mu, &
                                               ls_scf_env%sign_method, ls_scf_env%sign_order, matrix_mixing_old(ispin), &
                                               ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, nelectron_spin_real, &
-                                              ls_scf_env%eps_filter, ls_scf_env%sign_symmetric)
+                                              ls_scf_env%eps_filter, ls_scf_env%sign_symmetric, ls_scf_env%submatrix_sign_method)
                   CASE (ls_scf_tc2)
                      CALL density_matrix_tc2(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, &
                                              nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), &
@@ -980,7 +980,8 @@ SUBROUTINE post_scf_mu_scan(ls_scf_env)
             CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, ls_scf_env%sign_method, &
                                               ls_scf_env%sign_order, ls_scf_env%matrix_ks(ispin), &
                                               ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, &
-                                              ls_scf_env%eps_filter, ls_scf_env%sign_symmetric)
+                                              ls_scf_env%eps_filter, ls_scf_env%sign_symmetric, &
+                                              ls_scf_env%submatrix_sign_method)
          ENDDO
 
          t2 = m_walltime()

src/dm_ls_scf_create.F

@@ -30,7 +30,8 @@ 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_sign_submatrix, ls_scf_tc2, ls_scf_trs4
+        ls_scf_sign_ns, ls_scf_sign_proot, ls_scf_sign_submatrix, ls_scf_submatrix_sign_direct, &
+        ls_scf_submatrix_sign_ns, ls_scf_tc2, ls_scf_trs4
    USE input_enumeration_types,         ONLY: enum_i2c,&
                                               enumeration_type
    USE input_keyword_types,             ONLY: keyword_get,&
@@ -210,6 +211,7 @@ SUBROUTINE ls_scf_init_read_write_input(input, ls_scf_env, unit_nr)
       CALL section_vals_val_get(ls_scf_section, "PERFORM_MU_SCAN", l_val=ls_scf_env%perform_mu_scan)
       CALL section_vals_val_get(ls_scf_section, "PURIFICATION_METHOD", i_val=ls_scf_env%purification_method)
       CALL section_vals_val_get(ls_scf_section, "SIGN_METHOD", i_val=ls_scf_env%sign_method)
+      CALL section_vals_val_get(ls_scf_section, "SUBMATRIX_SIGN_METHOD", i_val=ls_scf_env%submatrix_sign_method)
       CALL section_vals_val_get(ls_scf_section, "SIGN_ORDER", i_val=ls_scf_env%sign_order)
       CALL section_vals_val_get(ls_scf_section, "SIGN_SYMMETRIC", l_val=ls_scf_env%sign_symmetric)
       CALL section_vals_val_get(ls_scf_section, "DYNAMIC_THRESHOLD", l_val=ls_scf_env%dynamic_threshold)
@@ -310,6 +312,12 @@ SUBROUTINE ls_scf_init_read_write_input(input, ls_scf_env, unit_nr)
       IF (ls_scf_env%curvy_steps .AND. .NOT. ls_scf_env%use_s_sqrt) &
          CPABORT("CURVY_STEPS requires the use of the sqrt inversion.")
 
+      ! verify requirements for direct submatrix sign method
+      IF (ls_scf_env%sign_method .EQ. ls_scf_sign_submatrix &
+          .AND. ls_scf_env%submatrix_sign_method .EQ. ls_scf_submatrix_sign_direct &
+          .AND. .NOT. ls_scf_env%sign_symmetric) &
+         CPABORT("DIRECT submatrix sign method requires SIGN_SYMMETRIC being set.")
+
       ! an undocumented feature ... allows for just doing the initial guess, no expensive stuff
       IF (ls_scf_env%max_scf < 0) THEN
          ls_scf_env%needs_s_inv = .FALSE.
@@ -392,6 +400,14 @@ SUBROUTINE ls_scf_init_read_write_input(input, ls_scf_env, unit_nr)
                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")
+               SELECT CASE (ls_scf_env%submatrix_sign_method)
+               CASE (ls_scf_submatrix_sign_ns)
+                  WRITE (unit_nr, '(T2,A,T61,A20)') "Submatrix sign method:", ADJUSTR("newton schulz")
+               CASE (ls_scf_submatrix_sign_direct)
+                  WRITE (unit_nr, '(T2,A,T61,A20)') "Submatrix sign method:", ADJUSTR("direct")
+               CASE DEFAULT
+                  CPABORT("Unkown submatrix sign method.")
+               END SELECT
             CASE DEFAULT
                CPABORT("Unkown sign method.")
             END SELECT

src/dm_ls_scf_methods.F

@@ -370,12 +370,13 @@ END SUBROUTINE apply_matrix_preconditioner
 !> \param nelectron ...
 !> \param threshold ...
 !> \param sign_symmetric ...
+!> \param submatrix_sign_method ...
 !> \par History
 !>       2010.10 created [Joost VandeVondele]
 !> \author Joost VandeVondele
 ! **************************************************************************************************
    SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, &
-                                  matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric)
+                                  matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric, submatrix_sign_method)
 
       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_p
       REAL(KIND=dp), INTENT(INOUT)                       :: mu
@@ -385,6 +386,7 @@ SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order,
       INTEGER, INTENT(IN)                                :: nelectron
       REAL(KIND=dp), INTENT(IN)                          :: threshold
       LOGICAL                                            :: sign_symmetric
+      INTEGER                                            :: submatrix_sign_method
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_sign', &
          routineP = moduleN//':'//routineN
@@ -418,7 +420,7 @@ SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order,
 
          CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, &
                                            matrix_ks, matrix_s, matrix_s_inv, threshold, &
-                                           sign_symmetric)
+                                           sign_symmetric, submatrix_sign_method)
          IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,I2,1X,F13.9,1X,F15.9)') &
             "Density matrix:  iter, mu, trace error: ", iter, mu, trace - nelectron
 
@@ -455,12 +457,13 @@ END SUBROUTINE density_matrix_sign
 !> \param matrix_s_inv ...
 !> \param threshold ...
 !> \param sign_symmetric ...
+!> \param submatrix_sign_method ...
 !> \par History
 !>       2010.10 created [Joost VandeVondele]
 !> \author Joost VandeVondele
 ! **************************************************************************************************
    SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, matrix_ks, &
-                                           matrix_s, matrix_s_inv, threshold, sign_symmetric)
+                                           matrix_s, matrix_s_inv, threshold, sign_symmetric, submatrix_sign_method)
 
       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_p
       REAL(KIND=dp), INTENT(OUT)                         :: trace
@@ -469,6 +472,7 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o
       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_ks, matrix_s, matrix_s_inv
       REAL(KIND=dp), INTENT(IN)                          :: threshold
       LOGICAL                                            :: sign_symmetric
+      INTEGER                                            :: submatrix_sign_method
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_sign_fixed_mu', &
          routineP = moduleN//':'//routineN
@@ -516,7 +520,7 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o
          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)
+            CALL matrix_sign_submatrix(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, submatrix_sign_method)
          CASE DEFAULT
             CPABORT("Unkown sign method.")
          END SELECT
@@ -537,7 +541,7 @@ SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_o
          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)
+            CALL matrix_sign_submatrix(matrix_sign, matrix_sinv_ks, threshold, sign_order, submatrix_sign_method)
          CASE DEFAULT
             CPABORT("Unkown sign method.")
          END SELECT
@@ -1285,4 +1289,3 @@ FUNCTION estimate_steps(homo, lumo, threshold) RESULT(steps)
    END FUNCTION estimate_steps
 
 END MODULE dm_ls_scf_methods
-

src/dm_ls_scf_types.F

@@ -104,6 +104,7 @@ MODULE dm_ls_scf_types
       INTEGER :: sign_method
       INTEGER :: sign_order
       LOGICAL :: sign_symmetric
+      INTEGER :: submatrix_sign_method
       INTEGER :: s_sqrt_method
       INTEGER :: s_sqrt_order
 
@@ -220,4 +221,3 @@ SUBROUTINE ls_scf_release(ls_scf_env)
    END SUBROUTINE ls_scf_release
 
 END MODULE dm_ls_scf_types
-

src/energy_corrections.F

@@ -1543,7 +1543,7 @@ SUBROUTINE ec_ls_solver(qs_env, ls_method, matrix_ks, matrix_s, matrix_p, matrix
             CALL density_matrix_sign(matrix_p(ispin, 1)%matrix, mu_spin(ispin), fixed_mu, &
                                      sign_method, sign_order, matrix_ks(ispin, 1)%matrix, &
                                      matrix_s(1, 1)%matrix, matrix_s_inv, nelectron, eps_filter, &
-                                     sign_symmetric=.FALSE.)
+                                     sign_symmetric=.FALSE., submatrix_sign_method=0)
          CASE (ec_matrix_trs4)
             CALL density_matrix_trs4(matrix_p(ispin, 1)%matrix, matrix_ks(ispin, 1)%matrix, &
                                      matrix_s_sqrt_inv, nelectron, eps_filter, &

src/input_constants.F

@@ -887,6 +887,8 @@ MODULE input_constants
 
    INTEGER, PARAMETER, PUBLIC               :: ls_scf_sign_ns = 1, ls_scf_sign_proot = 2, ls_scf_sign_submatrix = 3
 
+   INTEGER, PARAMETER, PUBLIC               :: ls_scf_submatrix_sign_ns = 1, ls_scf_submatrix_sign_direct = 2
+
    INTEGER, PARAMETER, PUBLIC               :: ls_scf_line_search_3point = 3, &
                                                ls_scf_line_search_3point_2d = 6
    ! some ZMP paramenters

src/input_cp2k_ls.F

@@ -21,7 +21,8 @@ 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_sign_submatrix, ls_scf_tc2, ls_scf_trs4
+        ls_scf_sign_proot, ls_scf_sign_submatrix, ls_scf_submatrix_sign_direct, &
+        ls_scf_submatrix_sign_ns, ls_scf_tc2, ls_scf_trs4
    USE input_keyword_types,             ONLY: keyword_create,&
                                               keyword_release,&
                                               keyword_type
@@ -216,6 +217,17 @@ SUBROUTINE create_ls_scf_section(section)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
 
+      CALL keyword_create(keyword, __LOCATION__, name="SUBMATRIX_SIGN_METHOD", &
+                          description="Method used for the computation of the sign matrix of all submatrices.", &
+                          usage="SUBMATRIX_SIGN_METHOD NEWTONSCHULZ", &
+                          default_i_val=ls_scf_submatrix_sign_ns, &
+                          enum_c_vals=s2a("NEWTONSCHULZ", "DIRECT"), &
+                          enum_desc=s2a("Newton-Schulz iteration.", &
+                                        "Directly calculate all eigenvalues."), &
+                          enum_i_vals=(/ls_scf_submatrix_sign_ns, ls_scf_submatrix_sign_direct/))
+      CALL section_add_keyword(section, keyword)
+      CALL keyword_release(keyword)
+
       CALL keyword_create(keyword, __LOCATION__, name="SIGN_ORDER", &
                           description="Order of the method used for the computation of the sign matrix.", &
                           usage="SIGN_ORDER 2", &

src/iterate_matrix.F

@@ -23,6 +23,8 @@ MODULE iterate_matrix
         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 input_constants,                 ONLY: ls_scf_submatrix_sign_direct,&
+                                              ls_scf_submatrix_sign_ns
    USE kinds,                           ONLY: dp,&
                                               int_8
    USE machine,                         ONLY: m_flush,&
@@ -1080,7 +1082,7 @@ END SUBROUTINE matrix_sign_proot
 !> \param matrix ...
 !> \param threshold ...
 !> \param sign_order ...
-!> \author Michael Lass
+!> \author Michael Lass, Robert Schade
 ! **************************************************************************************************
    SUBROUTINE dense_matrix_sign_Newton_Schulz(matrix_sign, matrix, threshold, sign_order)
 
@@ -1153,22 +1155,88 @@ SUBROUTINE dense_matrix_sign_Newton_Schulz(matrix_sign, matrix, threshold, sign_
 
    END SUBROUTINE dense_matrix_sign_Newton_Schulz
 
+! **************************************************************************************************
+!> \brief Calculate the sign matrix by direct calculation of all eigenvalues and eigenvectors
+!> \param sm_sign ...
+!> \param sm ...
+!> \param N ...
+!> \author Michael Lass, Robert Schade
+! **************************************************************************************************
+   SUBROUTINE dense_matrix_sign_direct(sm_sign, sm, N)
+      INTEGER, INTENT(IN)                                :: N
+      REAL(KIND=dp), INTENT(IN)                          :: sm(N, N)
+      REAL(KIND=dp), INTENT(INOUT)                       :: sm_sign(N, N)
+
+      INTEGER                                            :: i, info, liwork, lwork
+      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: iwork
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: W, work
+      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: tmp, U
+
+      ALLOCATE (U(N, N), tmp(N, N))
+      ALLOCATE (W(N))
+
+      ! symmetrize sm
+      U(:, :) = 0.5*(sm + TRANSPOSE(sm))
+
+      ! probe optimal sizes for WORK and IWORK
+      LWORK = -1
+      LIWORK = -1
+      ALLOCATE (WORK(1))
+      ALLOCATE (IWORK(1))
+      CALL dsyevd('V', 'U', N, U, N, W, WORK, LWORK, IWORK, LIWORK, INFO)
+      LWORK = INT(WORK(1))
+      LIWORK = INT(IWORK(1))
+      DEALLOCATE (IWORK)
+      DEALLOCATE (WORK)
+
+      ! calculate eigenvalues and eigenvectors
+      ALLOCATE (WORK(LWORK))
+      ALLOCATE (IWORK(LIWORK))
+      CALL dsyevd('V', 'U', N, U, N, W, WORK, LWORK, IWORK, LIWORK, INFO)
+      DEALLOCATE (IWORK)
+      DEALLOCATE (WORK)
+
+      IF (INFO .NE. 0) PRINT *, "ERR: _syevd did not succeed", INFO
+
+      sm_sign = 0
+      DO i = 1, N
+         IF (W(i) > 0) THEN
+            sm_sign(i, i) = 1.0
+         ELSE IF (W(i) < 0) THEN
+            sm_sign(i, i) = -1.0
+         ELSE
+            sm_sign(i, i) = 0.0
+         ENDIF
+      ENDDO
+
+      ! Create matrix with eigenvalues in {-1,0,1} and eigenvectors of sm:
+      ! sm_sign = U * sm_sign * U.T
+      CALL dgemm('N', 'N', N, N, N, 1.0_dp, U, N, sm_sign, N, 0.0_dp, tmp, N)
+      CALL dgemm('N', 'T', N, N, N, 1.0_dp, tmp, N, U, N, 0.0_dp, sm_sign, N)
+
+      DEALLOCATE (W)
+      DEALLOCATE (tmp)
+      DEALLOCATE (U)
+   END SUBROUTINE dense_matrix_sign_direct
+
 ! **************************************************************************************************
 !> \brief Submatrix method
 !> \param matrix_sign ...
 !> \param matrix ...
 !> \param threshold ...
 !> \param sign_order ...
+!> \param submatrix_sign_method ...
 !> \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)
+   SUBROUTINE matrix_sign_submatrix(matrix_sign, matrix, threshold, sign_order, submatrix_sign_method)
 
       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_sign, matrix
       REAL(KIND=dp), INTENT(IN)                          :: threshold
       INTEGER, INTENT(IN), OPTIONAL                      :: sign_order
+      INTEGER, INTENT(IN)                                :: submatrix_sign_method
 
       CHARACTER(LEN=*), PARAMETER :: routineN = 'matrix_sign_submatrix', &
          routineP = moduleN//':'//routineN
@@ -1209,9 +1277,15 @@ SUBROUTINE matrix_sign_submatrix(matrix_sign, matrix, threshold, sign_order)
          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)
+         SELECT CASE (submatrix_sign_method)
+         CASE (ls_scf_submatrix_sign_ns)
+            CALL dense_matrix_sign_Newton_Schulz(sm_sign, sm, threshold, order)
+         CASE (ls_scf_submatrix_sign_direct)
+            CALL dense_matrix_sign_direct(sm_sign, sm, sm_size)
+         CASE DEFAULT
+            CPABORT("Unkown submatrix sign method.")
+         END SELECT
          CALL dissection%copy_resultcol(my_sms(i), sm_sign)
-
          DEALLOCATE (sm, sm_sign)
       ENDDO
       !$OMP END DO