1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief The methods which allow to analyze and manipulate the arnoldi procedure 8!> The main routine and this should eb the only public access point for the 9!> method 10!> \par History 11!> 2014.09 created [Florian Schiffmann] 12!> \author Florian Schiffmann 13! ************************************************************************************************** 14 15MODULE arnoldi_data_methods 16 USE arnoldi_types, ONLY: & 17 arnoldi_control_type, arnoldi_data_c_type, arnoldi_data_d_type, arnoldi_data_s_type, & 18 arnoldi_data_type, arnoldi_data_z_type, get_control, get_data_c, get_data_d, get_data_s, & 19 get_data_z, get_evals_c, get_evals_d, get_evals_s, get_evals_z, get_sel_ind, has_d_cmplx, & 20 has_d_real, has_s_cmplx, has_s_real, set_control, set_data_c, set_data_d, set_data_s, & 21 set_data_z 22 USE dbcsr_api, ONLY: & 23 dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_data_p, dbcsr_get_data_type, & 24 dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_mp_grid_setup, dbcsr_p_type, dbcsr_release, & 25 dbcsr_type, dbcsr_type_complex_4, dbcsr_type_complex_8, dbcsr_type_real_4, & 26 dbcsr_type_real_8, dbcsr_type_symmetric 27 USE dbcsr_vector, ONLY: create_col_vec_from_matrix 28 USE kinds, ONLY: real_4,& 29 real_8 30 USE util, ONLY: sort 31#include "../base/base_uses.f90" 32 33 IMPLICIT NONE 34 35 PRIVATE 36 37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_data_methods' 38 39 PUBLIC :: select_evals, get_selected_ritz_val, arnoldi_is_converged, & 40 arnoldi_data_type, get_nrestart, set_arnoldi_initial_vector, & 41 setup_arnoldi_data, deallocate_arnoldi_data, get_selected_ritz_vector 42 43CONTAINS 44 45! ************************************************************************************************** 46!> \brief This routine sets the environment for the arnoldi iteration and 47!> the krylov subspace creation. All simulation parameters have to be given 48!> at this stage so the rest can run fully automated 49!> In addition, this routine allocates the data necessary for 50!> \param arnoldi_data this type which gets filled with information and 51!> on output contains all information necessary to extract 52!> whatever the user desires 53!> \param matrix vector of matrices, only the first gets used to get some dimensions 54!> and parallel information needed later on 55!> \param max_iter maximum dimension of the krylov subspace 56!> \param threshold convergence threshold, this is used for both subspace and eigenval 57!> \param selection_crit integer defining according to which criterion the 58!> eigenvalues are selected for the subspace 59!> \param nval_request for some sel_crit useful, how many eV to select 60!> \param nrestarts ... 61!> \param generalized_ev ... 62!> \param iram ... 63! ************************************************************************************************** 64 SUBROUTINE setup_arnoldi_data(arnoldi_data, matrix, max_iter, threshold, selection_crit, & 65 nval_request, nrestarts, generalized_ev, iram) 66 TYPE(arnoldi_data_type) :: arnoldi_data 67 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 68 INTEGER :: max_iter 69 REAL(real_8) :: threshold 70 INTEGER :: selection_crit, nval_request, nrestarts 71 LOGICAL :: generalized_ev, iram 72 73 CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_arnoldi_data', & 74 routineP = moduleN//':'//routineN 75 76 CALL setup_arnoldi_control(arnoldi_data, matrix, max_iter, threshold, selection_crit, & 77 nval_request, nrestarts, generalized_ev, iram) 78 79 SELECT CASE (dbcsr_get_data_type (matrix (1)%matrix)) 80 CASE (dbcsr_type_real_8) 81 CALL setup_arnoldi_data_d(arnoldi_data, matrix, max_iter) 82 CASE (dbcsr_type_real_4) 83 CALL setup_arnoldi_data_s(arnoldi_data, matrix, max_iter) 84 CASE (dbcsr_type_complex_8) 85 CALL setup_arnoldi_data_z(arnoldi_data, matrix, max_iter) 86 CASE (dbcsr_type_complex_4) 87 CALL setup_arnoldi_data_c(arnoldi_data, matrix, max_iter) 88 END SELECT 89 90 END SUBROUTINE setup_arnoldi_data 91 92! ************************************************************************************************** 93!> \brief Creates the control type for arnoldi, see above for details 94!> \param arnoldi_data ... 95!> \param matrix ... 96!> \param max_iter ... 97!> \param threshold ... 98!> \param selection_crit ... 99!> \param nval_request ... 100!> \param nrestarts ... 101!> \param generalized_ev ... 102!> \param iram ... 103! ************************************************************************************************** 104 SUBROUTINE setup_arnoldi_control(arnoldi_data, matrix, max_iter, threshold, selection_crit, nval_request, & 105 nrestarts, generalized_ev, iram) 106 TYPE(arnoldi_data_type) :: arnoldi_data 107 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 108 INTEGER :: max_iter 109 REAL(real_8) :: threshold 110 INTEGER :: selection_crit, nval_request, nrestarts 111 LOGICAL :: generalized_ev, iram 112 113 CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_arnoldi_control', & 114 routineP = moduleN//':'//routineN 115 116 LOGICAL :: subgroups_defined 117 TYPE(arnoldi_control_type), POINTER :: control 118 TYPE(dbcsr_distribution_type) :: distri 119 120 ALLOCATE (control) 121! Fill the information which will later on control the behavior of the arnoldi method and allow synchronization 122 CALL dbcsr_get_info(matrix=matrix(1)%matrix, distribution=distri) 123 CALL dbcsr_mp_grid_setup(distri) 124 CALL dbcsr_distribution_get(distri, & 125 group=control%mp_group, & 126 mynode=control%myproc, & 127 subgroups_defined=subgroups_defined, & 128 prow_group=control%prow_group, & 129 pcol_group=control%pcol_group) 130 131 IF (.NOT. subgroups_defined) & 132 CPABORT("arnoldi only with subgroups") 133 134 control%symmetric = .FALSE. 135! Will need a fix for complex because there it has to be hermitian 136 IF (SIZE(matrix) == 1) & 137 control%symmetric = dbcsr_get_matrix_type(matrix(1)%matrix) == dbcsr_type_symmetric 138 139! Set the control parameters 140 control%max_iter = max_iter 141 control%current_step = 0 142 control%selection_crit = selection_crit 143 control%nval_req = nval_request 144 control%threshold = threshold 145 control%converged = .FALSE. 146 control%has_initial_vector = .FALSE. 147 control%iram = iram 148 control%nrestart = nrestarts 149 control%generalized_ev = generalized_ev 150 151 IF (control%nval_req > 1 .AND. control%nrestart > 0 .AND. .NOT. control%iram) & 152 CALL cp_abort(__LOCATION__, 'with more than one eigenvalue requested '// & 153 'internal restarting with a previous EVEC is a bad idea, set IRAM or nrestsart=0') 154 155! some checks for the generalized EV mode 156 IF (control%generalized_ev .AND. selection_crit == 1) & 157 CALL cp_abort(__LOCATION__, & 158 'generalized ev can only highest OR lowest EV') 159 IF (control%generalized_ev .AND. nval_request .NE. 1) & 160 CALL cp_abort(__LOCATION__, & 161 'generalized ev can only compute one EV at the time') 162 IF (control%generalized_ev .AND. control%nrestart == 0) & 163 CALL cp_abort(__LOCATION__, & 164 'outer loops are mandatory for generalized EV, set nrestart appropriatly') 165 IF (SIZE(matrix) .NE. 2 .AND. control%generalized_ev) & 166 CALL cp_abort(__LOCATION__, & 167 'generalized ev needs exactly two matrices as input (2nd is the metric)') 168 169 ALLOCATE (control%selected_ind(max_iter)) 170 CALL set_control(arnoldi_data, control) 171 172 END SUBROUTINE setup_arnoldi_control 173 174! ************************************************************************************************** 175!> \brief Deallocate the data in arnoldi_data 176!> \param arnoldi_data ... 177!> \param ind ... 178!> \param matrix ... 179!> \param vector ... 180! ************************************************************************************************** 181 SUBROUTINE get_selected_ritz_vector(arnoldi_data, ind, matrix, vector) 182 TYPE(arnoldi_data_type) :: arnoldi_data 183 INTEGER :: ind 184 TYPE(dbcsr_type) :: matrix, vector 185 186 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_selected_ritz_vector', & 187 routineP = moduleN//':'//routineN 188 189 IF (has_d_real(arnoldi_data)) CALL get_selected_ritz_vector_d(arnoldi_data, ind, matrix, vector) 190 IF (has_s_real(arnoldi_data)) CALL get_selected_ritz_vector_s(arnoldi_data, ind, matrix, vector) 191 IF (has_d_cmplx(arnoldi_data)) CALL get_selected_ritz_vector_z(arnoldi_data, ind, matrix, vector) 192 IF (has_s_cmplx(arnoldi_data)) CALL get_selected_ritz_vector_c(arnoldi_data, ind, matrix, vector) 193 194 END SUBROUTINE get_selected_ritz_vector 195 196! ************************************************************************************************** 197!> \brief Deallocate the data in arnoldi_data 198!> \param arnoldi_data ... 199! ************************************************************************************************** 200 SUBROUTINE deallocate_arnoldi_data(arnoldi_data) 201 TYPE(arnoldi_data_type) :: arnoldi_data 202 203 CHARACTER(LEN=*), PARAMETER :: routineN = 'deallocate_arnoldi_data', & 204 routineP = moduleN//':'//routineN 205 206 TYPE(arnoldi_control_type), POINTER :: control 207 208 IF (has_d_real(arnoldi_data)) CALL deallocate_arnoldi_data_d(arnoldi_data) 209 IF (has_s_real(arnoldi_data)) CALL deallocate_arnoldi_data_s(arnoldi_data) 210 IF (has_d_cmplx(arnoldi_data)) CALL deallocate_arnoldi_data_z(arnoldi_data) 211 IF (has_s_cmplx(arnoldi_data)) CALL deallocate_arnoldi_data_c(arnoldi_data) 212 213 control => get_control(arnoldi_data) 214 DEALLOCATE (control%selected_ind) 215 DEALLOCATE (control) 216 217 END SUBROUTINE deallocate_arnoldi_data 218 219! ************************************************************************************************** 220!> \brief perform the selection of eigenvalues, fills the selected_ind array 221!> \param arnoldi_data ... 222! ************************************************************************************************** 223 SUBROUTINE select_evals(arnoldi_data) 224 TYPE(arnoldi_data_type) :: arnoldi_data 225 226 CHARACTER(LEN=*), PARAMETER :: routineN = 'select_evals', routineP = moduleN//':'//routineN 227 228 IF (has_d_real(arnoldi_data)) CALL select_evals_d(arnoldi_data) 229 IF (has_s_real(arnoldi_data)) CALL select_evals_s(arnoldi_data) 230 IF (has_d_cmplx(arnoldi_data)) CALL select_evals_z(arnoldi_data) 231 IF (has_s_cmplx(arnoldi_data)) CALL select_evals_c(arnoldi_data) 232 233 END SUBROUTINE select_evals 234 235! ************************************************************************************************** 236!> \brief set a new selection type, if you notice you didn't like the initial one 237!> \param ar_data ... 238!> \param itype ... 239! ************************************************************************************************** 240 SUBROUTINE set_eval_selection(ar_data, itype) 241 TYPE(arnoldi_data_type) :: ar_data 242 INTEGER :: itype 243 244 TYPE(arnoldi_control_type), POINTER :: control 245 246 control => get_control(ar_data) 247 control%selection_crit = itype 248 END SUBROUTINE set_eval_selection 249 250! ************************************************************************************************** 251!> \brief returns the number of restarts allowed for arnoldi 252!> \param arnoldi_data ... 253!> \return ... 254! ************************************************************************************************** 255 FUNCTION get_nrestart(arnoldi_data) RESULT(nrestart) 256 TYPE(arnoldi_data_type) :: arnoldi_data 257 INTEGER :: nrestart 258 259 TYPE(arnoldi_control_type), POINTER :: control 260 261 control => get_control(arnoldi_data) 262 nrestart = control%nrestart 263 264 END FUNCTION get_nrestart 265 266! ************************************************************************************************** 267!> \brief get the number of eigenvalues matching the search criterion 268!> \param ar_data ... 269!> \return ... 270! ************************************************************************************************** 271 FUNCTION get_nval_out(ar_data) RESULT(nval_out) 272 TYPE(arnoldi_data_type) :: ar_data 273 INTEGER :: nval_out 274 275 TYPE(arnoldi_control_type), POINTER :: control 276 277 control => get_control(ar_data) 278 nval_out = control%nval_out 279 280 END FUNCTION get_nval_out 281 282! ************************************************************************************************** 283!> \brief get the dimension of the krylov space. Can be less than max_iter 284!> if subspace converged early 285!> \param ar_data ... 286!> \return ... 287! ************************************************************************************************** 288 FUNCTION get_subsp_size(ar_data) RESULT(current_step) 289 TYPE(arnoldi_data_type) :: ar_data 290 INTEGER :: current_step 291 292 TYPE(arnoldi_control_type), POINTER :: control 293 294 control => get_control(ar_data) 295 current_step = control%current_step 296 297 END FUNCTION get_subsp_size 298 299! ************************************************************************************************** 300!> \brief Find out whether the method with the current search criterion is converged 301!> \param ar_data ... 302!> \return ... 303! ************************************************************************************************** 304 FUNCTION arnoldi_is_converged(ar_data) RESULT(converged) 305 TYPE(arnoldi_data_type) :: ar_data 306 LOGICAL :: converged 307 308 TYPE(arnoldi_control_type), POINTER :: control 309 310 control => get_control(ar_data) 311 converged = control%converged 312 313 END FUNCTION 314 315! ************************************************************************************************** 316!> \brief get a single specific Ritz value from the set of selected 317!> \param ar_data ... 318!> \param ind ... 319!> \return ... 320! ************************************************************************************************** 321 FUNCTION get_selected_ritz_val(ar_data, ind) RESULT(eval_out) 322 TYPE(arnoldi_data_type) :: ar_data 323 INTEGER :: ind 324 COMPLEX(real_8) :: eval_out 325 326 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_selected_ritz_val', & 327 routineP = moduleN//':'//routineN 328 329 COMPLEX(real_4), DIMENSION(:), POINTER :: evals_s 330 COMPLEX(real_8), DIMENSION(:), POINTER :: evals_d 331 INTEGER :: ev_ind 332 INTEGER, DIMENSION(:), POINTER :: selected_ind 333 334 IF (ind .GT. get_nval_out(ar_data)) & 335 CPABORT('outside range of indexed evals') 336 337 selected_ind => get_sel_ind(ar_data) 338 ev_ind = selected_ind(ind) 339 IF (has_d_real(ar_data)) THEN 340 evals_d => get_evals_d(ar_data); eval_out = evals_d(ev_ind) 341 ELSE IF (has_s_real(ar_data)) THEN 342 evals_s => get_evals_s(ar_data); eval_out = CMPLX(evals_s(ev_ind), KIND=real_8) 343 ELSE IF (has_d_cmplx(ar_data)) THEN 344 evals_d => get_evals_z(ar_data); eval_out = evals_d(ev_ind) 345 ELSE IF (has_s_cmplx(ar_data)) THEN 346 evals_s => get_evals_c(ar_data); eval_out = CMPLX(evals_s(ev_ind), KIND=real_8) 347 END IF 348 349 END FUNCTION get_selected_ritz_val 350 351! ************************************************************************************************** 352!> \brief Get all Ritz values of the selected set. eval_out has to be allocated 353!> at least the size of get_neval_out() 354!> \param ar_data ... 355!> \param eval_out ... 356! ************************************************************************************************** 357 SUBROUTINE get_all_selected_ritz_val(ar_data, eval_out) 358 TYPE(arnoldi_data_type) :: ar_data 359 COMPLEX(real_8), DIMENSION(:) :: eval_out 360 361 CHARACTER(LEN=*), PARAMETER :: routineN = 'get_all_selected_ritz_val', & 362 routineP = moduleN//':'//routineN 363 364 COMPLEX(real_4), DIMENSION(:), POINTER :: evals_s 365 COMPLEX(real_8), DIMENSION(:), POINTER :: evals_d 366 INTEGER :: ev_ind, ind 367 INTEGER, DIMENSION(:), POINTER :: selected_ind 368 369 NULLIFY (evals_d, evals_s) 370 IF (SIZE(eval_out) .LT. get_nval_out(ar_data)) & 371 CPABORT('array for eval output too small') 372 selected_ind => get_sel_ind(ar_data) 373 374 IF (has_d_real(ar_data)) evals_d => get_evals_d(ar_data) 375 IF (has_s_real(ar_data)) evals_s => get_evals_s(ar_data) 376 IF (has_d_cmplx(ar_data)) evals_d => get_evals_d(ar_data) 377 IF (has_s_cmplx(ar_data)) evals_s => get_evals_c(ar_data) 378 379 DO ind = 1, get_nval_out(ar_data) 380 ev_ind = selected_ind(ind) 381 IF (ASSOCIATED(evals_d)) eval_out(ind) = evals_d(ev_ind) 382 IF (ASSOCIATED(evals_s)) eval_out(ind) = CMPLX(evals_s(ev_ind), KIND=real_8) 383 END DO 384 385 END SUBROUTINE get_all_selected_ritz_val 386 387! ************************************************************************************************** 388!> \brief ... 389!> \param ar_data ... 390!> \param vector ... 391! ************************************************************************************************** 392 SUBROUTINE set_arnoldi_initial_vector(ar_data, vector) 393 TYPE(arnoldi_data_type) :: ar_data 394 TYPE(dbcsr_type) :: vector 395 396 TYPE(arnoldi_control_type), POINTER :: control 397 398 control => get_control(ar_data) 399 control%has_initial_vector = .TRUE. 400 401 IF (has_d_real(ar_data)) CALL set_initial_vector_d(ar_data, vector) 402 IF (has_s_real(ar_data)) CALL set_initial_vector_s(ar_data, vector) 403 IF (has_d_cmplx(ar_data)) CALL set_initial_vector_z(ar_data, vector) 404 IF (has_s_cmplx(ar_data)) CALL set_initial_vector_c(ar_data, vector) 405 406 END SUBROUTINE set_arnoldi_initial_vector 407 408#:include 'arnoldi_data_selection.f90' 409#:include 'arnoldi_data_manipulation.f90' 410END MODULE arnoldi_data_methods 411 412