LCOV - code coverage report
Current view: top level - src - preconditioner_solvers.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 78.5 % 107 84
Test Date: 2025-07-25 12:55:17 Functions: 85.7 % 7 6

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief solves the preconditioner, contains to utility function for
      10              : !>        fm<->dbcsr transfers, should be moved soon
      11              : !> \par History
      12              : !>      - [UB] 2009-05-13 Adding stable approximate inverse (full and sparse)
      13              : !> \author Joost VandeVondele (09.2002)
      14              : ! **************************************************************************************************
      15              : MODULE preconditioner_solvers
      16              :    USE arnoldi_api,                     ONLY: arnoldi_env_type,&
      17              :                                               arnoldi_ev,&
      18              :                                               deallocate_arnoldi_env,&
      19              :                                               get_selected_ritz_val,&
      20              :                                               setup_arnoldi_env
      21              :    USE bibliography,                    ONLY: Schiffmann2015,&
      22              :                                               cite_reference
      23              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      24              :    USE cp_dbcsr_api,                    ONLY: &
      25              :         dbcsr_create, dbcsr_filter, dbcsr_get_info, dbcsr_get_occupation, dbcsr_init_p, &
      26              :         dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry
      27              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      28              :                                               copy_fm_to_dbcsr
      29              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_uplo_to_full
      30              :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose,&
      31              :                                               cp_fm_cholesky_invert
      32              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      33              :                                               cp_fm_struct_release,&
      34              :                                               cp_fm_struct_type
      35              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      36              :                                               cp_fm_release,&
      37              :                                               cp_fm_set_all,&
      38              :                                               cp_fm_type
      39              :    USE input_constants,                 ONLY: ot_precond_solver_default,&
      40              :                                               ot_precond_solver_direct,&
      41              :                                               ot_precond_solver_inv_chol,&
      42              :                                               ot_precond_solver_update
      43              :    USE iterate_matrix,                  ONLY: invert_Hotelling
      44              :    USE kinds,                           ONLY: dp
      45              :    USE message_passing,                 ONLY: mp_para_env_type
      46              :    USE preconditioner_types,            ONLY: preconditioner_type
      47              : #include "./base/base_uses.f90"
      48              : 
      49              :    IMPLICIT NONE
      50              : 
      51              :    PRIVATE
      52              : 
      53              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner_solvers'
      54              : 
      55              :    PUBLIC :: solve_preconditioner, transfer_fm_to_dbcsr, transfer_dbcsr_to_fm
      56              : 
      57              : CONTAINS
      58              : 
      59              : ! **************************************************************************************************
      60              : !> \brief ...
      61              : !> \param my_solver_type ...
      62              : !> \param preconditioner_env ...
      63              : !> \param matrix_s ...
      64              : !> \param matrix_h ...
      65              : ! **************************************************************************************************
      66         8500 :    SUBROUTINE solve_preconditioner(my_solver_type, preconditioner_env, matrix_s, &
      67              :                                    matrix_h)
      68              :       INTEGER                                            :: my_solver_type
      69              :       TYPE(preconditioner_type)                          :: preconditioner_env
      70              :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
      71              :       TYPE(dbcsr_type), POINTER                          :: matrix_h
      72              : 
      73              :       REAL(dp)                                           :: occ_matrix
      74              : 
      75              : ! here comes the solver
      76              : 
      77        13916 :       SELECT CASE (my_solver_type)
      78              :       CASE (ot_precond_solver_inv_chol)
      79              :          !
      80              :          ! compute the full inverse
      81         5416 :          preconditioner_env%solver = ot_precond_solver_inv_chol
      82         5416 :          CALL make_full_inverse_cholesky(preconditioner_env, matrix_s)
      83              :       CASE (ot_precond_solver_direct)
      84              :          !
      85              :          ! prepare for the direct solver
      86            0 :          preconditioner_env%solver = ot_precond_solver_direct
      87            0 :          CALL make_full_fact_cholesky(preconditioner_env, matrix_s)
      88              :       CASE (ot_precond_solver_update)
      89              :          !
      90              :          ! uses an update of the full inverse (needs to be computed the first time)
      91              :          ! make sure preconditioner_env is not destroyed in between
      92            6 :          occ_matrix = 1.0_dp
      93            6 :          IF (ASSOCIATED(preconditioner_env%sparse_matrix)) THEN
      94            6 :             IF (preconditioner_env%condition_num < 0.0_dp) &
      95            2 :                CALL estimate_cond_num(preconditioner_env%sparse_matrix, preconditioner_env%condition_num)
      96              :             CALL dbcsr_filter(preconditioner_env%sparse_matrix, &
      97            6 :                               1.0_dp/preconditioner_env%condition_num*0.01_dp)
      98            6 :             occ_matrix = dbcsr_get_occupation(preconditioner_env%sparse_matrix)
      99              :          END IF
     100              :          ! check whether we are in the first step and if it is a good idea to use cholesky (matrix sparsity)
     101            6 :          IF (preconditioner_env%solver .NE. ot_precond_solver_update .AND. occ_matrix > 0.5_dp) THEN
     102            2 :             preconditioner_env%solver = ot_precond_solver_update
     103            2 :             CALL make_full_inverse_cholesky(preconditioner_env, matrix_s)
     104              :          ELSE
     105            4 :             preconditioner_env%solver = ot_precond_solver_update
     106            4 :             CALL make_inverse_update(preconditioner_env, matrix_h)
     107              :          END IF
     108              :       CASE (ot_precond_solver_default)
     109         3078 :          preconditioner_env%solver = ot_precond_solver_default
     110              :       CASE DEFAULT
     111              :          !
     112         8500 :          CPABORT("Doesn't know this type of solver")
     113              :       END SELECT
     114              : 
     115         8500 :    END SUBROUTINE solve_preconditioner
     116              : 
     117              : ! **************************************************************************************************
     118              : !> \brief Compute the inverse using cholseky factorization
     119              : !> \param preconditioner_env ...
     120              : !> \param matrix_s ...
     121              : ! **************************************************************************************************
     122        16254 :    SUBROUTINE make_full_inverse_cholesky(preconditioner_env, matrix_s)
     123              : 
     124              :       TYPE(preconditioner_type)                          :: preconditioner_env
     125              :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
     126              : 
     127              :       CHARACTER(len=*), PARAMETER :: routineN = 'make_full_inverse_cholesky'
     128              : 
     129              :       INTEGER                                            :: handle, info
     130              :       TYPE(cp_fm_type)                                   :: fm_work
     131              :       TYPE(cp_fm_type), POINTER                          :: fm
     132              : 
     133         5418 :       CALL timeset(routineN, handle)
     134              : 
     135              :       ! Maybe we will get a sparse Cholesky at a given point then this can go,
     136              :       ! if stuff was stored in fm anyway this simple returns
     137              :       CALL transfer_dbcsr_to_fm(preconditioner_env%sparse_matrix, preconditioner_env%fm, &
     138         5418 :                                 preconditioner_env%para_env, preconditioner_env%ctxt)
     139         5418 :       fm => preconditioner_env%fm
     140              : 
     141         5418 :       CALL cp_fm_create(fm_work, fm%matrix_struct, name="fm_work")
     142              :       !
     143              :       ! compute the inverse of SPD matrix fm using the Cholesky factorization
     144         5418 :       CALL cp_fm_cholesky_decompose(fm, info_out=info)
     145              : 
     146              :       !
     147              :       ! if fm not SPD we go with the overlap matrix
     148         5418 :       IF (info /= 0) THEN
     149              :          !
     150              :          ! just the overlap matrix
     151            0 :          IF (PRESENT(matrix_s)) THEN
     152            0 :             CALL copy_dbcsr_to_fm(matrix_s, fm)
     153            0 :             CALL cp_fm_cholesky_decompose(fm)
     154              :          ELSE
     155            0 :             CALL cp_fm_set_all(fm, alpha=0._dp, beta=1._dp)
     156              :          END IF
     157              :       END IF
     158         5418 :       CALL cp_fm_cholesky_invert(fm)
     159              : 
     160         5418 :       CALL cp_fm_uplo_to_full(fm, fm_work)
     161         5418 :       CALL cp_fm_release(fm_work)
     162              : 
     163         5418 :       CALL timestop(handle)
     164              : 
     165         5418 :    END SUBROUTINE make_full_inverse_cholesky
     166              : 
     167              : ! **************************************************************************************************
     168              : !> \brief Only perform the factorization, can be used later to solve the linear
     169              : !>        system on the fly
     170              : !> \param preconditioner_env ...
     171              : !> \param matrix_s ...
     172              : ! **************************************************************************************************
     173            0 :    SUBROUTINE make_full_fact_cholesky(preconditioner_env, matrix_s)
     174              : 
     175              :       TYPE(preconditioner_type)                          :: preconditioner_env
     176              :       TYPE(dbcsr_type), OPTIONAL, POINTER                :: matrix_s
     177              : 
     178              :       CHARACTER(len=*), PARAMETER :: routineN = 'make_full_fact_cholesky'
     179              : 
     180              :       INTEGER                                            :: handle, info_out
     181              :       TYPE(cp_fm_type), POINTER                          :: fm
     182              : 
     183            0 :       CALL timeset(routineN, handle)
     184              : 
     185              :       ! Maybe we will get a sparse Cholesky at a given point then this can go,
     186              :       ! if stuff was stored in fm anyway this simple returns
     187              :       CALL transfer_dbcsr_to_fm(preconditioner_env%sparse_matrix, preconditioner_env%fm, &
     188            0 :                                 preconditioner_env%para_env, preconditioner_env%ctxt)
     189              : 
     190            0 :       fm => preconditioner_env%fm
     191              :       !
     192              :       ! compute the inverse of SPD matrix fm using the Cholesky factorization
     193            0 :       CALL cp_fm_cholesky_decompose(fm, info_out=info_out)
     194              :       !
     195              :       ! if fm not SPD we go with the overlap matrix
     196            0 :       IF (info_out .NE. 0) THEN
     197              :          !
     198              :          ! just the overlap matrix
     199            0 :          IF (PRESENT(matrix_s)) THEN
     200            0 :             CALL copy_dbcsr_to_fm(matrix_s, fm)
     201            0 :             CALL cp_fm_cholesky_decompose(fm)
     202              :          ELSE
     203            0 :             CALL cp_fm_set_all(fm, alpha=0._dp, beta=1._dp)
     204              :          END IF
     205              :       END IF
     206              : 
     207            0 :       CALL timestop(handle)
     208              : 
     209            0 :    END SUBROUTINE make_full_fact_cholesky
     210              : 
     211              : ! **************************************************************************************************
     212              : !> \brief computes an approximate inverse using Hotelling iterations
     213              : !> \param preconditioner_env ...
     214              : !> \param matrix_h as S is not always present this is a safe template for the transfer
     215              : ! **************************************************************************************************
     216            4 :    SUBROUTINE make_inverse_update(preconditioner_env, matrix_h)
     217              :       TYPE(preconditioner_type)                          :: preconditioner_env
     218              :       TYPE(dbcsr_type), POINTER                          :: matrix_h
     219              : 
     220              :       CHARACTER(len=*), PARAMETER :: routineN = 'make_inverse_update'
     221              : 
     222              :       INTEGER                                            :: handle
     223              :       LOGICAL                                            :: use_guess
     224              :       REAL(KIND=dp)                                      :: filter_eps
     225              : 
     226            4 :       CALL timeset(routineN, handle)
     227            4 :       use_guess = .TRUE.
     228              :       !
     229              :       ! uses an update of the full inverse (needs to be computed the first time)
     230              :       ! make sure preconditioner_env is not destroyed in between
     231              : 
     232            4 :       CALL cite_reference(Schiffmann2015)
     233              : 
     234              :       ! Maybe I gonna add a fm Hotelling, ... for now the same as above make sure we are dbcsr
     235            4 :       CALL transfer_fm_to_dbcsr(preconditioner_env%fm, preconditioner_env%sparse_matrix, matrix_h)
     236            4 :       IF (.NOT. ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN
     237            0 :          use_guess = .FALSE.
     238            0 :          CALL dbcsr_init_p(preconditioner_env%dbcsr_matrix)
     239              :          CALL dbcsr_create(preconditioner_env%dbcsr_matrix, "prec_dbcsr", &
     240            0 :                            template=matrix_h, matrix_type=dbcsr_type_no_symmetry)
     241              :       END IF
     242              : 
     243              :       ! Try to get a reasonbale guess for the filtering threshold
     244            4 :       filter_eps = 1.0_dp/preconditioner_env%condition_num*0.1_dp
     245              : 
     246              :       ! Aggressive filtering on the initial guess is needed to avoid fill ins and retain sparsity
     247            4 :       CALL dbcsr_filter(preconditioner_env%dbcsr_matrix, filter_eps*100.0_dp)
     248              :       ! We don't need a high accuracy for the inverse so 0.4 is reasonable for convergence
     249              :       CALL invert_Hotelling(preconditioner_env%dbcsr_matrix, preconditioner_env%sparse_matrix, filter_eps*10.0_dp, &
     250            4 :                             use_inv_as_guess=use_guess, norm_convergence=0.4_dp, filter_eps=filter_eps)
     251              : 
     252            4 :       CALL timestop(handle)
     253              : 
     254            4 :    END SUBROUTINE make_inverse_update
     255              : 
     256              : ! **************************************************************************************************
     257              : !> \brief computes an approximation to the condition number of a matrix using
     258              : !>        arnoldi iterations
     259              : !> \param matrix ...
     260              : !> \param cond_num ...
     261              : ! **************************************************************************************************
     262            2 :    SUBROUTINE estimate_cond_num(matrix, cond_num)
     263              :       TYPE(dbcsr_type), POINTER                          :: matrix
     264              :       REAL(KIND=dp)                                      :: cond_num
     265              : 
     266              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'estimate_cond_num'
     267              : 
     268              :       INTEGER                                            :: handle
     269              :       REAL(KIND=dp)                                      :: max_ev, min_ev
     270              :       TYPE(arnoldi_env_type)                             :: arnoldi_env
     271            2 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrices
     272              : 
     273            2 :       CALL timeset(routineN, handle)
     274              : 
     275              :       ! its better to do 2 calculations as the maximum should quickly converge and the minimum won't need iram
     276            4 :       ALLOCATE (matrices(1))
     277            2 :       matrices(1)%matrix => matrix
     278              :       ! compute the minimum ev
     279              :       CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=5.0E-4_dp, selection_crit=2, &
     280            2 :                              nval_request=1, nrestarts=15, generalized_ev=.FALSE., iram=.FALSE.)
     281            2 :       CALL arnoldi_ev(matrices, arnoldi_env)
     282            2 :       max_ev = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
     283            2 :       CALL deallocate_arnoldi_env(arnoldi_env)
     284              : 
     285              :       CALL setup_arnoldi_env(arnoldi_env, matrices, max_iter=20, threshold=5.0E-4_dp, selection_crit=3, &
     286            2 :                              nval_request=1, nrestarts=15, generalized_ev=.FALSE., iram=.FALSE.)
     287            2 :       CALL arnoldi_ev(matrices, arnoldi_env)
     288            2 :       min_ev = REAL(get_selected_ritz_val(arnoldi_env, 1), dp)
     289            2 :       CALL deallocate_arnoldi_env(arnoldi_env)
     290              : 
     291            2 :       cond_num = max_ev/min_ev
     292            2 :       DEALLOCATE (matrices)
     293              : 
     294            2 :       CALL timestop(handle)
     295            2 :    END SUBROUTINE estimate_cond_num
     296              : 
     297              : ! **************************************************************************************************
     298              : !> \brief transfers a full matrix to a dbcsr
     299              : !> \param fm_matrix a full matrix gets deallocated in the end
     300              : !> \param dbcsr_matrix a dbcsr matrix, gets create from a template
     301              : !> \param template_mat the template which is used for the structure
     302              : ! **************************************************************************************************
     303         7068 :    SUBROUTINE transfer_fm_to_dbcsr(fm_matrix, dbcsr_matrix, template_mat)
     304              : 
     305              :       TYPE(cp_fm_type), POINTER                          :: fm_matrix
     306              :       TYPE(dbcsr_type), POINTER                          :: dbcsr_matrix, template_mat
     307              : 
     308              :       CHARACTER(len=*), PARAMETER :: routineN = 'transfer_fm_to_dbcsr'
     309              : 
     310              :       INTEGER                                            :: handle
     311              : 
     312         7068 :       CALL timeset(routineN, handle)
     313         7068 :       IF (ASSOCIATED(fm_matrix)) THEN
     314         7056 :          IF (.NOT. ASSOCIATED(dbcsr_matrix)) THEN
     315         4000 :             CALL dbcsr_init_p(dbcsr_matrix)
     316              :             CALL dbcsr_create(dbcsr_matrix, template=template_mat, &
     317              :                               name="preconditioner_env%dbcsr_matrix", &
     318         4000 :                               matrix_type=dbcsr_type_no_symmetry)
     319              :          END IF
     320         7056 :          CALL copy_fm_to_dbcsr(fm_matrix, dbcsr_matrix)
     321         7056 :          CALL cp_fm_release(fm_matrix)
     322         7056 :          DEALLOCATE (fm_matrix)
     323              :          NULLIFY (fm_matrix)
     324              :       END IF
     325              : 
     326         7068 :       CALL timestop(handle)
     327              : 
     328         7068 :    END SUBROUTINE transfer_fm_to_dbcsr
     329              : 
     330              : ! **************************************************************************************************
     331              : !> \brief transfers a dbcsr to a full matrix
     332              : !> \param dbcsr_matrix a dbcsr matrix, gets deallocated at the end
     333              : !> \param fm_matrix a full matrix gets created if not yet done
     334              : !> \param para_env the para_env
     335              : !> \param context the blacs context
     336              : ! **************************************************************************************************
     337         6858 :    SUBROUTINE transfer_dbcsr_to_fm(dbcsr_matrix, fm_matrix, para_env, context)
     338              : 
     339              :       TYPE(dbcsr_type), POINTER                          :: dbcsr_matrix
     340              :       TYPE(cp_fm_type), POINTER                          :: fm_matrix
     341              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     342              :       TYPE(cp_blacs_env_type), POINTER                   :: context
     343              : 
     344              :       CHARACTER(len=*), PARAMETER :: routineN = 'transfer_dbcsr_to_fm'
     345              : 
     346              :       INTEGER                                            :: handle, n
     347              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     348              : 
     349         6858 :       CALL timeset(routineN, handle)
     350         6858 :       IF (ASSOCIATED(dbcsr_matrix)) THEN
     351         5418 :          NULLIFY (fm_struct_tmp)
     352              : 
     353         5418 :          IF (ASSOCIATED(fm_matrix)) THEN
     354            0 :             CALL cp_fm_release(fm_matrix)
     355            0 :             DEALLOCATE (fm_matrix)
     356              :          END IF
     357              : 
     358         5418 :          CALL dbcsr_get_info(dbcsr_matrix, nfullrows_total=n)
     359              :          CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n, ncol_global=n, &
     360         5418 :                                   context=context, para_env=para_env)
     361         5418 :          ALLOCATE (fm_matrix)
     362         5418 :          CALL cp_fm_create(fm_matrix, fm_struct_tmp)
     363         5418 :          CALL cp_fm_struct_release(fm_struct_tmp)
     364              : 
     365         5418 :          CALL copy_dbcsr_to_fm(dbcsr_matrix, fm_matrix)
     366         5418 :          CALL dbcsr_release(dbcsr_matrix)
     367         5418 :          DEALLOCATE (dbcsr_matrix)
     368              :       END IF
     369              : 
     370         6858 :       CALL timestop(handle)
     371              : 
     372         6858 :    END SUBROUTINE transfer_dbcsr_to_fm
     373              : 
     374              : END MODULE preconditioner_solvers
        

Generated by: LCOV version 2.0-1