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