|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00013 MODULE preconditioner 00014 USE array_types, ONLY: array_data,& 00015 array_i1d_obj 00016 USE cp_dbcsr_interface, ONLY: & 00017 cp_create_bl_distribution, cp_dbcsr_add, cp_dbcsr_add_on_diag, & 00018 cp_dbcsr_btriu, cp_dbcsr_checksum, cp_dbcsr_col_block_sizes, & 00019 cp_dbcsr_copy, cp_dbcsr_create, cp_dbcsr_distribution, & 00020 cp_dbcsr_distribution_release, cp_dbcsr_filter, cp_dbcsr_finalize, & 00021 cp_dbcsr_get_block, cp_dbcsr_get_block_p, cp_dbcsr_get_data_size, & 00022 cp_dbcsr_get_data_type, cp_dbcsr_get_info, cp_dbcsr_get_matrix_type, & 00023 cp_dbcsr_get_num_blocks, cp_dbcsr_get_occupation, & 00024 cp_dbcsr_get_stored_coordinates, cp_dbcsr_init, cp_dbcsr_init_p, & 00025 cp_dbcsr_iterator_blocks_left, cp_dbcsr_iterator_next_block, & 00026 cp_dbcsr_iterator_start, cp_dbcsr_iterator_stop, cp_dbcsr_multiply, & 00027 cp_dbcsr_nblkcols_local, cp_dbcsr_nblkrows_local, & 00028 cp_dbcsr_nblkrows_total, cp_dbcsr_norm, cp_dbcsr_print, & 00029 cp_dbcsr_put_block, cp_dbcsr_release, cp_dbcsr_release_p, & 00030 cp_dbcsr_row_block_sizes, cp_dbcsr_trace, cp_dbcsr_verify_matrix, & 00031 cp_dbcsr_work_create 00032 USE cp_dbcsr_operations, ONLY: & 00033 copy_dbcsr_to_fm, copy_fm_to_dbcsr, cp_dbcsr_copy_vec, & 00034 cp_dbcsr_deallocate_matrix, cp_dbcsr_from_fm, cp_dbcsr_multiply_vec, & 00035 cp_dbcsr_pack_vec, cp_dbcsr_sm_fm_multiply, cp_dbcsr_unpack_vec, & 00036 packed_vec_bcast, packed_vec_bif_tech2, packed_vec_build_u, & 00037 packed_vec_ini, packed_vec_scale 00038 USE cp_dbcsr_types, ONLY: cp_dbcsr_iterator,& 00039 cp_dbcsr_p_type,& 00040 cp_dbcsr_type 00041 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,& 00042 cp_fm_gemm,& 00043 cp_fm_scale_and_add,& 00044 cp_fm_trace,& 00045 cp_fm_upper_to_full 00046 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,& 00047 cp_fm_cholesky_invert,& 00048 cp_fm_cholesky_reduce,& 00049 cp_fm_cholesky_restore 00050 USE cp_fm_diag, ONLY: choose_eigv_solver 00051 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 00052 cp_fm_struct_release,& 00053 cp_fm_struct_type 00054 USE cp_fm_types, ONLY: cp_fm_create,& 00055 cp_fm_get_diag,& 00056 cp_fm_get_info,& 00057 cp_fm_release,& 00058 cp_fm_set_all,& 00059 cp_fm_to_fm,& 00060 cp_fm_type 00061 USE dbcsr_block_operations, ONLY: block_add_on_diag,& 00062 block_chol_inv 00063 USE dbcsr_error_handling, ONLY: dbcsr_error_type 00064 USE dbcsr_methods, ONLY: & 00065 dbcsr_distribution_mp, dbcsr_distribution_new, & 00066 dbcsr_distribution_row_dist, dbcsr_mp_group, dbcsr_mp_mynode, & 00067 dbcsr_mp_mypcol, dbcsr_mp_myprow, dbcsr_mp_npcols, dbcsr_mp_nprows, & 00068 dbcsr_mp_numnodes 00069 USE dbcsr_operations, ONLY: dbcsr_scale_mat 00070 USE dbcsr_types, ONLY: dbcsr_distribution_obj,& 00071 dbcsr_mp_obj,& 00072 dbcsr_norm_frobenius,& 00073 dbcsr_type_no_symmetry,& 00074 dbcsr_type_real_4,& 00075 dbcsr_type_real_default 00076 USE input_constants, ONLY: & 00077 ot_precond_full_all, ot_precond_full_kinetic, ot_precond_full_single, & 00078 ot_precond_full_single_inverse, ot_precond_none, ot_precond_s_inverse, & 00079 ot_precond_solver_default, ot_precond_solver_direct, & 00080 ot_precond_solver_inv_chol, ot_precond_solver_sainv, & 00081 ot_precond_sparse_diag, ot_precond_sparse_kinetic, & 00082 ot_precond_sparse_kinetic_sainv 00083 USE kinds, ONLY: dp,& 00084 dp_size,& 00085 int_size,& 00086 sp 00087 USE message_passing, ONLY: mp_bcast,& 00088 mp_max 00089 USE preconditioner_types, ONLY: destroy_preconditioner,& 00090 init_preconditioner,& 00091 preconditioner_p_type,& 00092 preconditioner_type 00093 USE qs_environment_types, ONLY: qs_environment_type 00094 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues 00095 USE qs_mo_types, ONLY: get_mo_set,& 00096 mo_set_p_type,& 00097 mo_set_type 00098 USE termination, ONLY: stop_memory,& 00099 stop_program 00100 USE timings, ONLY: timeset,& 00101 timestop 00102 #include "cp_common_uses.h" 00103 00104 IMPLICIT NONE 00105 00106 PRIVATE 00107 00108 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'preconditioner' 00109 00110 PUBLIC :: make_preconditioner, restart_preconditioner 00111 PUBLIC :: apply_preconditioner, prepare_preconditioner, apply_full_all 00112 00113 PRIVATE :: make_sparse_diag, make_full_single, make_local_block 00114 00115 INTERFACE apply_preconditioner 00116 MODULE PROCEDURE apply_preconditioner_dbcsr 00117 MODULE PROCEDURE apply_preconditioner_fm 00118 END INTERFACE 00119 00120 INTERFACE apply_solve_lin_system 00121 MODULE PROCEDURE apply_solve_lin_system_dbcsr 00122 MODULE PROCEDURE apply_solve_lin_system_fm 00123 END INTERFACE 00124 00125 ! ***************************************************************************** 00126 00127 CONTAINS 00128 00129 ! ***************************************************************************** 00130 00131 ! creates a preconditioner for the system (H-energy_homo S) 00132 ! this preconditioner is (must be) symmetric positive definite. 00133 ! currently uses a atom-block-diagonal form 00134 ! each block will be .... 00135 ! might overwrite matrix_h, matrix_t 00136 00137 ! ***************************************************************************** 00138 SUBROUTINE make_preconditioner(preconditioner_env, precon_type, solver_type, matrix_h, matrix_s, & 00139 matrix_t, mo_set, energy_gap, mixed_precision, convert_precond_to_dbcsr, error) 00140 00141 TYPE(preconditioner_type) :: preconditioner_env 00142 INTEGER, INTENT(IN) :: precon_type, solver_type 00143 TYPE(cp_dbcsr_type), POINTER :: matrix_h 00144 TYPE(cp_dbcsr_type), OPTIONAL, POINTER :: matrix_s, matrix_t 00145 TYPE(mo_set_type), POINTER :: mo_set 00146 REAL(KIND=dp) :: energy_gap 00147 LOGICAL, INTENT(IN), OPTIONAL :: mixed_precision, 00148 convert_precond_to_dbcsr 00149 TYPE(cp_error_type), INTENT(inout) :: error 00150 00151 CHARACTER(len=*), PARAMETER :: routineN = 'make_preconditioner', 00152 routineP = moduleN//':'//routineN 00153 00154 INTEGER :: handle, istat, k, 00155 my_solver_type 00156 LOGICAL :: failure, my_convert_precond_to_dbcsr, my_mixed_precision, 00157 needs_full_spectrum, needs_homo, use_mo_coeff_b 00158 REAL(KIND=dp) :: energy_homo 00159 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues_ot 00160 TYPE(cp_dbcsr_type), POINTER :: mo_coeff_b 00161 TYPE(cp_fm_type), POINTER :: mo_coeff 00162 00163 failure=.FALSE. 00164 00165 CALL timeset(routineN,handle) 00166 00167 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, mo_coeff_b=mo_coeff_b) 00168 use_mo_coeff_b = mo_set%use_mo_coeff_b 00169 00170 CALL cp_fm_get_info(mo_coeff,ncol_global=k,error=error) 00171 00172 my_convert_precond_to_dbcsr = .FALSE. 00173 IF(PRESENT(convert_precond_to_dbcsr)) my_convert_precond_to_dbcsr = convert_precond_to_dbcsr 00174 00175 my_mixed_precision = .FALSE. 00176 IF(PRESENT(mixed_precision)) my_mixed_precision = mixed_precision 00177 IF(my_mixed_precision) THEN 00178 SELECT CASE(precon_type) 00179 CASE(ot_precond_full_kinetic,ot_precond_full_single_inverse) 00180 !supported 00181 CASE DEFAULT 00182 CALL stop_program(routineN,moduleN,__LINE__,& 00183 "This precond with mixed precision is not supported yet") 00184 END SELECT 00185 ENDIF 00186 00187 needs_full_spectrum=.FALSE. 00188 needs_homo=.FALSE. 00189 00190 SELECT CASE(precon_type) 00191 CASE (ot_precond_full_single_inverse,ot_precond_full_all) 00192 needs_full_spectrum=.TRUE. 00193 CASE (ot_precond_sparse_diag,ot_precond_full_single) 00194 needs_homo=.TRUE. 00195 ! XXXX to be removed if homo estimate only is implemented 00196 needs_full_spectrum=.TRUE. 00197 CASE (ot_precond_full_kinetic,ot_precond_s_inverse,ot_precond_sparse_kinetic) 00198 ! these should be happy without an estimate for the homo energy 00199 ! preconditioning can not depend on an absolute eigenvalue, only on eigenvalue differences 00200 CASE DEFAULT 00201 CALL stop_program(routineN,moduleN,__LINE__,& 00202 "The preconditioner is unknown ...") 00203 END SELECT 00204 00205 energy_homo=0.0_dp 00206 IF (needs_full_spectrum) THEN 00207 ALLOCATE(eigenvalues_ot(k),STAT=istat) 00208 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00209 ! XXXXXXXXXXXXXXXX do not touch the initial MOs, could be harmful for either 00210 ! the case of non-equivalent MOs but also for the derivate 00211 ! we could already have all eigenvalues e.g. full_all and we could skip this 00212 ! to be optimised later. 00213 ! one flaw is that not all SCF methods (i.e. that go over mo_derivs directly) 00214 ! have a 'valid' matrix_h... (we even don't know what evals are in that case) 00215 IF(use_mo_coeff_b) THEN 00216 CALL calculate_subspace_eigenvalues(mo_coeff_b,matrix_h,& 00217 eigenvalues_ot, do_rotation = .FALSE.,& 00218 para_env=mo_coeff%matrix_struct%para_env,& 00219 blacs_env=mo_coeff%matrix_struct%context,error=error) 00220 ELSE 00221 CALL calculate_subspace_eigenvalues(mo_coeff,matrix_h,& 00222 eigenvalues_ot, do_rotation = .FALSE.,error=error) 00223 ENDIF 00224 IF (k>0) energy_homo=eigenvalues_ot(k) 00225 ELSE 00226 IF (needs_homo) THEN 00227 CALL stop_program(routineN,moduleN,__LINE__,& 00228 "Not yet implemented") 00229 ENDIF 00230 ENDIF 00231 00232 my_solver_type = solver_type 00233 SELECT CASE (precon_type) 00234 CASE (ot_precond_sparse_diag) 00235 preconditioner_env%in_use=ot_precond_sparse_diag 00236 IF (my_solver_type.NE.ot_precond_solver_default) THEN 00237 CALL stop_program(routineN,moduleN,__LINE__,& 00238 "Only PRECOND_SOLVER DEFAULT for the moment") 00239 END IF 00240 IF ( PRESENT(matrix_s) ) THEN 00241 CALL make_sparse_diag(preconditioner_env, matrix_h, matrix_s, & 00242 energy_homo, energy_gap ,error=error) 00243 ELSE 00244 CALL make_sparse_diag_ortho(preconditioner_env, matrix_h, & 00245 energy_homo, energy_gap ,error=error) 00246 END IF 00247 CASE (ot_precond_full_single) 00248 preconditioner_env%in_use=ot_precond_full_single 00249 IF(my_solver_type.NE.ot_precond_solver_default) THEN 00250 CALL stop_program(routineN,moduleN,__LINE__,& 00251 "Only PRECOND_SOLVER DEFAULT for the moment") 00252 ENDIF 00253 IF ( PRESENT(matrix_s) ) THEN 00254 CALL make_full_single(preconditioner_env, preconditioner_env%fm,& 00255 matrix_h, matrix_s, energy_homo, energy_gap ,error=error) 00256 ELSE 00257 CALL make_full_single_ortho(preconditioner_env, preconditioner_env%fm,& 00258 matrix_h, energy_homo, energy_gap ,error=error) 00259 END IF 00260 CASE (ot_precond_s_inverse) 00261 preconditioner_env%in_use=ot_precond_s_inverse 00262 !IF(my_solver_type.NE.ot_precond_solver_default) THEN 00263 ! CALL stop_program(routineN,moduleN,__LINE__,"Only PRECOND_SOLVER DEFAULT for the moment") 00264 !ENDIF 00265 IF(my_solver_type.EQ.ot_precond_solver_default) my_solver_type=ot_precond_solver_inv_chol 00266 IF ( PRESENT(matrix_s) ) THEN 00267 CALL make_full_s_inverse(preconditioner_env, matrix_h, matrix_s, error=error) 00268 ELSE 00269 CALL stop_program(routineN,moduleN,__LINE__,& 00270 "Type for S=1 not implemented") 00271 END IF 00272 CASE (ot_precond_full_kinetic) 00273 preconditioner_env%in_use=ot_precond_full_kinetic 00274 IF(my_solver_type.EQ.ot_precond_solver_default) my_solver_type=ot_precond_solver_inv_chol 00275 IF ( PRESENT(matrix_s) .AND. PRESENT(matrix_t) ) THEN 00276 CALL make_full_kinetic(preconditioner_env, preconditioner_env%fm,& 00277 & matrix_t, matrix_s, energy_gap, & 00278 & my_mixed_precision, error=error) 00279 ELSE 00280 CALL stop_program(routineN,moduleN,__LINE__,& 00281 "Type for S=1 not implemented") 00282 ENDIF 00283 CASE (ot_precond_full_single_inverse) 00284 preconditioner_env%in_use=ot_precond_full_single_inverse 00285 IF(my_solver_type.EQ.ot_precond_solver_default) my_solver_type=ot_precond_solver_inv_chol 00286 IF(use_mo_coeff_b) THEN 00287 CALL copy_dbcsr_to_fm(mo_coeff_b,mo_coeff,error=error) 00288 ENDIF 00289 IF ( PRESENT(matrix_s) ) THEN 00290 CALL make_full_single_inverse(preconditioner_env,mo_coeff,matrix_h, matrix_s, & 00291 eigenvalues_ot, energy_gap, my_mixed_precision, & 00292 error=error) 00293 ELSE 00294 CALL make_full_single_inverse_ortho(preconditioner_env,mo_coeff,matrix_h, & 00295 eigenvalues_ot, energy_gap,error=error) 00296 END IF 00297 CASE (ot_precond_full_all) 00298 preconditioner_env%in_use=ot_precond_full_all 00299 IF(my_solver_type.NE.ot_precond_solver_default) THEN 00300 CALL stop_program(routineN,moduleN,__LINE__,& 00301 "Only PRECOND_SOLVER DEFAULT for the moment") 00302 ENDIF 00303 IF(use_mo_coeff_b) THEN 00304 CALL copy_dbcsr_to_fm(mo_coeff_b,mo_coeff,error=error) 00305 ENDIF 00306 IF ( PRESENT(matrix_s) ) THEN 00307 CALL make_full_all(preconditioner_env,mo_coeff,matrix_h, matrix_s, & 00308 eigenvalues_ot, energy_gap,error=error) 00309 ELSE 00310 CALL make_full_all_ortho(preconditioner_env,mo_coeff,matrix_h, & 00311 eigenvalues_ot, energy_gap,error=error) 00312 END IF 00313 CASE(ot_precond_sparse_kinetic) 00314 preconditioner_env%in_use=ot_precond_sparse_kinetic 00315 IF(my_solver_type.NE.ot_precond_solver_default) THEN 00316 CALL stop_program(routineN,moduleN,__LINE__,& 00317 "Only PRECOND_SOLVER DEFAULT for the moment") 00318 ENDIF 00319 CALL make_sparse_kinetic(preconditioner_env, matrix_t, matrix_s, energy_gap ,error=error) 00320 CALL make_diag_inner_precond(preconditioner_env, preconditioner_env%sparse_matrix,& 00321 preconditioner_env%sparse_matrix_inner, error=error) 00322 CASE DEFAULT 00323 CALL stop_program(routineN,moduleN,__LINE__,"Type not implemented") 00324 END SELECT 00325 ! 00326 ! here comes the solver 00327 SELECT CASE(my_solver_type) 00328 CASE (ot_precond_solver_inv_chol) 00329 ! 00330 ! compute the full inverse 00331 preconditioner_env%solver=ot_precond_solver_inv_chol 00332 CALL make_full_inverse_cholesky(preconditioner_env, preconditioner_env%fm, matrix_s, & 00333 & my_mixed_precision, error=error) 00334 CASE (ot_precond_solver_sainv) 00335 ! 00336 preconditioner_env%solver=ot_precond_solver_sainv 00337 CALL make_sparse_inverse_bif(preconditioner_env, preconditioner_env%fm, matrix_s, & 00338 & error=error) 00339 CASE (ot_precond_solver_direct) 00340 ! 00341 ! prepare for the direct solver 00342 preconditioner_env%solver=ot_precond_solver_direct 00343 CALL make_full_fact_cholesky(preconditioner_env, preconditioner_env%fm, matrix_s, & 00344 & error) 00345 CASE (ot_precond_solver_default) 00346 preconditioner_env%solver=ot_precond_solver_default 00347 CASE DEFAULT 00348 ! 00349 CALL stop_program(routineN,moduleN,__LINE__,"Doesn't know this type of solver") 00350 END SELECT 00351 00352 00353 IF(my_convert_precond_to_dbcsr) THEN 00354 IF(ASSOCIATED(preconditioner_env%dbcsr_matrix)) THEN 00355 CALL cp_dbcsr_release_p(preconditioner_env%dbcsr_matrix, error=error) 00356 ENDIF 00357 CALL cp_dbcsr_init_p(preconditioner_env%dbcsr_matrix,error=error) 00358 IF(my_mixed_precision) THEN 00359 CALL cp_dbcsr_create(preconditioner_env%dbcsr_matrix, "preconditioner_env%dbcsr_matrix", & 00360 cp_dbcsr_distribution(matrix_h), dbcsr_type_no_symmetry,& 00361 cp_dbcsr_row_block_sizes(matrix_h), cp_dbcsr_col_block_sizes(matrix_h), & 00362 0, 0, dbcsr_type_real_4, error=error) 00363 ELSE 00364 CALL cp_dbcsr_create(preconditioner_env%dbcsr_matrix, "preconditioner_env%dbcsr_matrix", & 00365 cp_dbcsr_distribution(matrix_h), dbcsr_type_no_symmetry,& 00366 cp_dbcsr_row_block_sizes(matrix_h), cp_dbcsr_col_block_sizes(matrix_h), & 00367 0, 0, dbcsr_type_real_default, error=error) 00368 ENDIF 00369 IF(ASSOCIATED(preconditioner_env%fm)) THEN 00370 CALL copy_fm_to_dbcsr(preconditioner_env%fm,preconditioner_env%dbcsr_matrix, error=error) 00371 ELSEIF(ASSOCIATED(preconditioner_env%sparse_matrix)) THEN 00372 CALL cp_dbcsr_copy(preconditioner_env%dbcsr_matrix, preconditioner_env%sparse_matrix, error=error) 00373 ELSE 00374 CALL stop_program(routineN,moduleN,__LINE__,"Something is wrong") 00375 ENDIF 00376 ENDIF 00377 00378 00379 IF (needs_full_spectrum) THEN 00380 DEALLOCATE(eigenvalues_ot,STAT=istat) 00381 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 00382 ENDIF 00383 00384 CALL timestop(handle) 00385 00386 END SUBROUTINE make_preconditioner 00387 00388 00389 ! applies a previously created preconditioner to a full matrix 00390 ! ***************************************************************************** 00391 SUBROUTINE apply_preconditioner_fm(preconditioner_env, matrix_in, matrix_out, error) 00392 00393 TYPE(preconditioner_type) :: preconditioner_env 00394 TYPE(cp_fm_type), POINTER :: matrix_in, matrix_out 00395 TYPE(cp_error_type), INTENT(inout) :: error 00396 00397 CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_fm', 00398 routineP = moduleN//':'//routineN 00399 00400 INTEGER :: handle 00401 00402 CALL timeset(routineN,handle) 00403 00404 SELECT CASE (preconditioner_env%in_use) 00405 CASE (0) 00406 CALL stop_program(routineN,moduleN,__LINE__,"No preconditioner in use") 00407 CASE (ot_precond_sparse_diag) 00408 CALL apply_sparse_diag(preconditioner_env, matrix_in, matrix_out,error=error) 00409 CASE (ot_precond_full_single) 00410 CALL apply_full_single(preconditioner_env, matrix_in, matrix_out,error=error) 00411 CASE (ot_precond_full_all) 00412 CALL apply_full_all(preconditioner_env, matrix_in, matrix_out,error=error) 00413 CASE(ot_precond_sparse_kinetic) 00414 CALL apply_solve_lin_system(preconditioner_env, matrix_in, matrix_out, error=error) 00415 CASE (ot_precond_sparse_kinetic_sainv) 00416 SELECT CASE (preconditioner_env%solver) 00417 CASE (ot_precond_solver_sainv) 00418 CALL apply_sparse_single(preconditioner_env, matrix_in, matrix_out,error=error) 00419 !CALL apply_full_single(preconditioner_env, matrix_in, matrix_out,error=error) 00420 CASE DEFAULT 00421 CALL stop_program(routineN,moduleN,__LINE__,"Solver not implemented") 00422 END SELECT 00423 CASE(ot_precond_full_kinetic,ot_precond_full_single_inverse,ot_precond_s_inverse) 00424 SELECT CASE (preconditioner_env%solver) 00425 CASE(ot_precond_solver_inv_chol,ot_precond_solver_sainv) 00426 CALL apply_full_single(preconditioner_env, matrix_in, matrix_out,error=error) 00427 CASE(ot_precond_solver_direct) 00428 CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out,error=error) 00429 CASE DEFAULT 00430 CALL stop_program(routineN,moduleN,__LINE__,"Solver not implemented") 00431 END SELECT 00432 CASE DEFAULT 00433 CALL stop_program(routineN,moduleN,__LINE__,"Unknown preconditioner") 00434 END SELECT 00435 00436 CALL timestop(handle) 00437 00438 END SUBROUTINE apply_preconditioner_fm 00439 00440 SUBROUTINE apply_preconditioner_dbcsr(preconditioner_env, matrix_in, matrix_out, error) 00441 00442 TYPE(preconditioner_type) :: preconditioner_env 00443 TYPE(cp_dbcsr_type) :: matrix_in, matrix_out 00444 TYPE(cp_error_type), INTENT(inout) :: error 00445 00446 CHARACTER(len=*), PARAMETER :: routineN = 'apply_preconditioner_dbcsr', 00447 routineP = moduleN//':'//routineN 00448 00449 INTEGER :: handle 00450 00451 CALL timeset(routineN,handle) 00452 00453 SELECT CASE (preconditioner_env%in_use) 00454 CASE (0) 00455 CALL stop_program(routineN,moduleN,__LINE__,"No preconditioner in use") 00456 CASE (ot_precond_sparse_diag,ot_precond_full_single) 00457 CALL apply_single(preconditioner_env, matrix_in, matrix_out,error=error) 00458 CASE (ot_precond_full_all) 00459 CALL apply_all(preconditioner_env, matrix_in, matrix_out,error=error) 00460 CASE(ot_precond_sparse_kinetic) 00461 CALL apply_solve_lin_system(preconditioner_env, matrix_in, matrix_out, error=error) 00462 CASE (ot_precond_sparse_kinetic_sainv) 00463 SELECT CASE (preconditioner_env%solver) 00464 CASE (ot_precond_solver_sainv) 00465 CALL apply_single(preconditioner_env, matrix_in, matrix_out,error=error) 00466 CASE DEFAULT 00467 CALL stop_program(routineN,moduleN,__LINE__,"Wrong solver") 00468 END SELECT 00469 CASE(ot_precond_full_kinetic,ot_precond_full_single_inverse,ot_precond_s_inverse) 00470 SELECT CASE (preconditioner_env%solver) 00471 CASE(ot_precond_solver_inv_chol,ot_precond_solver_sainv) 00472 CALL apply_single(preconditioner_env, matrix_in, matrix_out,error=error) 00473 CASE(ot_precond_solver_direct) 00474 CALL stop_program(routineN,moduleN,__LINE__,"Apply_full_direct not supported with ot") 00475 !CALL apply_full_direct(preconditioner_env, matrix_in, matrix_out,error=error) 00476 CASE DEFAULT 00477 CALL stop_program(routineN,moduleN,__LINE__,"Wrong solver") 00478 END SELECT 00479 CASE DEFAULT 00480 CALL stop_program(routineN,moduleN,__LINE__,"Wrong preconditioner") 00481 END SELECT 00482 00483 CALL timestop(handle) 00484 00485 END SUBROUTINE apply_preconditioner_dbcsr 00486 00487 ! different types of preconditioner come here 00488 ! a sparse block diagonal approximation 00489 ! ***************************************************************************** 00490 SUBROUTINE apply_sparse_diag(preconditioner_env, matrix_in, matrix_out, error) 00491 00492 TYPE(preconditioner_type) :: preconditioner_env 00493 TYPE(cp_fm_type), POINTER :: matrix_in, matrix_out 00494 TYPE(cp_error_type), INTENT(inout) :: error 00495 00496 INTEGER :: k 00497 00498 CALL cp_fm_get_info(matrix_in,ncol_global=k,error=error) 00499 CALL cp_dbcsr_sm_fm_multiply(preconditioner_env%sparse_matrix,matrix_in, & 00500 matrix_out,k,error=error) 00501 00502 END SUBROUTINE apply_sparse_diag 00503 00504 ! ***************************************************************************** 00505 SUBROUTINE make_sparse_diag(preconditioner_env, matrix_h, matrix_s, & 00506 energy_homo, energy_gap, error) 00507 00508 TYPE(preconditioner_type) :: preconditioner_env 00509 TYPE(cp_dbcsr_type), POINTER :: matrix_h, matrix_s 00510 REAL(KIND=dp) :: energy_homo, energy_gap 00511 TYPE(cp_error_type), INTENT(inout) :: error 00512 00513 CHARACTER(len=*), PARAMETER :: routineN = 'make_sparse_diag', 00514 routineP = moduleN//':'//routineN 00515 00516 INTEGER :: blk, handle, iblock_col, 00517 iblock_row, n, nblocks 00518 LOGICAL :: found 00519 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_h, block_pre, block_s 00520 TYPE(cp_dbcsr_iterator) :: iter 00521 00522 CALL timeset(routineN,handle) 00523 00524 IF (ASSOCIATED(preconditioner_env%sparse_matrix)) THEN 00525 CALL cp_dbcsr_deallocate_matrix(preconditioner_env%sparse_matrix,error=error) 00526 NULLIFY(preconditioner_env%sparse_matrix) 00527 ENDIF 00528 00529 ALLOCATE(preconditioner_env%sparse_matrix) 00530 CALL cp_dbcsr_init(preconditioner_env%sparse_matrix,error=error) 00531 CALL cp_dbcsr_create(preconditioner_env%sparse_matrix, ' PRECONDITIONER ', & 00532 cp_dbcsr_distribution (matrix_h),& 00533 cp_dbcsr_get_matrix_type (matrix_h), cp_dbcsr_row_block_sizes(matrix_h),& 00534 cp_dbcsr_col_block_sizes(matrix_h), 0, 0, error=error) 00535 00536 CALL cp_dbcsr_get_info(matrix_h,nfullrows_total=nblocks) 00537 00538 CALL cp_dbcsr_iterator_start(iter, matrix_h) 00539 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 00540 CALL cp_dbcsr_iterator_next_block(iter, iblock_row, iblock_col, block_h, blk) 00541 CALL cp_dbcsr_get_block_p(matrix=matrix_s,& 00542 row=iblock_row,col=iblock_col,BLOCK=block_s,found=found) 00543 00544 IF(.NOT.ASSOCIATED(block_s))CALL stop_program(routineN,moduleN,__LINE__,& 00545 ".NOT.ASSOCIATED(block_s)") 00546 00547 IF (iblock_col .EQ. iblock_row) THEN 00548 n=SIZE(block_s,1) 00549 00550 ALLOCATE(block_pre(n,n)) 00551 00552 CALL make_local_block(block_h,block_s,block_pre, & 00553 energy_homo,energy_gap,3) 00554 00555 CALL cp_dbcsr_put_block(matrix=preconditioner_env%sparse_matrix,& 00556 row=iblock_row,& 00557 col=iblock_col,& 00558 block=block_pre) 00559 00560 DEALLOCATE(block_pre) 00561 00562 ENDIF 00563 ENDDO 00564 00565 CALL cp_dbcsr_iterator_stop(iter) 00566 CALL cp_dbcsr_finalize(preconditioner_env%sparse_matrix,error=error) 00567 00568 CALL timestop(handle) 00569 00570 END SUBROUTINE make_sparse_diag 00571 ! ***************************************************************************** 00572 SUBROUTINE make_local_block(block_h,block_s,block_pre, & 00573 energy_homo,energy_gap, TYPE) 00574 00575 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_h, block_s, block_pre 00576 REAL(KIND=dp) :: energy_homo, energy_gap 00577 INTEGER :: TYPE 00578 00579 CHARACTER(len=*), PARAMETER :: routineN = 'make_local_block', 00580 routineP = moduleN//':'//routineN 00581 00582 INTEGER :: i, info, istat, liwork, 00583 lwork, n 00584 INTEGER, DIMENSION(:), POINTER :: iwork 00585 REAL(KIND=dp), DIMENSION(:), POINTER :: evals, work 00586 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_buf1, block_chol, 00587 block_evec 00588 00589 n=SIZE(block_s,1) 00590 lwork=1+6*n+2*n**2+50 00591 liwork=5*n+3 00592 ALLOCATE(block_chol(n,n)) 00593 ALLOCATE(block_evec(n,n)) 00594 ALLOCATE(block_buf1(n,n)) 00595 ALLOCATE(evals(n)) 00596 ALLOCATE(work(lwork),STAT=istat) 00597 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00598 "work",dp_size*lwork) 00599 ALLOCATE(iwork(liwork),STAT=istat) 00600 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 00601 "iwork",int_size*liwork) 00602 block_pre(:,:)=0.0_dp 00603 SELECT CASE (TYPE) 00604 CASE(1) 00605 DO i=1,n 00606 block_pre(i,i)=1.0_dp 00607 ENDDO 00608 CASE(2) 00609 DO i=1,n 00610 block_pre(i,i)=1.0_dp/MAX(energy_gap,block_h(i,i)-energy_homo) 00611 ENDDO 00612 CASE(3) 00613 ! more difficult constuct something like S^-0.5 K^T CASE(2) K^T S^-0.5 00614 block_chol(:,:)=block_s(:,:) 00615 block_evec(:,:)=block_h(:,:) 00616 CALL DPOTRF('U',n,block_chol(1,1),n,info) 00617 IF (info.ne.0) CALL stop_program(routineN,moduleN,__LINE__,"Error in dpotrf") 00618 CALL DSYGST(1,'U',n,block_evec(1,1),n,block_chol(1,1),n,info) 00619 IF (info.ne.0) CALL stop_program(routineN,moduleN,__LINE__,"Error in dsygst") 00620 CALL DSYEVD('V','U', n, block_evec(1,1), n, evals(1), work(1), lwork, & 00621 iwork(1), liwork, info) 00622 IF (info.NE.0) CALL stop_program(routineN,moduleN,__LINE__,"Error in dsyevd") 00623 block_pre(:,:)=0.0_dp 00624 DO i=1,n 00625 block_pre(i,i)=1.0_dp/MAX(evals(i)-energy_homo,energy_gap) 00626 ENDDO 00627 ! K = V E V ^ T 00628 CALL DGEMM('N','N',n,n,n,1.0_dp,block_evec(1,1),n,block_pre(1,1),n, & 00629 0.0_dp,block_buf1(1,1),n) 00630 CALL DGEMM('N','T',n,n,n,1.0_dp,block_buf1(1,1),n,block_evec(1,1),n, & 00631 0.0_dp,block_pre(1,1),n) 00632 ! inv(U) K inv(U)^T 00633 CALL DTRSM('L','U','N','N',n,n,1.0_dp,block_chol(1,1),n,block_pre(1,1),n) 00634 CALL DTRSM('R','U','T','N',n,n,1.0_dp,block_chol(1,1),n,block_pre(1,1),n) 00635 CASE(4) 00636 block_chol(:,:)=block_s(:,:) 00637 CALL DPOTRF('U',n,block_chol(1,1),n,info) 00638 block_pre(:,:)=0.0_dp 00639 DO i=1,n 00640 block_pre(i,i)=1.0_dp 00641 ENDDO 00642 CALL DTRSM('L','U','N','N',n,n,1.0_dp,block_chol(1,1),n,block_pre(1,1),n) 00643 CALL DTRSM('R','U','T','N',n,n,1.0_dp,block_chol(1,1),n,block_pre(1,1),n) 00644 CASE(5) ! like 3 but using s^-0.5 instead of the cholesky decomposition, and not transforming back 00645 block_evec(:,:)=block_s(:,:) 00646 CALL DSYEVD('V','U', n, block_evec(1,1), n, evals(1), work(1), lwork, & 00647 iwork(1), liwork, info) 00648 IF (info.ne.0) CALL stop_program(routineN,moduleN,__LINE__,"Error in dsyevd S") 00649 block_pre(:,:)=0.0_dp 00650 DO i=1,n 00651 block_pre(i,i)=1.0_dp/SQRT(evals(i)) 00652 ENDDO 00653 ! block_pre is s^-0.5 00654 CALL DGEMM('N','N',n,n,n,1.0_dp,block_evec(1,1),n,block_pre(1,1),n, & 00655 0.0_dp,block_buf1(1,1),n) 00656 CALL DGEMM('N','T',n,n,n,1.0_dp,block_buf1(1,1),n,block_evec(1,1),n, & 00657 0.0_dp,block_pre(1,1),n) 00658 ! transform H 00659 block_evec(:,:)=block_h(:,:) 00660 CALL DGEMM('N','N',n,n,n,1.0_dp,block_evec(1,1),n,block_pre(1,1),n, & 00661 0.0_dp,block_buf1(1,1),n) 00662 CALL DGEMM('N','N',n,n,n,1.0_dp,block_pre(1,1),n,block_buf1(1,1),n, & 00663 0.0_dp,block_evec(1,1),n) 00664 ! get evals and evecs 00665 CALL DSYEVD('V','U', n, block_evec(1,1), n, evals(1), work(1), lwork, & 00666 iwork(1), liwork, info) 00667 IF (info.ne.0) CALL stop_program(routineN,moduleN,__LINE__,"Error in dsyevd H") 00668 block_pre(:,:)=0.0_dp 00669 DO i=1,n 00670 block_pre(i,i)=1.0_dp/MAX(evals(i)-energy_homo,energy_gap) 00671 ENDDO 00672 CALL DGEMM('N','N',n,n,n,1.0_dp,block_evec(1,1),n,block_pre(1,1),n, & 00673 0.0_dp,block_buf1(1,1),n) 00674 CALL DGEMM('N','T',n,n,n,1.0_dp,block_buf1(1,1),n,block_evec(1,1),n, & 00675 0.0_dp,block_pre(1,1),n) 00676 CASE(6) ! like 3 not doing any transformation with s before or after (supposedly done by the caller) 00677 block_evec(:,:)=block_h(:,:) 00678 ! get evals and evecs 00679 CALL DSYEVD('V','U', n, block_evec(1,1), n, evals(1), work(1), lwork, & 00680 iwork(1), liwork, info) 00681 IF (info.ne.0) CALL stop_program(routineN,moduleN,__LINE__,"Error in dsyevd H") 00682 block_pre(:,:)=0.0_dp 00683 DO i=1,n 00684 block_pre(i,i)=1.0_dp/MAX(evals(i),energy_gap) 00685 ENDDO 00686 CALL DGEMM('N','N',n,n,n,1.0_dp,block_evec(1,1),n,block_pre(1,1),n, & 00687 0.0_dp,block_buf1(1,1),n) 00688 CALL DGEMM('N','T',n,n,n,1.0_dp,block_buf1(1,1),n,block_evec(1,1),n, & 00689 0.0_dp,block_pre(1,1),n) 00690 00691 END SELECT 00692 00693 DEALLOCATE (iwork,STAT=istat) 00694 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"iwork") 00695 DEALLOCATE (work,STAT=istat) 00696 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"work") 00697 DEALLOCATE (block_chol,STAT=istat) 00698 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"block_chol") 00699 DEALLOCATE (block_evec,STAT=istat) 00700 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"block_evec") 00701 DEALLOCATE (block_buf1,STAT=istat) 00702 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"block_buf1") 00703 DEALLOCATE (evals,STAT=istat) 00704 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,"evals") 00705 00706 END SUBROUTINE make_local_block 00707 ! different types of preconditioner come here 00708 ! a full matrix preconditioner 00709 ! ***************************************************************************** 00710 SUBROUTINE apply_full_single(preconditioner_env, matrix_in, matrix_out,error) 00711 00712 TYPE(preconditioner_type) :: preconditioner_env 00713 TYPE(cp_fm_type), POINTER :: matrix_in, matrix_out 00714 TYPE(cp_error_type), INTENT(inout) :: error 00715 00716 CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_single', 00717 routineP = moduleN//':'//routineN 00718 00719 INTEGER :: handle, k, n 00720 00721 CALL timeset(routineN,handle) 00722 00723 CALL cp_fm_get_info(matrix_in,nrow_global=n,ncol_global=k,error=error) 00724 CALL cp_fm_gemm('N','N',n,k,n,1.0_dp,preconditioner_env%fm, & 00725 matrix_in,0.0_dp,matrix_out,error=error) 00726 CALL timestop(handle) 00727 00728 END SUBROUTINE apply_full_single 00729 00730 SUBROUTINE apply_single(preconditioner_env, matrix_in, matrix_out,error) 00731 00732 TYPE(preconditioner_type) :: preconditioner_env 00733 TYPE(cp_dbcsr_type) :: matrix_in, matrix_out 00734 TYPE(cp_error_type), INTENT(inout) :: error 00735 00736 CHARACTER(len=*), PARAMETER :: routineN = 'apply_single', 00737 routineP = moduleN//':'//routineN 00738 00739 INTEGER :: handle, n 00740 00741 CALL timeset(routineN,handle) 00742 00743 IF (.NOT.ASSOCIATED(preconditioner_env%dbcsr_matrix)) & 00744 CALL stop_program(routineN,moduleN,__LINE__,& 00745 "NOT ASSOCIATED preconditioner_env%dbcsr_matrix") 00746 CALL cp_dbcsr_multiply('N','N',1.0_dp,preconditioner_env%dbcsr_matrix,matrix_in,& 00747 0.0_dp,matrix_out,& 00748 left_set=preconditioner_env%dbcsr_matrix%matrix%m%predistributed,error=error) 00749 00750 CALL timestop(handle) 00751 00752 END SUBROUTINE apply_single 00753 00754 SUBROUTINE apply_full_direct(preconditioner_env, matrix_in, matrix_out,error) 00755 00756 TYPE(preconditioner_type) :: preconditioner_env 00757 TYPE(cp_fm_type), POINTER :: matrix_in, matrix_out 00758 TYPE(cp_error_type), INTENT(inout) :: error 00759 00760 CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_direct', 00761 routineP = moduleN//':'//routineN 00762 00763 INTEGER :: handle, k, n 00764 TYPE(cp_fm_type), POINTER :: work 00765 00766 CALL timeset(routineN,handle) 00767 00768 CALL cp_fm_get_info(matrix_in,nrow_global=n,ncol_global=k,error=error) 00769 CALL cp_fm_create(work,matrix_in%matrix_struct,name="apply_full_single",& 00770 use_sp=matrix_in%use_sp,error=error) 00771 CALL cp_fm_cholesky_restore(matrix_in,k,preconditioner_env%fm,work,& 00772 & "SOLVE",transa="T",error=error) 00773 CALL cp_fm_cholesky_restore(work,k,preconditioner_env%fm,matrix_out,& 00774 & "SOLVE",transa="N",error=error) 00775 CALL cp_fm_release(work,error=error) 00776 00777 CALL timestop(handle) 00778 00779 END SUBROUTINE apply_full_direct 00780 00781 00782 ! ***************************************************************************** 00791 SUBROUTINE make_full_single_inverse(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, & 00792 energy_gap, mixed_precision, error) 00793 00794 TYPE(preconditioner_type) :: preconditioner_env 00795 TYPE(cp_fm_type), POINTER :: matrix_c0 00796 TYPE(cp_dbcsr_type), POINTER :: matrix_h, matrix_s 00797 REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals 00798 REAL(KIND=dp) :: energy_gap 00799 LOGICAL, INTENT(IN) :: mixed_precision 00800 TYPE(cp_error_type), INTENT(inout) :: error 00801 00802 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_inverse', 00803 routineP = moduleN//':'//routineN 00804 REAL(KIND=dp), PARAMETER :: eval_shift = 5.0_dp , 00805 fudge_factor = 2.0_dp 00806 00807 INTEGER :: handle, k, n 00808 REAL(KIND=dp) :: error_estimate, 00809 preconditioner_shift 00810 REAL(KIND=dp), DIMENSION(:), POINTER :: c0_shift, diag, evals 00811 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00812 TYPE(cp_fm_type), POINTER :: matrix_hc0, matrix_hc0_sp, matrix_ptr, 00813 matrix_s1, matrix_sc0, matrix_sc0_sp, matrix_shc0, matrix_tmp, 00814 matrix_tmp2 00815 00816 CALL timeset(routineN,handle) 00817 00818 ! arbitrary upshift of the occupied evals 00819 ! fudge factor for taking the error estimate into account 00820 00821 CALL cp_fm_get_info(matrix_c0,nrow_global=n,ncol_global=k,error=error) 00822 00823 IF (ASSOCIATED(preconditioner_env%fm)) CALL cp_fm_release(preconditioner_env%fm,error=error) 00824 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=n,ncol_global=n, & 00825 context=preconditioner_env%ctxt, & 00826 para_env=preconditioner_env%para_env,error=error) 00827 CALL cp_fm_create(preconditioner_env%fm,fm_struct_tmp,name="preconditioner_env%fm",& 00828 use_sp=mixed_precision,error=error) 00829 CALL cp_fm_create(matrix_tmp,fm_struct_tmp, name="preconditioner matrix_tmp",& 00830 use_sp=mixed_precision,error=error) 00831 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00832 00833 ! first try to get a ritz error estimate 00834 ! 0) cholesky decompose the overlap matrix, if this fails the basis is singular, 00835 ! more than EPS_DEFAULT 00836 CALL copy_dbcsr_to_fm(matrix_s,matrix_tmp,error=error) 00837 CALL cp_fm_cholesky_decompose(matrix_tmp,error=error) 00838 00839 ! get the error estimate 00840 CALL cp_fm_create(matrix_sc0,matrix_c0%matrix_struct,name="sc0",error=error) 00841 CALL cp_dbcsr_sm_fm_multiply(matrix_s,matrix_c0,matrix_sc0,k,error=error) 00842 CALL cp_fm_create(matrix_hc0,matrix_c0%matrix_struct,name="hc0",error=error) 00843 CALL cp_dbcsr_sm_fm_multiply(matrix_h,matrix_c0,matrix_hc0,k,error=error) 00844 00845 IF(mixed_precision) THEN 00846 CALL cp_fm_create(matrix_sc0_sp,matrix_c0%matrix_struct,name="sc0_sp",& 00847 use_sp=mixed_precision,error=error) 00848 CALL cp_fm_create(matrix_hc0_sp,matrix_c0%matrix_struct,name="hc0_sp",& 00849 use_sp=mixed_precision,error=error) 00850 CALL cp_fm_to_fm(matrix_sc0,matrix_sc0_sp,error=error) 00851 CALL cp_fm_to_fm(matrix_hc0,matrix_hc0_sp,error=error) 00852 matrix_ptr => matrix_sc0; matrix_sc0 => matrix_sc0_sp; matrix_sc0_sp => matrix_ptr 00853 CALL cp_fm_release(matrix_sc0_sp,error=error) 00854 matrix_ptr => matrix_hc0; matrix_hc0 => matrix_hc0_sp; matrix_hc0_sp => matrix_ptr 00855 CALL cp_fm_release(matrix_hc0_sp,error=error) 00856 ENDIF 00857 00858 CALL cp_fm_create(matrix_shc0,matrix_c0%matrix_struct,name="shc0",& 00859 use_sp=mixed_precision,error=error) 00860 CALL cp_fm_cholesky_restore(matrix_hc0,k,matrix_tmp,matrix_shc0,"SOLVE",transa="T",error=error) 00861 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=k,ncol_global=k, & 00862 context=preconditioner_env%ctxt, & 00863 para_env=preconditioner_env%para_env,error=error) 00864 CALL cp_fm_create(matrix_s1,fm_struct_tmp,name="matrix_s1",use_sp=mixed_precision,error=error) 00865 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00866 ! since we only use diagonal elements this is a bit of a waste 00867 CALL cp_fm_gemm('T','N',k,k,n,1.0_dp,matrix_shc0,matrix_shc0,0.0_dp,matrix_s1,error=error) 00868 ALLOCATE(diag(k)) 00869 CALL cp_fm_get_diag(matrix_s1,diag,error=error) 00870 error_estimate=MAXVAL(SQRT(ABS(diag-c0_evals**2))) 00871 DEALLOCATE(diag) 00872 CALL cp_fm_release(matrix_s1,error=error) 00873 CALL cp_fm_release(matrix_shc0,error=error) 00874 CALL cp_fm_release(matrix_hc0,error=error) 00875 00876 ! shift up the occupied subspace eigenvalues 00877 ALLOCATE(c0_shift(k)) 00878 c0_shift=SQRT(-(c0_evals-c0_evals(k))+eval_shift) 00879 CALL cp_fm_column_scale(matrix_sc0,c0_shift) 00880 CALL cp_fm_gemm('N','T',n,n,k,1.0_dp,matrix_sc0,matrix_sc0,0.0_dp,preconditioner_env%fm,error=error) 00881 CALL cp_fm_release(matrix_sc0,error=error) 00882 DEALLOCATE(c0_shift) 00883 00884 ! get H added to the shift 00885 CALL copy_dbcsr_to_fm(matrix_h,matrix_tmp,error=error) 00886 CALL cp_fm_scale_and_add(1.0_dp,preconditioner_env%fm,1.0_dp,matrix_tmp,error=error) 00887 00888 ! preconditioner shift, we target the middle of the occupied spectrum, and taking into account the error_estimate 00889 ! write(6,*) "Error estimate = ",error_estimate 00890 preconditioner_shift=-(MINVAL(c0_evals)+ MAXVAL(c0_evals))/2.0_dp + & 00891 error_estimate*fudge_factor 00892 CALL copy_dbcsr_to_fm(matrix_s,matrix_tmp,error=error) 00893 CALL cp_fm_scale_and_add(1.0_dp,preconditioner_env%fm,preconditioner_shift,matrix_tmp,error=error) 00894 00895 ! check evals 00896 IF (.FALSE.) THEN 00897 CALL cp_fm_to_fm(preconditioner_env%fm,matrix_tmp,error=error) 00898 CALL cp_fm_create(matrix_tmp2,matrix_tmp%matrix_struct,name="matrix_tmp2",error=error) 00899 ALLOCATE(evals(n)) 00900 CALL choose_eigv_solver(matrix_tmp,matrix_tmp2,evals,error=error) 00901 00902 CALL cp_fm_release(matrix_tmp2,error=error) 00903 WRITE(6,*) "evals ",evals 00904 DEALLOCATE(evals) 00905 ENDIF 00906 00907 CALL cp_fm_release(matrix_tmp,error=error) 00908 CALL timestop(handle) 00909 00910 END SUBROUTINE make_full_single_inverse 00911 00912 ! ***************************************************************************** 00926 SUBROUTINE make_full_all(preconditioner_env, matrix_c0, matrix_h, matrix_s, c0_evals, energy_gap, error) 00927 00928 TYPE(preconditioner_type) :: preconditioner_env 00929 TYPE(cp_fm_type), POINTER :: matrix_c0 00930 TYPE(cp_dbcsr_type), POINTER :: matrix_h, matrix_s 00931 REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals 00932 REAL(KIND=dp) :: energy_gap 00933 TYPE(cp_error_type), INTENT(inout) :: error 00934 00935 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all', 00936 routineP = moduleN//':'//routineN 00937 REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, 00938 lambda_base = 10.0_dp 00939 00940 INTEGER :: handle, k, n 00941 REAL(KIND=dp) :: error_estimate, lambda 00942 REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals 00943 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 00944 TYPE(cp_fm_type), POINTER :: matrix_hc0, matrix_left, matrix_pre, 00945 matrix_s1, matrix_s2, matrix_sc0, matrix_shc0, matrix_tmp, ortho 00946 00947 CALL timeset(routineN,handle) 00948 00949 IF (ASSOCIATED(preconditioner_env%fm)) CALL cp_fm_release(preconditioner_env%fm,error) 00950 CALL cp_fm_get_info(matrix_c0,nrow_global=n,ncol_global=k,error=error) 00951 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=n,ncol_global=n, & 00952 context=preconditioner_env%ctxt, & 00953 para_env=preconditioner_env%para_env,error=error) 00954 CALL cp_fm_create(preconditioner_env%fm,fm_struct_tmp,name="preconditioner_env%fm",error=error) 00955 matrix_pre=>preconditioner_env%fm 00956 CALL cp_fm_create(ortho,fm_struct_tmp,name="ortho",error=error) 00957 CALL cp_fm_create(matrix_tmp,fm_struct_tmp,name="matrix_tmp",error=error) 00958 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00959 ALLOCATE(preconditioner_env%full_evals(n)) 00960 ALLOCATE(preconditioner_env%occ_evals(k)) 00961 00962 ! 0) cholesky decompose the overlap matrix, if this fails the basis is singular, 00963 ! more than EPS_DEFAULT 00964 CALL copy_dbcsr_to_fm(matrix_s,ortho,error=error) 00965 CALL cp_fm_cholesky_decompose(ortho,error=error) 00966 00967 ! 1) Construct a new H matrix, which has the current C0 as eigenvectors, 00968 ! possibly shifted by an amount lambda, 00969 ! and the same spectrum as the original H matrix in the space orthogonal to the C0 00970 ! with P=C0 C0 ^ T 00971 ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS) 00972 ! we exploit that the C0 are already the ritz states of H 00973 CALL cp_fm_create(matrix_sc0,matrix_c0%matrix_struct,name="sc0",error=error) 00974 CALL cp_dbcsr_sm_fm_multiply(matrix_s,matrix_c0,matrix_sc0,k,error=error) 00975 CALL cp_fm_create(matrix_hc0,matrix_c0%matrix_struct,name="hc0",error=error) 00976 CALL cp_dbcsr_sm_fm_multiply(matrix_h,matrix_c0,matrix_hc0,k,error=error) 00977 00978 ! An aside, try to estimate the error on the ritz values, we'll need it later on 00979 CALL cp_fm_create(matrix_shc0,matrix_c0%matrix_struct,name="shc0",error=error) 00980 CALL cp_fm_cholesky_restore(matrix_hc0,k,ortho,matrix_shc0,"SOLVE",transa="T",error=error) 00981 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=k,ncol_global=k, & 00982 context=preconditioner_env%ctxt, & 00983 para_env=preconditioner_env%para_env,error=error) 00984 CALL cp_fm_create(matrix_s1,fm_struct_tmp,name="matrix_s1",error=error) 00985 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 00986 ! since we only use diagonal elements this is a bit of a waste 00987 CALL cp_fm_gemm('T','N',k,k,n,1.0_dp,matrix_shc0,matrix_shc0,0.0_dp,matrix_s1,error=error) 00988 ALLOCATE(diag(k)) 00989 CALL cp_fm_get_diag(matrix_s1,diag,error=error) 00990 error_estimate=MAXVAL(SQRT(ABS(diag-c0_evals**2))) 00991 DEALLOCATE(diag) 00992 CALL cp_fm_release(matrix_s1,error=error) 00993 CALL cp_fm_release(matrix_shc0,error=error) 00994 ! we'll only use the energy gap, if our estimate of the error on the eigenvalues 00995 ! is small enough. A large error combined with a small energy gap would otherwise lead to 00996 ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition 00997 ! aggressively 00998 preconditioner_env%energy_gap= MAX(energy_gap,error_estimate*fudge_factor) 00999 CALL copy_dbcsr_to_fm(matrix_h,matrix_tmp,error=error) 01000 CALL cp_fm_upper_to_full(matrix_tmp,matrix_pre,error=error) 01001 ! tmp = H ( 1 - PS ) 01002 CALL cp_fm_gemm('N','T',n,n,k,-1.0_dp,matrix_hc0,matrix_sc0,1.0_dp,matrix_tmp,error=error) 01003 01004 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=k,ncol_global=n, & 01005 context=preconditioner_env%ctxt, & 01006 para_env=preconditioner_env%para_env,error=error) 01007 CALL cp_fm_create(matrix_left,fm_struct_tmp,name="matrix_left",error=error) 01008 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01009 CALL cp_fm_gemm('T','N',k,n,n,1.0_dp,matrix_c0,matrix_tmp,0.0_dp,matrix_left,error=error) 01010 ! tmp = (1 - PS)^T H (1-PS) 01011 CALL cp_fm_gemm('N','N',n,n,k,-1.0_dp,matrix_sc0,matrix_left,1.0_dp,matrix_tmp,error=error) 01012 CALL cp_fm_release(matrix_left,error=error) 01013 01014 ALLOCATE(shifted_evals(k)) 01015 lambda = lambda_base + error_estimate 01016 shifted_evals=c0_evals - lambda 01017 CALL cp_fm_to_fm(matrix_sc0,matrix_hc0,error=error) 01018 CALL cp_fm_column_scale(matrix_hc0,shifted_evals) 01019 CALL cp_fm_gemm('N','T',n,n,k,1.0_dp,matrix_hc0,matrix_sc0,1.0_dp,matrix_tmp,error=error) 01020 01021 ! 2) diagonalize this operator 01022 CALL cp_fm_cholesky_reduce(matrix_tmp,ortho,error=error) 01023 CALL choose_eigv_solver(matrix_tmp,matrix_pre,preconditioner_env%full_evals,error=error) 01024 CALL cp_fm_cholesky_restore(matrix_pre,n,ortho,matrix_tmp,"SOLVE",error=error) 01025 CALL cp_fm_to_fm(matrix_tmp,matrix_pre,error=error) 01026 01027 ! test that the subspace remained conserved 01028 IF (.FALSE.) THEN 01029 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=k,ncol_global=k, & 01030 context=preconditioner_env%ctxt, & 01031 para_env=preconditioner_env%para_env,error=error) 01032 CALL cp_fm_create(matrix_s1,fm_struct_tmp,name="matrix_s1",error=error) 01033 CALL cp_fm_create(matrix_s2,fm_struct_tmp,name="matrix_s2",error=error) 01034 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01035 ALLOCATE(norms(k)) 01036 CALL cp_fm_gemm('T','N',k,k,n,1.0_dp,matrix_sc0,matrix_tmp,0.0_dp,matrix_s1,error=error) 01037 CALL choose_eigv_solver(matrix_s1,matrix_s2,norms,error=error) 01038 WRITE(6,*) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms)-1.0_dp)) 01039 DEALLOCATE(norms) 01040 CALL cp_fm_release(matrix_s1,error=error) 01041 CALL cp_fm_release(matrix_s2,error=error) 01042 ENDIF 01043 01044 ! 3) replace the lowest k evals and evecs with what they should be 01045 preconditioner_env%occ_evals=c0_evals 01046 ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all) 01047 preconditioner_env%full_evals(1:k)=c0_evals 01048 CALL cp_fm_to_fm(matrix_c0,matrix_pre,k,1,1) 01049 01050 CALL cp_fm_release(matrix_sc0,error=error) 01051 CALL cp_fm_release(matrix_hc0,error=error) 01052 CALL cp_fm_release(ortho,error=error) 01053 CALL cp_fm_release(matrix_tmp,error=error) 01054 DEALLOCATE(shifted_evals) 01055 CALL timestop(handle) 01056 01057 END SUBROUTINE make_full_all 01058 01059 ! 01060 ! the corresponding apply_full_all uses the decomposed form to apply the preconditioner 01061 ! 01062 01063 ! ***************************************************************************** 01064 SUBROUTINE apply_full_all(preconditioner_env, matrix_in, matrix_out, error) 01065 01066 TYPE(preconditioner_type) :: preconditioner_env 01067 TYPE(cp_fm_type), POINTER :: matrix_in, matrix_out 01068 TYPE(cp_error_type), INTENT(inout) :: error 01069 01070 CHARACTER(len=*), PARAMETER :: routineN = 'apply_full_all', 01071 routineP = moduleN//':'//routineN 01072 01073 INTEGER :: handle, i, j, k, n, 01074 ncol_local, nrow_local 01075 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices 01076 REAL(KIND=dp) :: dum 01077 REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_data 01078 TYPE(cp_fm_type), POINTER :: matrix_tmp 01079 01080 CALL timeset(routineN,handle) 01081 01082 CALL cp_fm_get_info(matrix_in,nrow_global=n,ncol_global=k,error=error) 01083 01084 CALL cp_fm_create(matrix_tmp,matrix_in%matrix_struct,name="apply_full_all",error=error) 01085 CALL cp_fm_get_info(matrix_tmp, nrow_local=nrow_local, ncol_local=ncol_local, & 01086 row_indices=row_indices, col_indices=col_indices, local_data=local_data,error=error) 01087 01088 ! 01089 CALL cp_fm_gemm('T','N',n,k,n,1.0_dp,preconditioner_env%fm, & 01090 matrix_in,0.0_dp,matrix_tmp,error=error) 01091 01092 ! do the right scaling 01093 DO j=1,ncol_local 01094 DO i=1,nrow_local 01095 dum=1.0_dp/MAX(preconditioner_env%energy_gap, & 01096 preconditioner_env%full_evals(row_indices(i))-preconditioner_env%occ_evals(col_indices(j))) 01097 local_data(i,j)=local_data(i,j)*dum 01098 ENDDO 01099 ENDDO 01100 01101 ! mult back 01102 CALL cp_fm_gemm('N','N',n,k,n,1.0_dp,preconditioner_env%fm, & 01103 matrix_tmp,0.0_dp,matrix_out,error=error) 01104 01105 CALL cp_fm_release(matrix_tmp,error=error) 01106 01107 CALL timestop(handle) 01108 01109 END SUBROUTINE apply_full_all 01110 01111 SUBROUTINE apply_all(preconditioner_env, matrix_in, matrix_out, error) 01112 01113 TYPE(preconditioner_type) :: preconditioner_env 01114 TYPE(cp_dbcsr_type) :: matrix_in, matrix_out 01115 TYPE(cp_error_type), INTENT(inout) :: error 01116 01117 CHARACTER(len=*), PARAMETER :: routineN = 'apply_all', 01118 routineP = moduleN//':'//routineN 01119 01120 INTEGER :: col, col_offset, col_size, 01121 handle, i, j, row, 01122 row_offset, row_size 01123 REAL(KIND=dp) :: dum 01124 REAL(KIND=dp), DIMENSION(:, :), POINTER :: DATA 01125 TYPE(cp_dbcsr_iterator) :: iter 01126 TYPE(cp_dbcsr_type) :: matrix_tmp 01127 01128 CALL timeset(routineN,handle) 01129 01130 CALL cp_dbcsr_init(matrix_tmp,error=error) 01131 CALL cp_dbcsr_copy(matrix_tmp,matrix_in,name="apply_full_all",error=error) 01132 CALL cp_dbcsr_multiply('T','N',1.0_dp,preconditioner_env%dbcsr_matrix, & 01133 matrix_in,0.0_dp,matrix_tmp,error=error) 01134 ! do the right scaling 01135 CALL cp_dbcsr_iterator_start(iter, matrix_tmp) 01136 DO WHILE (cp_dbcsr_iterator_blocks_left (iter)) 01137 CALL cp_dbcsr_iterator_next_block(iter, row, col, DATA, & 01138 row_size=row_size, col_size=col_size, & 01139 row_offset=row_offset, col_offset=col_offset) 01140 DO j=1,col_size 01141 DO i=1,row_size 01142 dum=1.0_dp/MAX(preconditioner_env%energy_gap, & 01143 preconditioner_env%full_evals( row_offset+i-1 )& 01144 -preconditioner_env%occ_evals( col_offset+j-1 )) 01145 DATA(i,j)=DATA(i,j)*dum 01146 ENDDO 01147 ENDDO 01148 ENDDO 01149 CALL cp_dbcsr_iterator_stop(iter) 01150 ! mult back 01151 CALL cp_dbcsr_multiply('N','N',1.0_dp,preconditioner_env%dbcsr_matrix, & 01152 matrix_tmp,0.0_dp,matrix_out,error=error) 01153 CALL cp_dbcsr_release(matrix_tmp, error=error) 01154 CALL timestop(handle) 01155 01156 END SUBROUTINE apply_all 01157 01158 ! ***************************************************************************** 01159 SUBROUTINE make_full_single(preconditioner_env, fm, matrix_h, matrix_s, & 01160 energy_homo, energy_gap , error) 01161 01162 TYPE(preconditioner_type) :: preconditioner_env 01163 TYPE(cp_fm_type), POINTER :: fm 01164 TYPE(cp_dbcsr_type), POINTER :: matrix_h, matrix_s 01165 REAL(KIND=dp) :: energy_homo, energy_gap 01166 TYPE(cp_error_type), INTENT(inout) :: error 01167 01168 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single', 01169 routineP = moduleN//':'//routineN 01170 01171 INTEGER :: handle, i, n 01172 REAL(KIND=dp), DIMENSION(:), POINTER :: evals 01173 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 01174 TYPE(cp_fm_type), POINTER :: fm_h, fm_s 01175 01176 ! --- 01177 ! --- 01178 01179 CALL timeset(routineN,handle) 01180 01181 NULLIFY(fm_h,fm_s,fm_struct_tmp,evals) 01182 01183 IF (ASSOCIATED(fm)) THEN 01184 CALL cp_fm_release(fm,error=error) 01185 ENDIF 01186 CALL cp_dbcsr_get_info(matrix_h,nfullrows_total=n) 01187 ALLOCATE(evals(n)) 01188 01189 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n,ncol_global=n,& 01190 context=preconditioner_env%ctxt, & 01191 para_env=preconditioner_env%para_env,error=error) 01192 CALL cp_fm_create(fm,fm_struct_tmp, name="preconditioner",error=error) 01193 CALL cp_fm_create(fm_h,fm_struct_tmp, name="fm_h",error=error) 01194 CALL cp_fm_create(fm_s,fm_struct_tmp, name="fm_s",error=error) 01195 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01196 01197 CALL copy_dbcsr_to_fm(matrix_h,fm_h,error=error) 01198 CALL copy_dbcsr_to_fm(matrix_s,fm_s,error=error) 01199 CALL cp_fm_cholesky_decompose(fm_s,error=error) 01200 CALL cp_fm_cholesky_reduce(fm_h,fm_s,error=error) 01201 CALL choose_eigv_solver(fm_h,fm,evals,error=error) 01202 CALL cp_fm_cholesky_restore(fm,n,fm_s,fm_h,"SOLVE",error=error) 01203 DO i=1,n 01204 evals(i)=1.0_dp/MAX(evals(i)-energy_homo,energy_gap) 01205 ENDDO 01206 CALL cp_fm_to_fm(fm_h,fm,error=error) 01207 CALL cp_fm_column_scale(fm,evals) 01208 CALL cp_fm_gemm('N','T',n,n,n,1.0_dp,fm,fm_h,0.0_dp,fm_s,error=error) 01209 CALL cp_fm_to_fm(fm_s,fm,error=error) 01210 01211 DEALLOCATE(evals) 01212 CALL cp_fm_release(fm_h,error=error) 01213 CALL cp_fm_release(fm_s,error=error) 01214 01215 CALL timestop(handle) 01216 01217 END SUBROUTINE make_full_single 01218 01219 ! different types of preconditioner come here 01220 ! ***************************************************************************** 01221 SUBROUTINE make_full_s_inverse(preconditioner_env, matrix_h, matrix_s, error) 01222 01223 TYPE(preconditioner_type) :: preconditioner_env 01224 TYPE(cp_dbcsr_type), POINTER :: matrix_h, matrix_s 01225 TYPE(cp_error_type), INTENT(inout) :: error 01226 01227 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_s_inverse', 01228 routineP = moduleN//':'//routineN 01229 01230 INTEGER :: handle, n 01231 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 01232 01233 CALL timeset(routineN,handle) 01234 01235 NULLIFY(fm_struct_tmp) 01236 01237 IF (ASSOCIATED(preconditioner_env%fm)) THEN 01238 CALL cp_fm_release(preconditioner_env%fm,error=error) 01239 ENDIF 01240 CALL cp_dbcsr_get_info(matrix_h,nfullrows_total=n) 01241 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n,ncol_global=n,& 01242 context=preconditioner_env%ctxt, & 01243 para_env=preconditioner_env%para_env,error=error) 01244 CALL cp_fm_create(preconditioner_env%fm,fm_struct_tmp, name="preconditioner",error=error) 01245 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01246 01247 CALL copy_dbcsr_to_fm(matrix_s,preconditioner_env%fm,error=error) 01248 01249 CALL timestop(handle) 01250 01251 END SUBROUTINE make_full_s_inverse 01252 01253 ! ***************************************************************************** 01254 SUBROUTINE make_full_kinetic(preconditioner_env, fm, matrix_t, matrix_s, & 01255 energy_gap, mixed_precision, error) 01256 TYPE(preconditioner_type) :: preconditioner_env 01257 TYPE(cp_fm_type), POINTER :: fm 01258 TYPE(cp_dbcsr_type), POINTER :: matrix_t, matrix_s 01259 REAL(KIND=dp) :: energy_gap 01260 LOGICAL, INTENT(IN) :: mixed_precision 01261 TYPE(cp_error_type), INTENT(inout) :: error 01262 01263 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_kinetic', 01264 routineP = moduleN//':'//routineN 01265 01266 INTEGER :: handle, n 01267 LOGICAL :: failure 01268 REAL(KIND=dp) :: shift 01269 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 01270 01271 failure = .FALSE. 01272 CALL timeset(routineN,handle) 01273 01274 CPPrecondition(ASSOCIATED(matrix_t),cp_failure_level,routineP,error,failure) 01275 CPPrecondition(ASSOCIATED(matrix_s),cp_failure_level,routineP,error,failure) 01276 01277 NULLIFY(fm_struct_tmp) 01278 01279 IF (ASSOCIATED(fm)) THEN 01280 CALL cp_fm_release(fm,error=error) 01281 ENDIF 01282 CALL cp_dbcsr_get_info(matrix_t,nfullrows_total=n) 01283 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n,ncol_global=n,& 01284 context=preconditioner_env%ctxt, & 01285 para_env=preconditioner_env%para_env,error=error) 01286 CALL cp_fm_create(fm,fm_struct_tmp, name="preconditioner", use_sp=mixed_precision, error=error) 01287 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01288 01289 shift=MAX(0.0_dp,energy_gap) 01290 CALL cp_dbcsr_add(matrix_t,matrix_s,alpha_scalar=1.0_dp,beta_scalar=shift,error=error) 01291 CALL copy_dbcsr_to_fm(matrix_t,fm,error=error) 01292 CALL timestop(handle) 01293 01294 END SUBROUTINE make_full_kinetic 01295 01296 SUBROUTINE make_full_inverse_cholesky(preconditioner_env, fm, matrix_s, mixed_precision, & 01297 error) 01298 01299 TYPE(preconditioner_type) :: preconditioner_env 01300 TYPE(cp_fm_type), POINTER :: fm 01301 TYPE(cp_dbcsr_type), OPTIONAL, POINTER :: matrix_s 01302 LOGICAL, INTENT(IN) :: mixed_precision 01303 TYPE(cp_error_type), INTENT(inout) :: error 01304 01305 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_inverse_cholesky', 01306 routineP = moduleN//':'//routineN 01307 01308 INTEGER :: handle 01309 LOGICAL :: failure 01310 TYPE(cp_error_type) :: sub_error 01311 TYPE(cp_fm_type), POINTER :: fm_work 01312 01313 failure = .FALSE. 01314 01315 CALL timeset(routineN,handle) 01316 01317 CPPrecondition(ASSOCIATED(fm),cp_failure_level,routineP,error,failure) 01318 01319 NULLIFY(fm_work) 01320 01321 CALL cp_fm_create(fm_work,fm%matrix_struct,name="fm_work",use_sp=mixed_precision,error=error) 01322 ! 01323 ! compute the inverse of SPD matrix fm using the Cholesky factorization 01324 CALL cp_error_init(sub_error,template_error=error,stop_level=cp_fatal_level) 01325 01326 CALL cp_fm_cholesky_decompose(fm,error=sub_error) 01327 01328 failure = .FALSE. 01329 CALL cp_error_check(sub_error,failure) 01330 CALL cp_error_dealloc_ref(sub_error) 01331 ! 01332 ! if fm not SPD we go with the overlap matrix 01333 IF (failure) THEN 01334 ! 01335 ! just the overlap matrix 01336 IF(PRESENT(matrix_s)) THEN 01337 CALL copy_dbcsr_to_fm(matrix_s,fm,error=error) 01338 CALL cp_fm_cholesky_decompose(fm,error=error) 01339 ELSE 01340 CALL cp_fm_set_all(fm,alpha=0._dp,beta=1._dp,error=error) 01341 ENDIF 01342 ENDIF 01343 CALL cp_fm_cholesky_invert(fm,error=error) 01344 01345 CALL cp_fm_upper_to_full(fm,fm_work,error=error) 01346 CALL cp_fm_release(fm_work,error=error) 01347 01348 CALL timestop(handle) 01349 01350 END SUBROUTINE make_full_inverse_cholesky 01351 01352 SUBROUTINE make_full_fact_cholesky(preconditioner_env, fm, matrix_s, error) 01353 01354 TYPE(preconditioner_type) :: preconditioner_env 01355 TYPE(cp_fm_type), POINTER :: fm 01356 TYPE(cp_dbcsr_type), OPTIONAL, POINTER :: matrix_s 01357 TYPE(cp_error_type), INTENT(inout) :: error 01358 01359 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_fact_cholesky', 01360 routineP = moduleN//':'//routineN 01361 01362 INTEGER :: handle 01363 LOGICAL :: failure 01364 TYPE(cp_error_type) :: sub_error 01365 01366 failure = .FALSE. 01367 01368 CALL timeset(routineN,handle) 01369 01370 CPPrecondition(ASSOCIATED(fm),cp_failure_level,routineP,error,failure) 01371 ! 01372 ! compute the inverse of SPD matrix fm using the Cholesky factorization 01373 CALL cp_error_init(sub_error,template_error=error,stop_level=cp_fatal_level) 01374 CALL cp_fm_cholesky_decompose(fm,error=sub_error) 01375 failure = .FALSE. 01376 CALL cp_error_check(sub_error,failure) 01377 CALL cp_error_dealloc_ref(sub_error) 01378 ! 01379 ! if fm not SPD we go with the overlap matrix 01380 IF (failure) THEN 01381 ! 01382 ! just the overlap matrix 01383 IF(PRESENT(matrix_s)) THEN 01384 CALL copy_dbcsr_to_fm(matrix_s,fm,error=error) 01385 CALL cp_fm_cholesky_decompose(fm,error=error) 01386 ELSE 01387 CALL cp_fm_set_all(fm,alpha=0._dp,beta=1._dp,error=error) 01388 ENDIF 01389 ENDIF 01390 01391 CALL timestop(handle) 01392 01393 END SUBROUTINE make_full_fact_cholesky 01394 01395 ! ***************************************************************************** 01396 ! Preconditioners for orthogonal basis sets 01397 ! ***************************************************************************** 01398 SUBROUTINE make_sparse_diag_ortho(preconditioner_env, matrix_h, energy_homo, energy_gap, error) 01399 01400 TYPE(preconditioner_type) :: preconditioner_env 01401 TYPE(cp_dbcsr_type), POINTER :: matrix_h 01402 REAL(KIND=dp) :: energy_homo, energy_gap 01403 TYPE(cp_error_type), INTENT(inout) :: error 01404 01405 CHARACTER(len=*), PARAMETER :: routineN = 'make_sparse_diag_ortho', 01406 routineP = moduleN//':'//routineN 01407 01408 INTEGER :: blk, handle, i, iblock_col, 01409 iblock_row, n, nblocks 01410 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_h, block_pre, block_s 01411 TYPE(cp_dbcsr_iterator) :: iter 01412 01413 CALL timeset(routineN,handle) 01414 01415 IF (ASSOCIATED(preconditioner_env%sparse_matrix)) THEN 01416 CALL cp_dbcsr_deallocate_matrix(preconditioner_env%sparse_matrix,error=error) 01417 NULLIFY(preconditioner_env%sparse_matrix) 01418 ENDIF 01419 01420 ALLOCATE(preconditioner_env%sparse_matrix) 01421 CALL cp_dbcsr_init(preconditioner_env%sparse_matrix,error=error) 01422 CALL cp_dbcsr_create(preconditioner_env%sparse_matrix, ' PRECONDITIONER ', & 01423 cp_dbcsr_distribution (matrix_h), cp_dbcsr_get_matrix_type (matrix_h),& 01424 cp_dbcsr_row_block_sizes(matrix_h),& 01425 cp_dbcsr_col_block_sizes(matrix_h), 0, 0, error=error) 01426 01427 CALL cp_dbcsr_get_info(matrix_h,nfullrows_total=nblocks) 01428 01429 CALL cp_dbcsr_iterator_start(iter, matrix_h) 01430 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 01431 CALL cp_dbcsr_iterator_next_block(iter, iblock_row, iblock_col, block_h, blk) 01432 01433 IF (iblock_col .EQ. iblock_row) THEN 01434 n=SIZE(block_h,1) 01435 ALLOCATE(block_pre(n,n),block_s(n,n)) 01436 block_s(:,:)=0._dp 01437 DO i=1,n 01438 block_s(i,i)=1._dp 01439 END DO 01440 01441 CALL make_local_block(block_h,block_s,block_pre,energy_homo,energy_gap,3) 01442 01443 CALL cp_dbcsr_put_block(matrix=preconditioner_env%sparse_matrix,& 01444 row=iblock_row,& 01445 col=iblock_col,& 01446 BLOCK=block_pre) 01447 01448 DEALLOCATE(block_pre,block_s) 01449 01450 ENDIF 01451 ENDDO 01452 CALL cp_dbcsr_iterator_stop(iter) 01453 CALL cp_dbcsr_finalize(preconditioner_env%sparse_matrix,error=error) 01454 01455 CALL timestop(handle) 01456 01457 END SUBROUTINE make_sparse_diag_ortho 01458 01459 ! ***************************************************************************** 01460 01461 SUBROUTINE make_full_single_ortho(preconditioner_env, fm, matrix_h, & 01462 energy_homo, energy_gap , error) 01463 01464 TYPE(preconditioner_type) :: preconditioner_env 01465 TYPE(cp_fm_type), POINTER :: fm 01466 TYPE(cp_dbcsr_type), POINTER :: matrix_h 01467 REAL(KIND=dp) :: energy_homo, energy_gap 01468 TYPE(cp_error_type), INTENT(inout) :: error 01469 01470 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_single_ortho', 01471 routineP = moduleN//':'//routineN 01472 01473 INTEGER :: handle, i, n 01474 REAL(KIND=dp), DIMENSION(:), POINTER :: evals 01475 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 01476 TYPE(cp_fm_type), POINTER :: fm_h, fm_s 01477 01478 CALL timeset(routineN,handle) 01479 NULLIFY(fm_h,fm_s,fm_struct_tmp,evals) 01480 01481 IF (ASSOCIATED(fm)) THEN 01482 CALL cp_fm_release(fm,error=error) 01483 ENDIF 01484 CALL cp_dbcsr_get_info(matrix_h,nfullrows_total=n) 01485 ALLOCATE(evals(n)) 01486 01487 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n,ncol_global=n,& 01488 context=preconditioner_env%ctxt, & 01489 para_env=preconditioner_env%para_env,error=error) 01490 CALL cp_fm_create(fm,fm_struct_tmp, name="preconditioner",error=error) 01491 CALL cp_fm_create(fm_h,fm_struct_tmp, name="fm_h",error=error) 01492 CALL cp_fm_create(fm_s,fm_struct_tmp, name="fm_s",error=error) 01493 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01494 01495 CALL copy_dbcsr_to_fm(matrix_h,fm_h,error=error) 01496 01497 CALL choose_eigv_solver(fm_h,fm,evals,error=error) 01498 DO i=1,n 01499 evals(i)=1.0_dp/MAX(evals(i)-energy_homo,energy_gap) 01500 ENDDO 01501 CALL cp_fm_to_fm(fm,fm_h,error=error) 01502 CALL cp_fm_column_scale(fm,evals) 01503 CALL cp_fm_gemm('N','T',n,n,n,1.0_dp,fm,fm_h,0.0_dp,fm_s,error=error) 01504 CALL cp_fm_to_fm(fm_s,fm,error=error) 01505 01506 DEALLOCATE(evals) 01507 CALL cp_fm_release(fm_h,error=error) 01508 CALL cp_fm_release(fm_s,error=error) 01509 01510 CALL timestop(handle) 01511 01512 END SUBROUTINE make_full_single_ortho 01513 ! ***************************************************************************** 01514 SUBROUTINE make_full_single_inverse_ortho(preconditioner_env, matrix_c0, matrix_h, c0_evals, energy_gap,error) 01515 01516 TYPE(preconditioner_type) :: preconditioner_env 01517 TYPE(cp_fm_type), POINTER :: matrix_c0 01518 TYPE(cp_dbcsr_type), POINTER :: matrix_h 01519 REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals 01520 REAL(KIND=dp) :: energy_gap 01521 TYPE(cp_error_type), INTENT(inout) :: error 01522 01523 CHARACTER(len=*), PARAMETER :: 01524 routineN = 'make_full_single_inverse_ortho', 01525 routineP = moduleN//':'//routineN 01526 REAL(KIND=dp), PARAMETER :: eval_shift = 5.0_dp , 01527 fudge_factor = 2.0_dp 01528 01529 INTEGER :: handle, k, n 01530 REAL(KIND=dp) :: error_estimate, 01531 preconditioner_shift 01532 REAL(KIND=dp), DIMENSION(:), POINTER :: c0_shift, diag, evals 01533 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 01534 TYPE(cp_fm_type), POINTER :: matrix_hc0, matrix_s1, 01535 matrix_sc0, matrix_tmp, 01536 matrix_tmp2 01537 01538 ! arbitrary upshift of the occupied evals 01539 ! fudge factor for taking the error estimate into account 01540 01541 CALL timeset(routineN,handle) 01542 01543 CALL cp_fm_get_info(matrix_c0,nrow_global=n,ncol_global=k,error=error) 01544 01545 IF (ASSOCIATED(preconditioner_env%fm)) CALL cp_fm_release(preconditioner_env%fm,error=error) 01546 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=n,ncol_global=n, & 01547 context=preconditioner_env%ctxt, & 01548 para_env=preconditioner_env%para_env,error=error) 01549 CALL cp_fm_create(preconditioner_env%fm,fm_struct_tmp,name="preconditioner_env%fm",error=error) 01550 CALL cp_fm_create(matrix_tmp,fm_struct_tmp, name="preconditioner matrix_tmp",error=error) 01551 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01552 01553 ! get the error estimate 01554 CALL cp_fm_create(matrix_sc0,matrix_c0%matrix_struct,name="sc0",error=error) 01555 CALL cp_fm_to_fm(matrix_c0,matrix_sc0,error=error) 01556 CALL cp_fm_create(matrix_hc0,matrix_c0%matrix_struct,name="hc0",error=error) 01557 CALL cp_dbcsr_sm_fm_multiply(matrix_h,matrix_c0,matrix_hc0,k,error=error) 01558 01559 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=k,ncol_global=k, & 01560 context=preconditioner_env%ctxt, & 01561 para_env=preconditioner_env%para_env,error=error) 01562 CALL cp_fm_create(matrix_s1,fm_struct_tmp,name="matrix_s1",error=error) 01563 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01564 ! since we only use diagonal elements this is a bit of a waste 01565 CALL cp_fm_gemm('T','N',k,k,n,1.0_dp,matrix_hc0,matrix_hc0,0.0_dp,matrix_s1,error=error) 01566 ALLOCATE(diag(k)) 01567 CALL cp_fm_get_diag(matrix_s1,diag,error=error) 01568 error_estimate=MAXVAL(SQRT(ABS(diag-c0_evals**2))) 01569 DEALLOCATE(diag) 01570 CALL cp_fm_release(matrix_s1,error=error) 01571 CALL cp_fm_release(matrix_hc0,error=error) 01572 01573 ! shift up the occupied subspace eigenvalues 01574 ALLOCATE(c0_shift(k)) 01575 c0_shift=SQRT(-(c0_evals-c0_evals(k))+eval_shift) 01576 CALL cp_fm_column_scale(matrix_sc0,c0_shift) 01577 CALL cp_fm_gemm('N','T',n,n,k,1.0_dp,matrix_sc0,matrix_sc0,0.0_dp,preconditioner_env%fm,error=error) 01578 CALL cp_fm_release(matrix_sc0,error=error) 01579 DEALLOCATE(c0_shift) 01580 01581 ! get H added to the shift 01582 CALL copy_dbcsr_to_fm(matrix_h,matrix_tmp,error=error) 01583 CALL cp_fm_scale_and_add(1.0_dp,preconditioner_env%fm,1.0_dp,matrix_tmp,error=error) 01584 01585 ! preconditioner shift, we target the middle of the occupied spectrum, and taking into account the error_estimate 01586 ! write(6,*) "Error estimate = ",error_estimate 01587 preconditioner_shift=-(MINVAL(c0_evals)+ MAXVAL(c0_evals))/2.0_dp + & 01588 error_estimate*fudge_factor 01589 CALL cp_fm_set_all(matrix_tmp,alpha=0._dp,beta=1._dp,error=error) 01590 CALL cp_fm_scale_and_add(1.0_dp,preconditioner_env%fm,preconditioner_shift,matrix_tmp,error=error) 01591 ! check evals 01592 IF (.FALSE.) THEN 01593 CALL cp_fm_to_fm(preconditioner_env%fm,matrix_tmp,error=error) 01594 CALL cp_fm_create(matrix_tmp2,matrix_tmp%matrix_struct,name="matrix_tmp2",error=error) 01595 ALLOCATE(evals(n)) 01596 CALL choose_eigv_solver(matrix_tmp,matrix_tmp2,evals,error=error) 01597 CALL cp_fm_release(matrix_tmp2,error=error) 01598 WRITE(6,*) "evals ",evals 01599 DEALLOCATE(evals) 01600 ENDIF 01601 01602 CALL cp_fm_release(matrix_tmp,error=error) 01603 01604 CALL timestop(handle) 01605 01606 END SUBROUTINE make_full_single_inverse_ortho 01607 ! ***************************************************************************** 01608 SUBROUTINE make_full_all_ortho(preconditioner_env, matrix_c0, matrix_h, c0_evals, energy_gap, error) 01609 01610 TYPE(preconditioner_type) :: preconditioner_env 01611 TYPE(cp_fm_type), POINTER :: matrix_c0 01612 TYPE(cp_dbcsr_type), POINTER :: matrix_h 01613 REAL(KIND=dp), DIMENSION(:), POINTER :: c0_evals 01614 REAL(KIND=dp) :: energy_gap 01615 TYPE(cp_error_type), INTENT(inout) :: error 01616 01617 CHARACTER(len=*), PARAMETER :: routineN = 'make_full_all_ortho', 01618 routineP = moduleN//':'//routineN 01619 REAL(KIND=dp), PARAMETER :: fudge_factor = 0.25_dp, 01620 lambda_base = 10.0_dp 01621 01622 INTEGER :: handle, k, n 01623 REAL(KIND=dp) :: error_estimate, lambda 01624 REAL(KIND=dp), DIMENSION(:), POINTER :: diag, norms, shifted_evals 01625 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 01626 TYPE(cp_fm_type), POINTER :: matrix_hc0, matrix_left, 01627 matrix_pre, matrix_s1, 01628 matrix_s2, matrix_sc0, 01629 matrix_tmp 01630 01631 CALL timeset(routineN,handle) 01632 01633 IF (ASSOCIATED(preconditioner_env%fm)) CALL cp_fm_release(preconditioner_env%fm,error) 01634 CALL cp_fm_get_info(matrix_c0,nrow_global=n,ncol_global=k,error=error) 01635 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=n,ncol_global=n, & 01636 context=preconditioner_env%ctxt, & 01637 para_env=preconditioner_env%para_env,error=error) 01638 CALL cp_fm_create(preconditioner_env%fm,fm_struct_tmp,name="preconditioner_env%fm",error=error) 01639 matrix_pre=>preconditioner_env%fm 01640 CALL cp_fm_create(matrix_tmp,fm_struct_tmp,name="matrix_tmp",error=error) 01641 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01642 ALLOCATE(preconditioner_env%full_evals(n)) 01643 ALLOCATE(preconditioner_env%occ_evals(k)) 01644 01645 ! 1) Construct a new H matrix, which has the current C0 as eigenvectors, 01646 ! possibly shifted by an amount lambda, 01647 ! and the same spectrum as the original H matrix in the space orthogonal to the C0 01648 ! with P=C0 C0 ^ T 01649 ! (1 - PS)^T H (1-PS) + (PS)^T (H - lambda S ) (PS) 01650 ! we exploit that the C0 are already the ritz states of H 01651 CALL cp_fm_create(matrix_sc0,matrix_c0%matrix_struct,name="sc0",error=error) 01652 CALL cp_fm_to_fm(matrix_c0,matrix_sc0,error=error) 01653 CALL cp_fm_create(matrix_hc0,matrix_c0%matrix_struct,name="hc0",error=error) 01654 CALL cp_dbcsr_sm_fm_multiply(matrix_h,matrix_c0,matrix_hc0,k,error=error) 01655 01656 ! An aside, try to estimate the error on the ritz values, we'll need it later on 01657 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=k,ncol_global=k, & 01658 context=preconditioner_env%ctxt, & 01659 para_env=preconditioner_env%para_env,error=error) 01660 CALL cp_fm_create(matrix_s1,fm_struct_tmp,name="matrix_s1",error=error) 01661 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01662 ! since we only use diagonal elements this is a bit of a waste 01663 CALL cp_fm_gemm('T','N',k,k,n,1.0_dp,matrix_hc0,matrix_hc0,0.0_dp,matrix_s1,error=error) 01664 ALLOCATE(diag(k)) 01665 CALL cp_fm_get_diag(matrix_s1,diag,error=error) 01666 error_estimate=MAXVAL(SQRT(ABS(diag-c0_evals**2))) 01667 DEALLOCATE(diag) 01668 CALL cp_fm_release(matrix_s1,error=error) 01669 ! we'll only use the energy gap, if our estimate of the error on the eigenvalues 01670 ! is small enough. A large error combined with a small energy gap would otherwise lead to 01671 ! an aggressive but bad preconditioner. Only when the error is small (MD), we can precondition 01672 ! aggressively 01673 preconditioner_env%energy_gap= MAX(energy_gap,error_estimate*fudge_factor) 01674 01675 CALL copy_dbcsr_to_fm(matrix_h,matrix_tmp,error=error) 01676 CALL cp_fm_upper_to_full(matrix_tmp,matrix_pre,error=error) 01677 ! tmp = H ( 1 - PS ) 01678 CALL cp_fm_gemm('N','T',n,n,k,-1.0_dp,matrix_hc0,matrix_sc0,1.0_dp,matrix_tmp,error=error) 01679 01680 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=k,ncol_global=n, & 01681 context=preconditioner_env%ctxt, & 01682 para_env=preconditioner_env%para_env,error=error) 01683 CALL cp_fm_create(matrix_left,fm_struct_tmp,name="matrix_left",error=error) 01684 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01685 CALL cp_fm_gemm('T','N',k,n,n,1.0_dp,matrix_c0,matrix_tmp,0.0_dp,matrix_left,error=error) 01686 ! tmp = (1 - PS)^T H (1-PS) 01687 CALL cp_fm_gemm('N','N',n,n,k,-1.0_dp,matrix_sc0,matrix_left,1.0_dp,matrix_tmp,error=error) 01688 CALL cp_fm_release(matrix_left,error=error) 01689 01690 ALLOCATE(shifted_evals(k)) 01691 lambda = lambda_base + error_estimate 01692 shifted_evals=c0_evals - lambda 01693 CALL cp_fm_to_fm(matrix_sc0,matrix_hc0,error=error) 01694 CALL cp_fm_column_scale(matrix_hc0,shifted_evals) 01695 CALL cp_fm_gemm('N','T',n,n,k,1.0_dp,matrix_hc0,matrix_sc0,1.0_dp,matrix_tmp,error=error) 01696 01697 ! 2) diagonalize this operator 01698 CALL choose_eigv_solver(matrix_tmp,matrix_pre,preconditioner_env%full_evals,error=error) 01699 01700 01701 ! test that the subspace remained conserved 01702 IF (.FALSE.) THEN 01703 CALL cp_fm_to_fm(matrix_pre,matrix_tmp,error=error) 01704 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=k,ncol_global=k, & 01705 context=preconditioner_env%ctxt, & 01706 para_env=preconditioner_env%para_env,error=error) 01707 CALL cp_fm_create(matrix_s1,fm_struct_tmp,name="matrix_s1",error=error) 01708 CALL cp_fm_create(matrix_s2,fm_struct_tmp,name="matrix_s2",error=error) 01709 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01710 ALLOCATE(norms(k)) 01711 CALL cp_fm_gemm('T','N',k,k,n,1.0_dp,matrix_sc0,matrix_tmp,0.0_dp,matrix_s1,error=error) 01712 CALL choose_eigv_solver(matrix_s1,matrix_s2,norms,error=error) 01713 01714 WRITE(6,*) "matrix norm deviation (should be close to zero): ", MAXVAL(ABS(ABS(norms)-1.0_dp)) 01715 DEALLOCATE(norms) 01716 CALL cp_fm_release(matrix_s1,error=error) 01717 CALL cp_fm_release(matrix_s2,error=error) 01718 ENDIF 01719 01720 ! 3) replace the lowest k evals and evecs with what they should be 01721 preconditioner_env%occ_evals=c0_evals 01722 ! notice, this choice causes the preconditioner to be constant when applied to sc0 (see apply_full_all) 01723 preconditioner_env%full_evals(1:k)=c0_evals 01724 CALL cp_fm_to_fm(matrix_c0,matrix_pre,k,1,1) 01725 01726 CALL cp_fm_release(matrix_sc0,error=error) 01727 CALL cp_fm_release(matrix_hc0,error=error) 01728 CALL cp_fm_release(matrix_tmp,error=error) 01729 DEALLOCATE(shifted_evals) 01730 01731 CALL timestop(handle) 01732 01733 01734 END SUBROUTINE make_full_all_ortho 01735 ! ***************************************************************************** 01736 01737 SUBROUTINE make_diag_inner_precond(preconditioner_env, sm, sm_inner, error) 01738 ! 01739 TYPE(preconditioner_type) :: preconditioner_env 01740 TYPE(cp_dbcsr_type), POINTER :: sm, sm_inner 01741 TYPE(cp_error_type), INTENT(inout) :: error 01742 01743 CHARACTER(len=*), PARAMETER :: routineN = 'make_diag_inner_precond', 01744 routineP = moduleN//':'//routineN 01745 01746 INTEGER :: blk, handle, i, iblock_col, 01747 iblock_row, info, istat, 01748 liwork, lwork, n, nblocks 01749 INTEGER, DIMENSION(:), POINTER :: iwork 01750 LOGICAL :: failure 01751 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals, work 01752 REAL(KIND=dp), ALLOCATABLE, 01753 DIMENSION(:, :) :: block_evec 01754 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block_buf1, block_h, block_pre 01755 TYPE(cp_dbcsr_iterator) :: iter 01756 01757 CALL timeset(routineN,handle) 01758 failure = .FALSE. 01759 01760 CPPrecondition(ASSOCIATED(sm),cp_failure_level,routineP,error,failure) 01761 01762 IF(ASSOCIATED(sm_inner)) THEN 01763 CALL cp_dbcsr_deallocate_matrix(sm_inner,error=error) 01764 NULLIFY(sm_inner) 01765 ENDIF 01766 01767 ALLOCATE(sm_inner) 01768 CALL cp_dbcsr_init(sm_inner,error=error) 01769 CALL cp_dbcsr_create(sm_inner, ' PRECONDITIONER ', & 01770 cp_dbcsr_distribution (sm), cp_dbcsr_get_matrix_type (sm),& 01771 cp_dbcsr_row_block_sizes(sm),& 01772 cp_dbcsr_col_block_sizes(sm), 0, 0, error=error) 01773 01774 CALL cp_dbcsr_get_info(sm,nfullrows_total=nblocks) 01775 01776 CALL cp_dbcsr_iterator_start(iter, sm) 01777 DO WHILE (cp_dbcsr_iterator_blocks_left(iter)) 01778 CALL cp_dbcsr_iterator_next_block(iter, iblock_row, iblock_col, block_h, blk) 01779 01780 IF (iblock_col .EQ. iblock_row) THEN 01781 n = SIZE(block_h,1) 01782 ALLOCATE(block_pre(n,n),block_evec(n,n),block_buf1(n,n),evals(n)) 01783 lwork=1+6*n+2*n**2+50 01784 liwork=5*n+3 01785 ALLOCATE (work(lwork),STAT=istat) 01786 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01787 "work",dp_size*lwork) 01788 ALLOCATE (iwork(liwork),STAT=istat) 01789 IF (istat /= 0) CALL stop_memory(routineN,moduleN,__LINE__,& 01790 "iwork",int_size*liwork) 01791 ! 01792 ! A simple block diagonal preconditoner 01793 block_evec(:,:) = block_h(:,:) 01794 CALL DSYEVD('V','U', n, block_evec(1,1), n, evals(1), work(1), lwork, & 01795 iwork(1), liwork, info) 01796 IF (info.ne.0) CALL stop_program(routineN,moduleN,__LINE__,"DSYEVD H") 01797 block_pre(:,:)=0.0_dp 01798 DO i=1,n 01799 IF (evals(i).lt.1e-10_dp) CALL stop_program(routineN,moduleN,__LINE__,& 01800 "evals(i).lt.1e-10_dp") 01801 block_pre(i,i)=1.0_dp/evals(i) 01802 ENDDO 01803 CALL DGEMM('N','N',n,n,n,1.0_dp,block_evec(1,1),n,block_pre(1,1),n, & 01804 0.0_dp,block_buf1(1,1),n) 01805 CALL DGEMM('N','T',n,n,n,1.0_dp,block_buf1(1,1),n,block_evec(1,1),n, & 01806 0.0_dp,block_pre(1,1),n) 01807 ! 01808 ! add the block 01809 CALL cp_dbcsr_put_block(matrix=sm_inner,& 01810 row=iblock_row,& 01811 col=iblock_col,& 01812 BLOCK=block_pre) 01813 01814 DEALLOCATE(block_pre,block_evec,block_buf1,evals,work,iwork) 01815 01816 ENDIF 01817 ENDDO 01818 CALL cp_dbcsr_iterator_stop(iter) 01819 CALL cp_dbcsr_finalize(sm_inner,error=error) 01820 01821 CALL timestop(handle) 01822 01823 END SUBROUTINE make_diag_inner_precond 01824 01825 SUBROUTINE apply_solve_lin_system_fm(preconditioner_env, matrix_in, matrix_out, error) 01826 ! 01827 TYPE(preconditioner_type) :: preconditioner_env 01828 TYPE(cp_fm_type), POINTER :: matrix_in, matrix_out 01829 TYPE(cp_error_type), INTENT(inout) :: error 01830 01831 CHARACTER(len=*), PARAMETER :: routineN = 'apply_solve_lin_system_fm', 01832 routineP = moduleN//':'//routineN 01833 INTEGER, PARAMETER :: max_iter = 200, 01834 restart_every = 50 01835 REAL(dp), PARAMETER :: eps = 1.0e-2_dp 01836 01837 INTEGER :: iter, m, n, output_unit 01838 LOGICAL :: failure 01839 REAL(dp) :: alpha, beta, residual_ot, 01840 tr_pAp, tr_rr, tr_sr_new, 01841 tr_sr_old 01842 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 01843 TYPE(cp_fm_type), POINTER :: Ap, p, r, s 01844 TYPE(cp_logger_type), POINTER :: logger 01845 01846 failure = .FALSE. 01847 CPPrecondition(ASSOCIATED(preconditioner_env%sparse_matrix),cp_failure_level,routineP,error,failure) 01848 ! 01849 CALL cp_fm_get_info(matrix_in,nrow_global=n,ncol_global=m,error=error) 01850 CALL cp_fm_struct_create(fm_struct_tmp,nrow_global=n,ncol_global=m,& 01851 context=preconditioner_env%ctxt, & 01852 para_env=preconditioner_env%para_env,error=error) 01853 CALL cp_fm_create(r ,fm_struct_tmp,name="solve_r" ,error=error) 01854 CALL cp_fm_create(p ,fm_struct_tmp,name="solve_p" ,error=error) 01855 CALL cp_fm_create(s ,fm_struct_tmp,name="solve_s" ,error=error) 01856 CALL cp_fm_create(Ap,fm_struct_tmp,name="solve_Ap",error=error) 01857 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 01858 CALL cp_fm_set_all(p ,0.0_dp,error=error) 01859 ! 01860 logger => cp_error_get_logger(error) 01861 output_unit= cp_logger_get_default_io_unit(logger) 01862 ! 01863 ! residual 01864 CALL cp_fm_trace(matrix_in,matrix_in,residual_ot,error=error) 01865 residual_ot = SQRT(residual_ot/(REAL(n,dp)*REAL(m,dp))) 01866 ! 01867 ! r = b-A*x 01868 CALL cp_dbcsr_sm_fm_multiply(preconditioner_env%sparse_matrix,matrix_out,r,m,error=error) 01869 CALL cp_fm_scale_and_add(-1.0_dp,r,1.0_dp,matrix_in,error=error) 01870 CALL cp_fm_trace(r,r,tr_rr,error=error) 01871 IF(output_unit>0) THEN 01872 WRITE(output_unit,'(A,I3,2(A,E12.5))') 'apply_solve_lin_system:',0,': ||r_ls||=',& 01873 & SQRT(tr_rr/(REAL(n,dp)*REAL(m,dp))),' ||r_ot||=',residual_ot 01874 ENDIF 01875 ! 01876 ! let's go! 01877 DO iter = 1,max_iter 01878 ! 01879 ! s = M * r 01880 IF(ASSOCIATED(preconditioner_env%sparse_matrix_inner)) THEN 01881 CALL cp_dbcsr_sm_fm_multiply(preconditioner_env%sparse_matrix_inner,r,s,m,error=error) 01882 ELSE 01883 CALL cp_fm_to_fm(r,s,error=error) 01884 ENDIF 01885 CALL cp_fm_trace(s,r,tr_sr_new,error=error) 01886 ! 01887 ! beta (this might be better when O(N) beta=(s_k,r_k-r_{k-1})/(s_{k-1},r_{k-1})) 01888 IF(iter.EQ.1) THEN 01889 beta = 0.0_dp 01890 ELSE 01891 beta = tr_sr_new/tr_sr_old 01892 ENDIF 01893 ! 01894 ! p = r + beta * p 01895 CALL cp_fm_scale_and_add(beta,p,1.0_dp,s,error=error) 01896 ! 01897 ! Ap = A * p 01898 CALL cp_dbcsr_sm_fm_multiply(preconditioner_env%sparse_matrix,p,Ap,m,error=error) 01899 CALL cp_fm_trace(Ap,p,tr_pAp,error=error) 01900 ! 01901 ! alpha = Tr(s_k'*r_k)/Tr(p_k'*Ap_k) 01902 alpha = tr_sr_new/tr_pAp 01903 ! 01904 ! x = x + alpha * p 01905 CALL cp_fm_scale_and_add(1.0_dp,matrix_out,alpha,p,error=error) 01906 ! 01907 ! r = r - alpha * Ap or r = b - A * x 01908 IF(MOD(iter,restart_every).EQ.0) THEN 01909 CALL cp_dbcsr_sm_fm_multiply(preconditioner_env%sparse_matrix,matrix_out,r,m,error=error) 01910 CALL cp_fm_scale_and_add(-1.0_dp,r,1.0_dp,matrix_in,error=error) 01911 ELSE 01912 CALL cp_fm_scale_and_add(1.0_dp,r,-alpha,Ap,error=error) 01913 ENDIF 01914 tr_sr_old = tr_sr_new 01915 ! 01916 ! printing 01917 IF(MOD(iter,10).EQ.0.AND.output_unit>0) THEN 01918 WRITE(output_unit,'(A,I3,2(A,E12.5))') 'apply_solve_lin_system:',iter,': ||r_ls||=',& 01919 & SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp))),' ||r_ls||/||r_ot||=', 01920 SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp)))/residual_ot 01921 ENDIF 01922 ! 01923 ! exit 01924 IF(SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp)))/residual_ot.LT.eps) EXIT 01925 ! 01926 ENDDO 01927 ! 01928 IF(output_unit>0) THEN 01929 WRITE(output_unit,'(A,I3,2(A,E12.5))') 'apply_solve_lin_system:',iter,': ||r_ls||=',& 01930 & SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp))),' ||r_ls||/||r_ot||=', 01931 SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp)))/residual_ot 01932 ENDIF 01933 IF(iter.GT.max_iter) THEN 01934 CALL stop_program(routineN,moduleN,__LINE__,"Increase max_iter") 01935 ENDIF 01936 ! 01937 CALL cp_fm_release(r ,error=error) 01938 CALL cp_fm_release(p ,error=error) 01939 CALL cp_fm_release(s ,error=error) 01940 CALL cp_fm_release(Ap,error=error) 01941 END SUBROUTINE apply_solve_lin_system_fm 01942 01943 SUBROUTINE apply_solve_lin_system_dbcsr(preconditioner_env, matrix_in, matrix_out, error) 01944 ! 01945 TYPE(preconditioner_type) :: preconditioner_env 01946 TYPE(cp_dbcsr_type) :: matrix_in, matrix_out 01947 TYPE(cp_error_type), INTENT(inout) :: error 01948 01949 CHARACTER(len=*), PARAMETER :: routineN = 'apply_solve_lin_system_dbcsr', 01950 routineP = moduleN//':'//routineN 01951 INTEGER, PARAMETER :: max_iter = 200, 01952 restart_every = 50 01953 REAL(dp), PARAMETER :: eps = 1.0e-2_dp 01954 01955 INTEGER :: iter, m, n, output_unit 01956 LOGICAL :: failure 01957 REAL(dp) :: alpha, beta, residual_ot, 01958 tr_pAp, tr_rr, tr_sr_new, 01959 tr_sr_old 01960 TYPE(cp_dbcsr_type) :: Ap, p, r, s 01961 TYPE(cp_logger_type), POINTER :: logger 01962 01963 failure = .FALSE. 01964 CPPrecondition(ASSOCIATED(preconditioner_env%sparse_matrix),cp_failure_level,routineP,error,failure) 01965 ! 01966 CALL cp_dbcsr_get_info(matrix_in,nfullrows_total=n,nfullcols_total=m) 01967 CALL cp_dbcsr_init(r,error=error) 01968 CALL cp_dbcsr_init(p,error=error) 01969 CALL cp_dbcsr_init(s,error=error) 01970 CALL cp_dbcsr_init(Ap,error=error) 01971 CALL cp_dbcsr_copy(r,matrix_in,error=error) 01972 CALL cp_dbcsr_copy(p,matrix_in,error=error) 01973 CALL cp_dbcsr_copy(s,matrix_in,error=error) 01974 CALL cp_dbcsr_copy(Ap,matrix_in,error=error) 01975 ! 01976 logger => cp_error_get_logger(error) 01977 output_unit= cp_logger_get_default_io_unit(logger) 01978 ! 01979 ! residual 01980 CALL cp_dbcsr_trace(matrix_in,matrix_in,residual_ot,error=error) 01981 residual_ot = SQRT(residual_ot/(REAL(n,dp)*REAL(m,dp))) 01982 ! 01983 ! r = b-A*x 01984 CALL cp_dbcsr_multiply('N','N',1.0_dp,preconditioner_env%sparse_matrix,matrix_out,& 01985 0.0_dp,r,error=error) 01986 CALL cp_dbcsr_add(r,matrix_in,-1.0_dp,1.0_dp,error=error) 01987 CALL cp_dbcsr_trace(r,r,tr_rr,error=error) 01988 IF(output_unit>0) THEN 01989 WRITE(output_unit,'(A,I3,2(A,E12.5))') 'apply_solve_lin_system:',0,': ||r_ls||=',& 01990 & SQRT(tr_rr/(REAL(n,dp)*REAL(m,dp))),' ||r_ot||=',residual_ot 01991 ENDIF 01992 ! 01993 ! let's go! 01994 DO iter = 1,max_iter 01995 ! 01996 ! s = M * r 01997 IF(ASSOCIATED(preconditioner_env%sparse_matrix_inner)) THEN 01998 CALL cp_dbcsr_multiply('N','N',1.0_dp,preconditioner_env%sparse_matrix_inner,r,& 01999 0.0_dp,s,error=error) 02000 ELSE 02001 CALL cp_dbcsr_copy(s,r,error=error) 02002 ENDIF 02003 CALL cp_dbcsr_trace(s,r,tr_sr_new,error=error) 02004 ! 02005 ! beta (this might be better when O(N) beta=(s_k,r_k-r_{k-1})/(s_{k-1},r_{k-1})) 02006 IF(iter.EQ.1) THEN 02007 beta = 0.0_dp 02008 ELSE 02009 beta = tr_sr_new/tr_sr_old 02010 ENDIF 02011 ! 02012 ! p = r + beta * p 02013 CALL cp_dbcsr_add(p,s,beta,1.0_dp,error=error) 02014 ! 02015 ! Ap = A * p 02016 CALL cp_dbcsr_multiply('N','N',1.0_dp,preconditioner_env%sparse_matrix,p,& 02017 0.0_dp,Ap,error=error) 02018 CALL cp_dbcsr_trace(Ap,p,tr_pAp,error=error) 02019 ! 02020 ! alpha = Tr(s_k'*r_k)/Tr(p_k'*Ap_k) 02021 alpha = tr_sr_new/tr_pAp 02022 ! 02023 ! x = x + alpha * p 02024 CALL cp_dbcsr_add(matrix_out,p,1.0_dp,alpha,error=error) 02025 ! 02026 ! r = r - alpha * Ap or r = b - A * x 02027 IF(MOD(iter,restart_every).EQ.0) THEN 02028 CALL cp_dbcsr_multiply('N','N',1.0_dp,preconditioner_env%sparse_matrix,matrix_out,& 02029 0.0_dp,r,error=error) 02030 CALL cp_dbcsr_add(r,matrix_in,-1.0_dp,1.0_dp,error=error) 02031 ELSE 02032 CALL cp_dbcsr_add(r,Ap,1.0_dp,-alpha,error=error) 02033 ENDIF 02034 tr_sr_old = tr_sr_new 02035 ! 02036 ! printing 02037 IF(MOD(iter,10).EQ.0.AND.output_unit>0) THEN 02038 WRITE(output_unit,'(A,I3,2(A,E12.5))') 'apply_solve_lin_system:',iter,': ||r_ls||=',& 02039 & SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp))),' ||r_ls||/||r_ot||=', 02040 SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp)))/residual_ot 02041 ENDIF 02042 ! 02043 ! exit 02044 IF(SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp)))/residual_ot.LT.eps) EXIT 02045 ! 02046 ENDDO 02047 ! 02048 IF(output_unit>0) THEN 02049 WRITE(output_unit,'(A,I3,2(A,E12.5))') 'apply_solve_lin_system:',iter,': ||r_ls||=',& 02050 & SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp))),' ||r_ls||/||r_ot||=', 02051 SQRT(tr_sr_new/(REAL(n,dp)*REAL(m,dp)))/residual_ot 02052 ENDIF 02053 IF(iter.GT.max_iter) THEN 02054 CALL stop_program(routineN,moduleN,__LINE__,"increase max_iter") 02055 ENDIF 02056 ! 02057 CALL cp_dbcsr_release(r, error=error) 02058 CALL cp_dbcsr_release(p, error=error) 02059 CALL cp_dbcsr_release(s, error=error) 02060 CALL cp_dbcsr_release(Ap, error=error) 02061 END SUBROUTINE apply_solve_lin_system_dbcsr 02062 02063 ! ***************************************************************************** 02074 SUBROUTINE make_sparse_inverse_bif(preconditioner_env, fm, matrix_s, error) 02075 ! 02076 TYPE(preconditioner_type), INTENT(INOUT) :: preconditioner_env 02077 TYPE(cp_fm_type), POINTER :: fm 02078 TYPE(cp_dbcsr_type), POINTER :: matrix_s 02079 TYPE(cp_error_type), INTENT(INOUT) :: error 02080 02081 CHARACTER(len=*), PARAMETER :: routineN = 'make_sparse_inverse_bif', 02082 routineP = moduleN//':'//routineN 02083 REAL(KIND=dp), PARAMETER :: drop_v = 1.0E-5_dp, 02084 thresh_inner = 1.0E-2_dp, 02085 thresh_post = 1.0E-5_dp 02086 02087 INTEGER :: d_offset, handle, istat, k, k_s, max_blk_size, 02088 max_row_blk_size, max_vec_size, mp_group, mynode, mypcol, myprow, n, 02089 nblkcols_local, nblkrows_local, nblkrows_total, npcols, nprows, 02090 numnodes, proc_holds_diag_blk 02091 INTEGER, DIMENSION(:), POINTER :: row_blk_size 02092 LOGICAL :: dbg_print_chksum, 02093 dbg_print_matrix, failure, 02094 found, tr 02095 REAL(dp) :: chksum, occ_in, occ_out, s, 02096 t(100), trace 02097 REAL(KIND=dp), DIMENSION(:), POINTER :: blk, blk_p, pkd_d, pkd_inv_d, 02098 pkd_u, pkd_u_fac, pkd_v, 02099 pkd_v_fac, pkd_v_tmp 02100 TYPE(cp_dbcsr_type) :: mat_li, mat_li_scaled, 02101 mat_lt, mat_lt_scaled, mat_t, 02102 mat_tmp, mat_v 02103 TYPE(cp_dbcsr_type), POINTER :: prcnd 02104 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 02105 TYPE(cp_fm_type), POINTER :: fm_work1, fm_work2 02106 TYPE(dbcsr_error_type) :: dbcsr_error 02107 TYPE(dbcsr_mp_obj) :: mp_obj 02108 02109 t=0.0_dp 02110 dbg_print_matrix = .FALSE. 02111 dbg_print_chksum = .TRUE. 02112 ! 02113 02114 failure = .FALSE. 02115 02116 CALL timeset(routineN,handle) 02117 02118 CALL cpu_time(t(1)) 02119 02120 CPPrecondition(ASSOCIATED(fm),cp_failure_level,routineP,error,failure) 02121 02122 ! 02123 ! some infos 02124 CALL cp_fm_get_info(fm,nrow_global=n,error=error) 02125 ! 02126 ! allocate 02127 CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=n,ncol_global=n,& 02128 context=preconditioner_env%ctxt, & 02129 para_env=preconditioner_env%para_env,error=error) 02130 CALL cp_fm_create(fm_work1,fm_struct_tmp, name="preconditioner_work1",error=error) 02131 CALL cp_fm_create(fm_work2,fm_struct_tmp, name="preconditioner_work2",error=error) 02132 CALL cp_fm_struct_release(fm_struct_tmp,error=error) 02133 ! 02134 ! here fm should contain the matrix to inverse 02135 CALL cp_fm_to_fm(fm,fm_work1,error=error) 02136 CALL cp_fm_upper_to_full(fm_work1,fm_work2,error=error) 02137 ! 02138 ! 02139 ! copy matrices 02140 CALL cp_dbcsr_init (mat_t, error) 02141 CALL cp_dbcsr_create(mat_t, 'mat t', & 02142 cp_dbcsr_distribution(matrix_s), dbcsr_type_no_symmetry, cp_dbcsr_row_block_sizes(matrix_s),& 02143 cp_dbcsr_col_block_sizes(matrix_s), 0, 0,& 02144 cp_dbcsr_get_data_type(matrix_s), error=error) 02145 CALL copy_fm_to_dbcsr(fm_work1,mat_t, error=error) 02146 IF(.TRUE.)CALL cp_dbcsr_verify_matrix(mat_t, error=error) 02147 IF(dbg_print_matrix)CALL cp_dbcsr_print(mat_t,error=error) 02148 02149 occ_in = cp_dbcsr_get_occupation(mat_t) 02150 CALL cp_dbcsr_filter(mat_t,1e-5_dp,error=error) 02151 IF(.FALSE.)CALL cp_dbcsr_verify_matrix(mat_t, error=error) 02152 occ_out = cp_dbcsr_get_occupation(mat_t) 02153 02154 WRITE(*,*) routineN//' occ_in',occ_in,' occ_out',occ_out 02155 02156 CALL cp_dbcsr_init (mat_li,error=error) 02157 CALL cp_dbcsr_create(mat_li, 'matrix li', & 02158 cp_dbcsr_distribution(mat_t), dbcsr_type_no_symmetry, cp_dbcsr_row_block_sizes(mat_t),& 02159 cp_dbcsr_col_block_sizes(mat_t), 0, 0,& 02160 cp_dbcsr_get_data_type(mat_t), error=error) 02161 CALL cp_dbcsr_init (mat_v,error=error) 02162 CALL cp_dbcsr_create(mat_v, 'matrix v', & 02163 cp_dbcsr_distribution(mat_t), dbcsr_type_no_symmetry, cp_dbcsr_row_block_sizes(mat_t),& 02164 cp_dbcsr_col_block_sizes(mat_t), 0, 0,& 02165 cp_dbcsr_get_data_type(mat_t), error=error) 02166 02167 CALL cp_dbcsr_work_create(mat_li,work_mutable=.TRUE.,& 02168 nblks_guess=cp_dbcsr_get_num_blocks(mat_t),sizedata_guess=cp_dbcsr_get_data_size(mat_t),& 02169 error=error) 02170 CALL cp_dbcsr_work_create(mat_v,work_mutable=.TRUE.,& 02171 nblks_guess=cp_dbcsr_get_num_blocks(mat_t),sizedata_guess=cp_dbcsr_get_data_size(mat_t),& 02172 error=error) 02173 ! 02174 ! 02175 prcnd => preconditioner_env%sparse_matrix 02176 IF (ASSOCIATED(prcnd)) THEN 02177 CALL cp_dbcsr_deallocate_matrix(prcnd,error=error) 02178 NULLIFY(prcnd) 02179 ENDIF 02180 ALLOCATE(prcnd) 02181 CALL cp_dbcsr_init(prcnd,error=error) 02182 CALL cp_dbcsr_create(prcnd, ' PRECONDITIONER ', & 02183 cp_dbcsr_distribution(mat_t), dbcsr_type_no_symmetry, cp_dbcsr_row_block_sizes(mat_t),& 02184 cp_dbcsr_col_block_sizes(mat_t), cp_dbcsr_get_num_blocks(mat_t), cp_dbcsr_get_data_size(mat_t),& 02185 cp_dbcsr_get_data_type(mat_t), error=error) 02186 02187 ! some variables 02188 s = 1.0_dp ! s > 0 02189 02190 nblkrows_total = cp_dbcsr_nblkrows_total(mat_t) 02191 nblkrows_local = cp_dbcsr_nblkrows_local(mat_t) 02192 nblkcols_local = cp_dbcsr_nblkcols_local(mat_t) 02193 row_blk_size => array_data(cp_dbcsr_row_block_sizes(mat_t)) 02194 mp_obj = dbcsr_distribution_mp (cp_dbcsr_distribution (mat_t)) 02195 mp_group = dbcsr_mp_group (mp_obj) 02196 02197 max_row_blk_size = MAXVAL(row_blk_size) 02198 max_blk_size = max_row_blk_size**2 02199 max_vec_size = 100 * MAX(nblkrows_local,nblkcols_local)*max_blk_size+nblkrows_total+1 ! need to fix that 02200 CALL mp_max(max_vec_size,mp_group) 02201 02202 WRITE(*,*) 'nblkrows_local',nblkrows_local 02203 WRITE(*,*) 'max_blk_size',max_blk_size 02204 WRITE(*,*) 'nblkrows_local',nblkrows_local 02205 WRITE(*,*) 'nblkcols_local',nblkcols_local 02206 WRITE(*,*) 'nblkrows_total',nblkrows_total 02207 WRITE(*,*) 'max_vec_size',max_vec_size 02208 02209 ALLOCATE(pkd_u(max_vec_size), pkd_v(max_vec_size), pkd_v_fac(max_vec_size),& 02210 & pkd_u_fac(max_vec_size), pkd_d(max_vec_size), pkd_inv_d(max_vec_size), & 02211 & blk(max_blk_size), pkd_v_tmp(max_vec_size), STAT=istat) 02212 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02213 02214 pkd_u(:) = 0.0_dp 02215 pkd_v(:) = 0.0_dp 02216 pkd_d(:) = 0.0_dp 02217 pkd_v_fac(:) = 0.0_dp 02218 pkd_u_fac(:) = 0.0_dp 02219 pkd_inv_d(:) = 0.0_dp 02220 blk(:) = 0.0_dp 02221 02222 numnodes = dbcsr_mp_numnodes (mp_obj) 02223 mynode = dbcsr_mp_mynode (mp_obj) 02224 myprow = dbcsr_mp_myprow (mp_obj) 02225 mypcol = dbcsr_mp_mypcol (mp_obj) 02226 npcols = dbcsr_mp_npcols (mp_obj) 02227 nprows = dbcsr_mp_nprows (mp_obj) 02228 !blacs2mpi => dbcsr_mp_pgrid (mp_obj) 02229 02230 WRITE(*,*) 'numnodes',numnodes 02231 WRITE(*,*) 'mynode',mynode 02232 WRITE(*,*) 'myprow',myprow 02233 WRITE(*,*) 'mypcol',mypcol 02234 WRITE(*,*) 'npcols',npcols 02235 WRITE(*,*) 'nprows',nprows 02236 02237 ! 02238 !CALL dbcsr_get_stored_coordinates (matrix, r, c, tr, p) 02239 !p .EQ. dbcsr_mp_mynode (mp_obj) 02240 ! 02241 02242 ! 02243 CALL cpu_time(t(2)) 02244 WRITE(*,*) 't1',t(2)-t(1) 02245 ! 02246 ! let's bif ! 02247 d_offset = 1 02248 DO k = 1,nblkrows_total 02249 !write(*,*) 'enter loop k=',k 02250 tr = .FALSE. 02251 k_s = k 02252 CALL cp_dbcsr_get_stored_coordinates (mat_t, k_s, k_s, tr, proc_holds_diag_blk) 02253 !write(*,*) 'proc_holds_diag_blk',proc_holds_diag_blk, proc_holds_diag_blk .EQ. dbcsr_mp_mynode (mp_obj) 02254 02255 CALL cpu_time(t(3)) 02256 !------------------------------------------------------------------------ 02257 ! set v and u 02258 ! 02259 !v(k,1:n)=A(k,1:n); 02260 CALL cp_dbcsr_copy_vec(mat_v, mat_t, 'row', k, error=error) 02261 !chksum = cp_dbcsr_checksum(mat_v,error) 02262 !write(*,*) 'checksum1 checksum(mat_v)',chksum 02263 CALL cpu_time(t(4)) 02264 t(50) = t(50) + t(4)-t(3) 02265 ! 02266 !v(k,k)=v(k,k)-s; 02267 IF(proc_holds_diag_blk .EQ. dbcsr_mp_mynode (mp_obj)) THEN 02268 CALL cp_dbcsr_get_block_p(mat_v, k, k, blk_p, found) 02269 IF (.NOT.found) CALL stop_program(routineN,moduleN,__LINE__,& 02270 "Matrix block not found") 02271 CALL block_add_on_diag(row_blk_size(k), blk_p, -s) 02272 ENDIF 02273 02274 !chksum = cp_dbcsr_checksum(mat_v,error) 02275 !write(*,*) 'checksum2 checksum(mat_v)',chksum 02276 ! 02277 !u(k,k)=1; 02278 !CALL block_set_d(row_blk_size(k), row_blk_size(k), blk, 1.0_dp, 0.0_dp) 02279 !CALL cp_dbcsr_put_block(mat_li, blk, k, k, error=error)!PARA local put 02280 !old CALL block_set(row_blk_size(k), row_blk_size(k), blk, 1.0_dp, 0.0_dp) 02281 !old CALL cp_dbcsr_put_block(mat_li, blk, k, k, error=error)!PARA local put 02282 ! 02283 !------------------------------------------------------------------------ 02284 ! build the factor inv(d)*u*A 02285 ! 02286 CALL cpu_time(t(5)) 02287 t(51) = t(51) + t(5)-t(4) 02288 ! compute u(1:k-1,1:n) * A(1:n,k); 02289 CALL cp_dbcsr_multiply_vec(mat_li, mat_t, 1, k-1, k, pkd_v_fac, error=error)!PARA pkd_v_fac is hold by 02290 ! !PARA the same procs as A(1:n,k) 02291 !write(*,*) 'pkd_v_fac:',pkd_v_fac(1:int(pkd_v_fac(nblkrows_total+1))-1) 02292 !chksum = sum(abs( pkd_v_fac ( nblkrows_total+2: int(pkd_v_fac (nblkrows_total+1))-1) )) 02293 !write(*,*) 'checksum3.0',chksum 02294 ! 02295 CALL cpu_time(t(6)) 02296 t(52) = t(52) + t(6)-t(5) 02297 ! 02298 ! bcast pkd_v_fac 02299 !CALL packed_vec_bcast(pkd_v_fac, source, 'columnwise', mp_obj, error=error) !PARA replicate columnwise 02300 ! 02301 ! compute inv(d(i,i)) * uAk(i,k) / s , i=1,k-1 02302 !chksum = sum(abs( pkd_inv_d )) 02303 !write(*,*) 'checksum3.1',chksum 02304 CALL packed_vec_scale(1.0_dp/s, pkd_inv_d, pkd_v_fac, k-1, k, row_blk_size, 'left', error=error)!PARA local 02305 CALL cpu_time(t(7)) 02306 t(53) = t(53) + t(7)-t(6) 02307 !------------------------------------------------------------------------ 02308 ! build the factor inv(d)*v 02309 ! 02310 ! pack v(:,k) 02311 !old CALL cp_dbcsr_pack_vec(mat_v, k, pkd_u_fac, 'column', error=error) 02312 ! 02313 ! bcast pkd_u_fac 02314 !CALL packed_vec_bcast(pkd_u_fac, source, 'columnwise', mp_obj, error=error) !PARA replicate columnwise 02315 ! 02316 ! compute inv(d(i,i)) * v(i,k) / s , i=1,k-1 02317 !old CALL packed_vec_scale(1.0_dp/s, pkd_inv_d, pkd_u_fac, k-1, k, row_blk_size, 'left', error=error)!PARA local 02318 !------------------------------------------------------------------------ 02319 ! update v and u 02320 ! 02321 ! pack u(k,:) and v(k,:) 02322 CALL packed_vec_ini(pkd_v, nblkrows_total, error) 02323 02324 !CALL cp_dbcsr_pack_vec(mat_v, k, pkd_v, 'row', error=error) !PARA local 02325 CALL cpu_time(t(8)) 02326 t(54) = t(54) + t(8)-t(7) 02327 !old CALL cp_dbcsr_pack_vec(mat_li, k, pkd_u, 'row', error=error) !PARA local 02328 ! 02329 ! bcast the u and v 02330 !call packed_vec_bcast(pkd_v, -1, 'all', .false., nblkrows_total, row_blk_size(k)*row_blk_size, mp_obj, error) 02331 !chksum = sum(abs(pkd_v(nblkrows_total+2:int(pkd_v(nblkrows_total+1))-1))) 02332 !write(*,*) 'checksum3.2 pkd_v',chksum 02333 !chksum = sum(abs(pkd_v_fac(nblkrows_total+2:int(pkd_v_fac(nblkrows_total+1))-1))) 02334 !write(*,*) 'checksum3.2 pkd_v_fac',chksum 02335 02336 02337 !CALL packed_vec_bcast(pkd_v, source, 'columnwise', mp_obj, error=error) 02338 !CALL packed_vec_bcast(pkd_u, source, 'columnwise', mp_obj, error=error) 02339 ! 02340 ! the core 02341 !v(k,1:n) = v(k,1:n) - uAk(i,k)' * inv(d(i,i)) * v(i,1:n) / s; i = 1,k-1 02342 !u(k,1:n) = u(k,1:n) - v(i,k)' * inv(d(i,i)) * u(i,1:n) / s; i = 1,k-1 02343 !old CALL packed_vec_bif_tech(mat_v, mat_li, pkd_v_fac, pkd_u_fac, k, pkd_v, pkd_u, error=error)!PARA local 02344 CALL packed_vec_bif_tech2(mat_v, pkd_v_fac, k, pkd_v, error=error)!PARA local 02345 CALL cpu_time(t(9)) 02346 t(55) = t(55) + t(9)-t(8) 02347 !chksum = cp_dbcsr_checksum(mat_v,error) 02348 !write(*,*) 'checksum4.0 mat_v',chksum 02349 ! bcast the v 02350 02351 02352 !write(*,'(A,1000F12.6)') '1pkd_v',pkd_v(nblkrows_total+2:int(pkd_v(nblkrows_total+1))-1) 02353 CALL packed_vec_bcast(pkd_v, -1, 'all', .TRUE. ,nblkrows_total, & 02354 row_blk_size(k)*row_blk_size, mp_obj, error) 02355 02356 02357 CALL cp_dbcsr_unpack_vec(mat_v, k, pkd_v, 'row', do_sum=.TRUE., error=error) !PARA local 02358 CALL cp_dbcsr_pack_vec(mat_v, k, pkd_v, 'row', error=error) !PARA local 02359 CALL packed_vec_bcast(pkd_v, -1, 'all', .FALSE. ,nblkrows_total, & 02360 row_blk_size(k)*row_blk_size, mp_obj, error) 02361 02362 !write(*,*) '2pkd_v',pkd_v( 1 : int(pkd_v(nblkrows_total+1)) - 1 ) 02363 02364 !chksum = sum(abs( pkd_v ( nblkrows_total+2 : int(pkd_v(nblkrows_total+1))-1 ) )) 02365 !write(*,'(A,1000F12.6)') '2pkd_v',pkd_v(nblkrows_total+2:int(pkd_v(nblkrows_total+1))-1) 02366 !write(*,*) 'checksum4.1',chksum 02367 ! 02368 ! compute norm 02369 ! 02370 ! filter the packed v 02371 ! 02372 ! build u(k,:) = [-v(k,1:k-1)/s 1 0 ... 0] 02373 CALL packed_vec_build_u(pkd_u, pkd_v, k, nblkrows_total, s, row_blk_size, error=error) 02374 CALL cpu_time(t(10)) 02375 t(56) = t(56) + t(10)-t(9) 02376 ! 02377 ! reduce the u and v 02378 !CALL packed_vec_reduce(pkd_v, to, 'columnwise', vec_blk_size, n, mp_obj, error=error) 02379 !CALL packed_vec_reduce(pkd_u, to, 'columnwise', vec_blk_size, n, mp_obj, error=error) 02380 ! 02381 ! unpack u(k,:) and v(k,:) 02382 CALL cp_dbcsr_unpack_vec(mat_v, k, pkd_v, 'row', error=error) !PARA local 02383 CALL cp_dbcsr_unpack_vec(mat_li, k, pkd_u, 'row', error=error) !PARA local 02384 !chksum = cp_dbcsr_checksum(mat_v,error) 02385 !write(*,*) 'checksum5 checksum(mat_v)',chksum 02386 !chksum = cp_dbcsr_checksum(mat_li,error) 02387 !write(*,*) 'checksum6.0 checksum(mat_li)',chksum 02388 !write(*,*) 'checksum6.1',sum(abs(pkd_u)) 02389 !write(*,*) 'pkd_u',pkd_u 02390 CALL cpu_time(t(11)) 02391 t(57) = t(57) + t(11)-t(10) 02392 !------------------------------------------------------------------------ 02393 ! build d and inv(d) 02394 ! 02395 ! d(k,k) = v(k,k) / s + eye(k); 02396 IF(proc_holds_diag_blk .EQ. dbcsr_mp_mynode (mp_obj)) THEN 02397 CALL cp_dbcsr_get_block(mat_v, k, k, blk, found) 02398 IF (.NOT.found) CALL stop_program(routineN,moduleN,__LINE__,& 02399 "Matrix block not found" ) 02400 CALL dscal(row_blk_size(k)**2, 1.0_dp/s, blk, 1) 02401 CALL block_add_on_diag(row_blk_size(k), blk, 1.0_dp) 02402 CALL dcopy(row_blk_size(k)**2, blk, 1, pkd_d(d_offset), 1) 02403 ! 02404 ! inv(d(k,k)) 02405 CALL block_chol_inv(row_blk_size(k), blk) 02406 CALL dcopy(row_blk_size(k)**2, blk, 1, pkd_inv_d(d_offset), 1) 02407 ENDIF 02408 ! 02409 ! bcast d(k,k) and inv(d(k,k)) 02410 CALL mp_bcast(pkd_d(d_offset:d_offset+row_blk_size(k)**2-1), & 02411 proc_holds_diag_blk, dbcsr_mp_group (mp_obj)) 02412 CALL mp_bcast(pkd_inv_d(d_offset:d_offset+row_blk_size(k)**2-1), & 02413 proc_holds_diag_blk, dbcsr_mp_group (mp_obj)) 02414 ! 02415 ! update the block pointer 02416 d_offset = d_offset + row_blk_size(k)**2 02417 !------------------------------------------------------------------------ 02418 CALL cpu_time(t(12)) 02419 t(58) = t(58) + t(12)-t(11) 02420 ENDDO 02421 CALL cp_dbcsr_finalize(mat_v, error=error) 02422 CALL cp_dbcsr_finalize(mat_li, error=error) 02423 02424 !WRITE(*,*) 'tin1',t(50) 02425 !WRITE(*,*) 'tin2',t(51) 02426 !WRITE(*,*) 'tin3',t(52) 02427 !WRITE(*,*) 'tin4',t(53) 02428 !WRITE(*,*) 'tin5',t(54) 02429 !WRITE(*,*) 'tin6',t(55) 02430 !WRITE(*,*) 'tin7',t(56) 02431 !WRITE(*,*) 'tin8',t(57) 02432 !WRITE(*,*) 'tin9',t(58) 02433 ! 02434 ! D=d; 02435 ! Linv=u; 02436 ! Lt=inv(D)*triu(v+s*I); 02437 CALL cp_dbcsr_init (mat_lt,error=error) 02438 CALL cp_dbcsr_btriu(mat_lt, mat_v, error=error) 02439 !chksum = cp_dbcsr_checksum(mat_lt,error) 02440 !WRITE(*,*) 'checksum(lt)_1',chksum 02441 CALL cp_dbcsr_add_on_diag(mat_lt, s, error=error) 02442 !chksum = cp_dbcsr_checksum(mat_lt,error) 02443 !WRITE(*,*) 'checksum(lt)_2',chksum,SUM(ABS(pkd_inv_d)) 02444 CALL dbcsr_scale_mat(mat_lt%matrix, alpha_matrix=pkd_inv_d, side='left', error=dbcsr_error) 02445 !chksum = cp_dbcsr_checksum(mat_lt,error) 02446 !WRITE(*,*) 'checksum(lt)_3',chksum 02447 CALL cp_dbcsr_init (mat_li_scaled,error=error) 02448 CALL cp_dbcsr_copy(mat_li_scaled, mat_li, error=error) 02449 CALL dbcsr_scale_mat(mat_li_scaled%matrix, alpha_matrix=pkd_inv_d, side='left', error=dbcsr_error) 02450 02451 CALL cp_dbcsr_init (mat_lt_scaled,error=error) 02452 CALL cp_dbcsr_copy(mat_lt_scaled, mat_lt, error=error) 02453 CALL dbcsr_scale_mat(mat_lt_scaled%matrix, alpha_matrix=pkd_d, side='left', error=dbcsr_error) 02454 02455 occ_in = cp_dbcsr_get_occupation(mat_lt) 02456 occ_out = cp_dbcsr_get_occupation(mat_li) 02457 WRITE(*,*) routineN//' occ(Lt)',occ_in,' occ(Li)',occ_out 02458 02459 chksum = cp_dbcsr_checksum(mat_v,error=error) 02460 IF(dbg_print_chksum)WRITE(*,*) 'checksum(v)',chksum 02461 chksum = cp_dbcsr_checksum(mat_li,error=error) 02462 IF(dbg_print_chksum)WRITE(*,*) 'checksum(li)',chksum 02463 chksum = cp_dbcsr_checksum(mat_li_scaled,error=error) 02464 IF(dbg_print_chksum)WRITE(*,*) 'checksum(li_scaled)',chksum 02465 chksum = cp_dbcsr_checksum(mat_lt,error=error) 02466 IF(dbg_print_chksum)WRITE(*,*) 'checksum(lt)',chksum 02467 02468 02469 IF(dbg_print_matrix)CALL cp_dbcsr_print(mat_v,error=error) 02470 IF(dbg_print_matrix)CALL cp_dbcsr_print(mat_li,error=error,matlab_format=.TRUE.) 02471 IF(dbg_print_matrix)CALL cp_dbcsr_print(mat_li_scaled,error=error,matlab_format=.TRUE.) 02472 IF(dbg_print_matrix)CALL cp_dbcsr_print(mat_lt,error=error) 02473 02474 ! 02475 ! At this point we should have inv(L), L^t, D and inv(D) 02476 ! compute the precond or apply the factors 02477 ! 02478 ! C = A*B' 02479 IF(.FALSE.)CALL cp_dbcsr_verify_matrix(prcnd, error=error) 02480 IF(.FALSE.)CALL cp_dbcsr_verify_matrix(mat_li, error=error) 02481 IF(.FALSE.)CALL cp_dbcsr_verify_matrix(mat_t, error=error) 02482 IF(.FALSE.)CALL cp_dbcsr_verify_matrix(mat_li_scaled, error=error) 02483 CALL cp_dbcsr_multiply("T", "N", 1.0_dp, mat_li, mat_li_scaled, 0.0_dp, prcnd, error=error) 02484 !CALL cp_dbcsr_print(prcnd, error, matlab_format=.TRUE.) 02485 02486 chksum = cp_dbcsr_checksum(prcnd,error=error) 02487 IF(dbg_print_chksum)WRITE(*,*) 'checksum(prcnd)',chksum 02488 02489 02490 ! 02491 ! checks 02492 IF(.TRUE.) THEN 02493 CALL cp_dbcsr_trace(mat_lt, mat_li, trace, error=error) 02494 WRITE(*,*) 'trace(Lt^t*Li)/n-1=',trace / REAL(SUM(row_blk_size),dp) - 1.0_dp 02495 02496 CALL cp_dbcsr_trace(mat_li, mat_lt, trace, error=error) 02497 WRITE(*,*) 'trace(Li^t*Lt)/n-1=',trace / REAL(SUM(row_blk_size),dp) - 1.0_dp 02498 02499 CALL cp_dbcsr_init(mat_tmp, error=error) 02500 CALL cp_dbcsr_create (mat_tmp, "tmp mat_t", cp_dbcsr_distribution (mat_t), dbcsr_type_no_symmetry,& 02501 cp_dbcsr_row_block_sizes (mat_t), cp_dbcsr_col_block_sizes (mat_t),& 02502 error=error) 02503 CALL cp_dbcsr_multiply("N", "N", 1.0_dp, mat_li_scaled, mat_t, 0.0_dp, mat_tmp, error=error) 02504 CALL cp_dbcsr_trace(mat_li, mat_tmp, trace, error=error) 02505 WRITE(*,*) 'trace(Li^t*A*Li*D)/n-1=',trace / REAL(SUM(row_blk_size),dp) - 1.0_dp 02506 02507 ! 02508 CALL cp_dbcsr_multiply("N", "N", 1.0_dp, prcnd, mat_t, 0.0_dp, mat_tmp, error=error) 02509 CALL cp_dbcsr_norm(mat_tmp, which_norm=dbcsr_norm_frobenius, norm_scalar=trace, error=error) 02510 WRITE(*,*) 'norm(A*Li*D*li^t)/n-1 = ',trace**2 / REAL(SUM(row_blk_size),dp) - 1.0_dp 02511 02512 CALL cp_dbcsr_multiply("T", "T", 1.0_dp, mat_t, prcnd, 0.0_dp, mat_tmp, error=error) 02513 CALL cp_dbcsr_norm(mat_tmp, which_norm=dbcsr_norm_frobenius, norm_scalar=trace, error=error) 02514 WRITE(*,*) 'norm(A*Li*D*li^t)/n-1 = ',trace**2 / REAL(SUM(row_blk_size),dp) - 1.0_dp 02515 02516 CALL cp_dbcsr_multiply("N", "N", 1.0_dp, mat_t, prcnd, 0.0_dp, mat_tmp, error=error) 02517 CALL cp_dbcsr_norm(mat_tmp, which_norm=dbcsr_norm_frobenius, norm_scalar=trace, error=error) 02518 WRITE(*,*) 'norm(Li*D*li^t*A)/n-1 = ',trace**2 / REAL(SUM(row_blk_size),dp) - 1.0_dp 02519 ! 02520 !norm(Lt'*D*Lt-A,'fro') 02521 CALL cp_dbcsr_copy(mat_tmp, mat_t, error=error) 02522 CALL cp_dbcsr_multiply("T", "N", 1.0_dp, mat_lt, mat_lt_scaled, -1.0_dp, mat_tmp, error=error) 02523 CALL cp_dbcsr_norm(mat_tmp, which_norm=dbcsr_norm_frobenius, norm_scalar=trace, error=error) 02524 WRITE(*,*) 'norm(Lt^t*D*Lt-A) = ',trace 02525 ! 02526 !norm(Lt*Linv'-I,'fro') 02527 CALL cp_dbcsr_multiply("N", "T", 1.0_dp, mat_lt, mat_li, 0.0_dp, mat_tmp, error=error) 02528 CALL cp_dbcsr_norm(mat_tmp, which_norm=dbcsr_norm_frobenius, norm_scalar=trace, error=error) 02529 WRITE(*,*) 'norm(Lt*Linv^t)/n-1 = ',trace**2 / REAL(SUM(row_blk_size),dp) - 1.0_dp 02530 02531 ENDIF 02532 ! 02533 ! cleanup 02534 DEALLOCATE(pkd_u, pkd_v, pkd_v_fac, pkd_u_fac, pkd_d, pkd_inv_d, pkd_v_tmp, blk, STAT=istat) 02535 CPPostcondition(istat==0,cp_failure_level,routineP,error,failure) 02536 02537 CALL cp_dbcsr_release(mat_t, error=error) 02538 CALL cp_dbcsr_release(mat_v, error=error) 02539 CALL cp_dbcsr_release(mat_lt, error=error) 02540 CALL cp_dbcsr_release(mat_li, error=error) 02541 CALL cp_dbcsr_release(mat_li_scaled, error=error) 02542 CALL cp_dbcsr_release(mat_lt_scaled, error=error) 02543 CALL cp_dbcsr_release(mat_tmp, error=error) 02544 02545 CALL cp_fm_release(fm_work1,error=error) 02546 CALL cp_fm_release(fm_work2,error=error) 02547 02548 CALL copy_dbcsr_to_fm (prcnd, fm,error=error)!for the moment 02549 CALL cp_dbcsr_deallocate_matrix(prcnd,error=error)!for the moment 02550 02551 CALL cpu_time(t(13)) 02552 02553 WRITE(*,*) 'tend',t(13)-t(12) 02554 02555 CALL timestop(handle) 02556 02557 END SUBROUTINE make_sparse_inverse_bif 02558 02559 SUBROUTINE apply_sparse_single(preconditioner_env, matrix_in, matrix_out,& 02560 error) 02561 02562 TYPE(preconditioner_type) :: preconditioner_env 02563 TYPE(cp_fm_type), POINTER :: matrix_in, matrix_out 02564 TYPE(cp_error_type), INTENT(inout) :: error 02565 02566 CHARACTER(len=*), PARAMETER :: routineN = 'apply_sparse_single', 02567 routineP = moduleN//':'//routineN 02568 02569 INTEGER :: handle, k 02570 TYPE(array_i1d_obj) :: col_blk_size, col_dist 02571 TYPE(cp_dbcsr_type) :: matrix_result, matrix_right 02572 TYPE(dbcsr_distribution_obj) :: dist 02573 02574 CALL timeset(routineN,handle) 02575 02576 WRITE(*,*)' APPLYING SPARSE SINGLE KINETIC' 02577 CALL cp_fm_get_info(matrix_in, ncol_global=k, error=error) 02578 CALL cp_create_bl_distribution (col_dist, col_blk_size, k, & 02579 dbcsr_mp_npcols(dbcsr_distribution_mp(cp_dbcsr_distribution( preconditioner_env%dbcsr_matrix )))) 02580 CALL dbcsr_distribution_new (dist, dbcsr_distribution_mp (cp_dbcsr_distribution( preconditioner_env%dbcsr_matrix )),& 02581 dbcsr_distribution_row_dist(cp_dbcsr_distribution( preconditioner_env%dbcsr_matrix )), col_dist) 02582 CALL cp_dbcsr_init (matrix_right, error) 02583 CALL cp_dbcsr_from_fm (matrix_right, matrix_in, 0.0_dp, dist,& 02584 row_blk_size = cp_dbcsr_row_block_sizes(preconditioner_env%dbcsr_matrix),& 02585 col_blk_size = col_blk_size, error=error) 02586 WRITE(*,*)'Pre checksum:',cp_dbcsr_checksum (matrix_right, error=error) 02587 CALL cp_dbcsr_init (matrix_result, error) 02588 CALL cp_dbcsr_create(matrix_result, "result", dist, 'N',& 02589 cp_dbcsr_row_block_sizes (preconditioner_env%dbcsr_matrix),& 02590 cp_dbcsr_col_block_sizes (matrix_right), error=error) 02591 CALL cp_dbcsr_multiply("N", "N", 1.0_dp, preconditioner_env%dbcsr_matrix,& 02592 matrix_right, 0.0_dp, matrix_result, error=error) 02593 CALL cp_dbcsr_distribution_release (dist) 02594 CALL cp_dbcsr_release (matrix_right, error=error) 02595 WRITE(*,*)'Post checksum:',cp_dbcsr_checksum (matrix_result, error=error) 02596 CALL cp_dbcsr_release (matrix_result, error=error) 02597 WRITE(*,*)' DONE' 02598 02599 CALL timestop(handle) 02600 02601 END SUBROUTINE apply_sparse_single 02602 02603 SUBROUTINE make_sparse_kinetic(preconditioner_env, matrix_t, matrix_s, & 02604 & energy_gap, error) 02605 02606 TYPE(preconditioner_type) :: preconditioner_env 02607 TYPE(cp_dbcsr_type), POINTER :: matrix_t, matrix_s 02608 REAL(KIND=dp), INTENT(IN) :: energy_gap 02609 TYPE(cp_error_type), INTENT(inout) :: error 02610 02611 CHARACTER(len=*), PARAMETER :: routineN = 'make_sparse_kinetic', 02612 routineP = moduleN//':'//routineN 02613 02614 INTEGER :: handle, n 02615 LOGICAL :: failure 02616 REAL(KIND=dp) :: shift 02617 02618 failure = .FALSE. 02619 CALL timeset(routineN,handle) 02620 CPPrecondition(ASSOCIATED(matrix_t),cp_failure_level,routineP,error,failure) 02621 CPPrecondition(ASSOCIATED(matrix_s),cp_failure_level,routineP,error,failure) 02622 ! some infos 02623 CALL cp_dbcsr_get_info(matrix_t,nfullrows_total=n) 02624 ! allocate 02625 02626 IF (ASSOCIATED(preconditioner_env%sparse_matrix)) THEN 02627 CALL cp_dbcsr_deallocate_matrix(preconditioner_env%sparse_matrix,error=error) 02628 NULLIFY(preconditioner_env%sparse_matrix) 02629 ENDIF 02630 02631 ALLOCATE(preconditioner_env%sparse_matrix) 02632 CALL cp_dbcsr_init(preconditioner_env%sparse_matrix,error=error) 02633 CALL cp_dbcsr_create(preconditioner_env%sparse_matrix, ' PRECONDITIONER ', & 02634 cp_dbcsr_distribution (matrix_t),& 02635 cp_dbcsr_get_matrix_type (matrix_t), cp_dbcsr_row_block_sizes(matrix_t),& 02636 cp_dbcsr_col_block_sizes(matrix_t), 0, 0, error=error) 02637 02638 ! M = T - epsilon * S 02639 shift=MAX(0.0_dp,energy_gap) 02640 CALL cp_dbcsr_add(preconditioner_env%sparse_matrix,matrix_t,& 02641 alpha_scalar=0.0_dp,beta_scalar=1.0_dp,error=error) 02642 CALL cp_dbcsr_add(preconditioner_env%sparse_matrix,matrix_s,& 02643 alpha_scalar=1.0_dp,beta_scalar=shift,error=error) 02644 CALL timestop(handle) 02645 02646 END SUBROUTINE make_sparse_kinetic 02647 02648 SUBROUTINE restart_preconditioner(qs_env,preconditioner,prec_type,nspins,error) 02649 02650 TYPE(qs_environment_type), POINTER :: qs_env 02651 TYPE(preconditioner_p_type), 02652 DIMENSION(:), POINTER :: preconditioner 02653 INTEGER, INTENT(IN) :: prec_type, nspins 02654 TYPE(cp_error_type), INTENT(inout) :: error 02655 02656 CHARACTER(LEN=*), PARAMETER :: routineN = 'restart_preconditioner', 02657 routineP = moduleN//':'//routineN 02658 02659 INTEGER :: ispin, stat 02660 LOGICAL :: failure 02661 02662 IF (ASSOCIATED(preconditioner)) THEN 02663 SELECT CASE(prec_type) 02664 CASE(ot_precond_full_all,ot_precond_full_single, ot_precond_full_single_inverse) ! these depend on the ks matrix 02665 DO ispin=1,SIZE(preconditioner) 02666 CALL destroy_preconditioner(preconditioner(ispin)%preconditioner,error=error) 02667 DEALLOCATE(preconditioner(ispin)%preconditioner) 02668 ENDDO 02669 DEALLOCATE(preconditioner,stat=stat) 02670 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 02671 CASE(ot_precond_none,ot_precond_full_kinetic,ot_precond_s_inverse, & 02672 ot_precond_sparse_diag,ot_precond_sparse_kinetic) ! these are 'independent' 02673 ! do nothing 02674 CASE DEFAULT 02675 CPPrecondition(.FALSE.,cp_failure_level,routineP,error,failure) 02676 END SELECT 02677 END IF 02678 02679 ! add an OT preconditioner if none is present 02680 IF (.NOT.ASSOCIATED(preconditioner)) THEN 02681 SELECT CASE(prec_type) 02682 CASE(ot_precond_full_all,ot_precond_full_single_inverse) 02683 ALLOCATE(preconditioner(nspins), stat=stat) 02684 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02685 CASE DEFAULT 02686 ALLOCATE(preconditioner(1), stat=stat) 02687 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 02688 END SELECT 02689 DO ispin=1,SIZE(preconditioner) 02690 ALLOCATE(preconditioner(ispin)%preconditioner) 02691 CALL init_preconditioner(preconditioner(ispin)%preconditioner,& 02692 para_env=qs_env%para_env,& 02693 blacs_env=qs_env%blacs_env,error=error) 02694 ENDDO 02695 END IF 02696 02697 END SUBROUTINE restart_preconditioner 02698 02699 SUBROUTINE prepare_preconditioner(qs_env,mos,matrix_ks,matrix_s,& 02700 ot_preconditioner,prec_type,solver_type,& 02701 energy_gap,nspins,has_unit_metric,mixed_precision,error) 02702 02703 TYPE(qs_environment_type), POINTER :: qs_env 02704 TYPE(mo_set_p_type), DIMENSION(:), 02705 POINTER :: mos 02706 TYPE(cp_dbcsr_p_type), DIMENSION(:), 02707 POINTER :: matrix_ks, matrix_s 02708 TYPE(preconditioner_p_type), 02709 DIMENSION(:), POINTER :: ot_preconditioner 02710 INTEGER, INTENT(IN) :: prec_type, solver_type 02711 REAL(dp), INTENT(IN) :: energy_gap 02712 INTEGER, INTENT(IN) :: nspins 02713 LOGICAL, INTENT(IN), OPTIONAL :: has_unit_metric, 02714 mixed_precision 02715 TYPE(cp_error_type), INTENT(inout) :: error 02716 02717 CHARACTER(LEN=*), PARAMETER :: routineN = 'prepare_preconditioner', 02718 routineP = moduleN//':'//routineN 02719 02720 CHARACTER(LEN=80) :: msg 02721 INTEGER :: handle, ispin 02722 LOGICAL :: do_co_rotate, do_rotation, failure, my_has_unit_metric, 02723 my_mixed_precision, use_mo_coeff_b 02724 TYPE(cp_dbcsr_type), POINTER :: matrix_t, mo_coeff_b 02725 TYPE(cp_fm_type), POINTER :: mo_coeff 02726 02727 CALL timeset(routineN,handle) 02728 02729 failure = .FALSE. 02730 my_has_unit_metric = .FALSE. 02731 IF(PRESENT(has_unit_metric)) my_has_unit_metric = has_unit_metric 02732 my_mixed_precision = .FALSE. 02733 IF(PRESENT(mixed_precision)) my_mixed_precision = mixed_precision 02734 02735 NULLIFY(matrix_t, mo_coeff_b, mo_coeff) 02736 02737 IF(qs_env%dft_control%qs_control%semi_empirical .OR. qs_env%dft_control%qs_control%dftb) THEN 02738 IF(prec_type==ot_precond_full_kinetic) THEN 02739 msg="Full_kinetic not available for semi-empirical methods" 02740 CPErrorMessage(cp_failure_level,routineP,TRIM(msg),error) 02741 END IF 02742 matrix_t => matrix_s(1)%matrix 02743 ELSE 02744 CPPrecondition(.NOT. my_has_unit_metric,cp_failure_level,routineP,error,failure) 02745 matrix_t => qs_env%kinetic(1)%matrix 02746 END IF 02747 02748 SELECT CASE(prec_type) 02749 CASE(ot_precond_none) 02750 DO ispin = 1,SIZE(ot_preconditioner) 02751 ot_preconditioner(ispin)%preconditioner%in_use=0 02752 END DO 02753 CASE(ot_precond_full_all,ot_precond_full_single_inverse) 02754 do_co_rotate = ASSOCIATED(qs_env%mo_derivs) 02755 02756 DO ispin=1,nspins 02757 CALL get_mo_set(mo_set=mos(ispin)%mo_set,mo_coeff_b=mo_coeff_b,mo_coeff=mo_coeff) 02758 use_mo_coeff_b =mos(ispin)%mo_set%use_mo_coeff_b 02759 IF(use_mo_coeff_b .AND. do_co_rotate) THEN 02760 CALL calculate_subspace_eigenvalues(mo_coeff_b,matrix_ks(ispin)%matrix,& 02761 do_rotation = .TRUE., & 02762 co_rotate=qs_env%mo_derivs(ispin)%matrix,& 02763 para_env=qs_env%para_env,& 02764 blacs_env=qs_env%blacs_env,error=error) 02765 ELSEIF(use_mo_coeff_b) THEN 02766 CALL calculate_subspace_eigenvalues(mo_coeff_b,matrix_ks(ispin)%matrix,& 02767 do_rotation = .TRUE., & 02768 para_env=qs_env%para_env,& 02769 blacs_env=qs_env%blacs_env,error=error) 02770 ELSE 02771 CALL calculate_subspace_eigenvalues(mo_coeff,matrix_ks(ispin)%matrix,& 02772 do_rotation = .TRUE., error=error) 02773 END IF 02774 IF(my_has_unit_metric) THEN 02775 CALL make_preconditioner(ot_preconditioner(ispin)%preconditioner, & 02776 prec_type, & 02777 solver_type, & 02778 matrix_h=matrix_ks(ispin)%matrix, & 02779 mo_set=mos(ispin)%mo_set, & 02780 energy_gap=energy_gap, & 02781 convert_precond_to_dbcsr=.TRUE.,& 02782 error=error) 02783 ELSE 02784 CALL make_preconditioner(ot_preconditioner(ispin)%preconditioner, & 02785 prec_type, & 02786 solver_type, & 02787 matrix_h=matrix_ks(ispin)%matrix,& 02788 matrix_s=matrix_s(1)%matrix,& 02789 matrix_t=matrix_t, & 02790 mo_set=mos(ispin)%mo_set,& 02791 energy_gap=energy_gap,& 02792 mixed_precision=my_mixed_precision,& 02793 convert_precond_to_dbcsr=.TRUE.,& 02794 error=error) 02795 END IF 02796 ENDDO 02797 02798 CASE DEFAULT 02799 IF(my_has_unit_metric) THEN 02800 CALL make_preconditioner(ot_preconditioner(1)%preconditioner, & 02801 prec_type, & 02802 solver_type, & 02803 matrix_h=matrix_ks(1)%matrix,& 02804 mo_set=mos(1)%mo_set,& 02805 energy_gap=energy_gap,& 02806 convert_precond_to_dbcsr=.TRUE.,& 02807 error=error) 02808 ELSE 02809 CALL make_preconditioner(ot_preconditioner(1)%preconditioner, & 02810 prec_type, & 02811 solver_type, & 02812 matrix_h=matrix_ks(1)%matrix,& 02813 matrix_s=matrix_s(1)%matrix,& 02814 matrix_t=matrix_t, & 02815 mo_set=mos(1)%mo_set,& 02816 energy_gap=energy_gap,& 02817 mixed_precision=my_mixed_precision,& 02818 convert_precond_to_dbcsr=.TRUE.,& 02819 error=error) 02820 END IF 02821 02822 END SELECT 02823 02824 CALL timestop(handle) 02825 02826 END SUBROUTINE prepare_preconditioner 02827 02828 END MODULE preconditioner 02829
1.7.3