1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief arnoldi iteration using dbcsr 8!> \par History 9!> 2014.09 created [Florian Schiffmann] 10!> \author Florian Schiffmann 11! ************************************************************************************************** 12 13MODULE arnoldi_api 14 USE arnoldi_data_methods, ONLY: arnoldi_is_converged,& 15 deallocate_arnoldi_data,& 16 get_nrestart,& 17 get_selected_ritz_val,& 18 get_selected_ritz_vector,& 19 select_evals,& 20 set_arnoldi_initial_vector,& 21 setup_arnoldi_data 22 USE arnoldi_methods, ONLY: arnoldi_init,& 23 arnoldi_iram,& 24 build_subspace,& 25 compute_evals,& 26 gev_arnoldi_init,& 27 gev_build_subspace,& 28 gev_update_data 29 USE arnoldi_types, ONLY: arnoldi_control_type,& 30 arnoldi_data_type,& 31 get_control,& 32 m_x_v_vectors_type 33 USE dbcsr_api, ONLY: & 34 dbcsr_add, dbcsr_copy, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, & 35 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 36 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type 37 USE dbcsr_vector, ONLY: create_col_vec_from_matrix,& 38 create_replicated_col_vec_from_matrix,& 39 create_replicated_row_vec_from_matrix,& 40 dbcsr_matrix_colvec_multiply 41 USE kinds, ONLY: dp 42 USE message_passing, ONLY: mp_sum 43#include "../base/base_uses.f90" 44 45 IMPLICIT NONE 46 47 PRIVATE 48 49 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_api' 50 51 PUBLIC :: arnoldi_ev, arnoldi_extremal, arnoldi_conjugate_gradient 52 PUBLIC :: arnoldi_data_type, setup_arnoldi_data, deallocate_arnoldi_data 53 PUBLIC :: set_arnoldi_initial_vector 54 PUBLIC :: get_selected_ritz_val, get_selected_ritz_vector 55 56CONTAINS 57 58! ************************************************************************************************** 59!> \brief Driver routine for different arnoldi eigenvalue methods 60!> the selection which one is to be taken is made beforehand in the 61!> setup call passing the generalized_ev flag or not 62!> \param matrix ... 63!> \param arnoldi_data ... 64! ************************************************************************************************** 65 66 SUBROUTINE arnoldi_ev(matrix, arnoldi_data) 67 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 68 TYPE(arnoldi_data_type) :: arnoldi_data 69 70 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_ev', routineP = moduleN//':'//routineN 71 72 TYPE(arnoldi_control_type), POINTER :: control 73 74 control => get_control(arnoldi_data) 75 76 IF (control%generalized_ev) THEN 77 CALL arnoldi_generalized_ev(matrix, arnoldi_data) 78 ELSE 79 CALL arnoldi_normal_ev(matrix, arnoldi_data) 80 END IF 81 82 END SUBROUTINE arnoldi_ev 83 84! ************************************************************************************************** 85!> \brief The main routine for arnoldi method to compute ritz values 86!> vectors of a matrix. Can take multiple matrices to solve 87!> ( M(N)*...*M(2)*M(1) )*v=v*e. A, B, ... have to be merged in a array of pointers 88!> arnoldi data has to be create with the setup routine and 89!> will contain on input all necessary information to start/restart 90!> the calculation. On output it contains all data 91!> \param matrix a pointer array to dbcsr_matrices. Multiplication order is as 92!> described above 93!> \param arnoldi_data On input data_type contains all information to start/restart 94!> an arnoldi iteration 95!> On output all data areas are filled to allow arbitrary post 96!> processing of the created subspace 97!> arnoldi_data has to be created with setup_arnoldi_data 98! ************************************************************************************************** 99 SUBROUTINE arnoldi_normal_ev(matrix, arnoldi_data) 100 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 101 TYPE(arnoldi_data_type) :: arnoldi_data 102 103 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_normal_ev', & 104 routineP = moduleN//':'//routineN 105 106 INTEGER :: handle, i_loop, ncol_local, nrow_local 107 TYPE(arnoldi_control_type), POINTER :: control 108 TYPE(dbcsr_type), POINTER :: restart_vec 109 TYPE(m_x_v_vectors_type) :: vectors 110 111 NULLIFY (restart_vec) 112 CALL timeset(routineN, handle) 113 114!prepare the vector like matrives needed in the matrix vector products, they will be reused throughout the iterations 115 CALL create_col_vec_from_matrix(vectors%input_vec, matrix(1)%matrix, 1) 116 CALL dbcsr_copy(vectors%result_vec, vectors%input_vec) 117 CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix(1)%matrix, 1) 118 CALL create_replicated_row_vec_from_matrix(vectors%rep_row_vec, matrix(1)%matrix, 1) 119 120! Tells whether we have local data available on the processor (usually all in pcol 0 but even ther can be some without data) 121 control => get_control(arnoldi_data) 122 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local) 123 control%local_comp = ncol_local > 0 .AND. nrow_local > 0 124 125 DO i_loop = 0, get_nrestart(arnoldi_data) 126 127 IF (.NOT. control%iram .OR. i_loop == 0) THEN 128! perform the standard arnoldi, if restarts are requested use the first (only makes sense if 1ev is requested) 129 IF (ASSOCIATED(restart_vec)) CALL set_arnoldi_initial_vector(arnoldi_data, restart_vec) 130 CALL arnoldi_init(matrix, vectors, arnoldi_data) 131 ELSE 132! perform an implicit restart 133 CALL arnoldi_iram(arnoldi_data) 134 END IF 135 136! Generate the subspace 137 CALL build_subspace(matrix, vectors, arnoldi_data) 138 139! If we reached the maximum number of steps or the subspace converged we still need to get the eigenvalues 140 CALL compute_evals(arnoldi_data) 141 142! Select the evals according to user selection and keep them in arnoldi_data 143 CALL select_evals(arnoldi_data) 144 145! Prepare for a restart with the best eigenvector not needed in case of iram but who cares 146 IF (.NOT. ASSOCIATED(restart_vec)) ALLOCATE (restart_vec) 147 CALL get_selected_ritz_vector(arnoldi_data, 1, matrix(1)%matrix, restart_vec) 148 149! Check whether we can already go home 150 IF (control%converged) EXIT 151 END DO 152 153! Deallocated the work vectors 154 CALL dbcsr_release(vectors%input_vec) 155 CALL dbcsr_release(vectors%result_vec) 156 CALL dbcsr_release(vectors%rep_col_vec) 157 CALL dbcsr_release(vectors%rep_row_vec) 158 CALL dbcsr_release(restart_vec) 159 DEALLOCATE (restart_vec) 160 CALL timestop(handle) 161 162 END SUBROUTINE arnoldi_normal_ev 163 164! ************************************************************************************************** 165!> \brief The main routine for arnoldi method to compute the lowest ritz pair 166!> of a symmetric generalized eigenvalue problem. 167!> as input it takes a vector of matrices which for the GEV: 168!> M(1)*v=M(2)*v*lambda 169!> In other words, M(1) is the matrix and M(2) the metric 170!> This only works if the two matrices are symmetric in values 171!> (flag in dbcsr does not need to be set) 172!> \param matrix a pointer array to dbcsr_matrices. Multiplication order is as 173!> described above 174!> \param arnoldi_data On input data_type contains all information to start/restart 175!> an arnoldi iteration 176!> On output all data areas are filled to allow arbitrary post 177!> processing of the created subspace 178!> arnoldi_data has to be created with setup_arnoldi_data 179! ************************************************************************************************** 180 SUBROUTINE arnoldi_generalized_ev(matrix, arnoldi_data) 181 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 182 TYPE(arnoldi_data_type) :: arnoldi_data 183 184 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_generalized_ev', & 185 routineP = moduleN//':'//routineN 186 187 INTEGER :: handle, i_loop, ncol_local, nrow_local 188 TYPE(arnoldi_control_type), POINTER :: control 189 TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: matrix_arnoldi 190 TYPE(dbcsr_type), TARGET :: A_rho_B 191 TYPE(m_x_v_vectors_type) :: vectors 192 193 CALL timeset(routineN, handle) 194 ALLOCATE (matrix_arnoldi(2)) 195 ! this matrix will contain +/- A-rho*B 196 matrix_arnoldi(1)%matrix => A_rho_B 197 matrix_arnoldi(2)%matrix => matrix(2)%matrix 198 199!prepare the vector like matrives needed in the matrix vector products, they will be reused throughout the iterations 200 CALL create_col_vec_from_matrix(vectors%input_vec, matrix(1)%matrix, 1) 201 CALL dbcsr_copy(vectors%result_vec, vectors%input_vec) 202 CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix(1)%matrix, 1) 203 CALL create_replicated_row_vec_from_matrix(vectors%rep_row_vec, matrix(1)%matrix, 1) 204 205! Tells whether we have local data available on the processor (usually all in pcol 0 but even ther can be some without data) 206 control => get_control(arnoldi_data) 207 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local) 208 control%local_comp = ncol_local > 0 .AND. nrow_local > 0 209 210 DO i_loop = 0, get_nrestart(arnoldi_data) 211 IF (i_loop == 0) THEN 212! perform the standard arnoldi initialization with a random vector 213 CALL gev_arnoldi_init(matrix, matrix_arnoldi, vectors, arnoldi_data) 214 END IF 215 216! Generate the subspace 217 CALL gev_build_subspace(matrix_arnoldi, vectors, arnoldi_data) 218 219! If we reached the maximum number of steps or the subspace converged we still need to get the eigenvalues 220 CALL compute_evals(arnoldi_data) 221 222! Select the evals according to user selection and keep them in arnoldi_data 223 CALL select_evals(arnoldi_data) 224 225! update the matrices and compute the convergence 226 CALL gev_update_data(matrix, matrix_arnoldi, vectors, arnoldi_data) 227 228! Check whether we can already go home 229 IF (control%converged) EXIT 230 END DO 231 232! Deallocated the work vectors 233 CALL dbcsr_release(vectors%input_vec) 234 CALL dbcsr_release(vectors%result_vec) 235 CALL dbcsr_release(vectors%rep_col_vec) 236 CALL dbcsr_release(vectors%rep_row_vec) 237 CALL dbcsr_release(A_rho_B) 238 DEALLOCATE (matrix_arnoldi) 239 240 CALL timestop(handle) 241 242 END SUBROUTINE arnoldi_generalized_ev 243 244! ************************************************************************************************** 245!> \brief simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface 246!> this hides some of the power of the arnoldi routines (e.g. only min or max eval, generalized eval, ritz vectors, etc.), 247!> and does not allow for providing an initial guess of the ritz vector. 248!> \param matrix_a input mat 249!> \param max_ev estimated max eval 250!> \param min_ev estimated min eval 251!> \param converged Usually arnoldi is more accurate than claimed. 252!> \param threshold target precision 253!> \param max_iter max allowed iterations (will be rounded up) 254! ************************************************************************************************** 255 SUBROUTINE arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter) 256 TYPE(dbcsr_type), INTENT(INOUT), TARGET :: matrix_a 257 REAL(KIND=dp), INTENT(OUT) :: max_ev, min_ev 258 LOGICAL, INTENT(OUT) :: converged 259 REAL(KIND=dp), INTENT(IN) :: threshold 260 INTEGER, INTENT(IN) :: max_iter 261 262 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_extremal', & 263 routineP = moduleN//':'//routineN 264 265 INTEGER :: handle, max_iter_internal, nrestarts 266 TYPE(arnoldi_data_type) :: my_arnoldi 267 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices 268 269 CALL timeset(routineN, handle) 270 271 ! we go in chunks of max_iter_internal, and restart ater each of those. 272 ! at low threshold smaller values of max_iter_internal make sense 273 IF (.TRUE.) max_iter_internal = 16 274 IF (threshold <= 1.0E-3_dp) max_iter_internal = 32 275 IF (threshold <= 1.0E-4_dp) max_iter_internal = 64 276 277 ! the max number of iter will be (nrestarts+1)*max_iter_internal 278 nrestarts = max_iter/max_iter_internal 279 280 ALLOCATE (arnoldi_matrices(1)) 281 arnoldi_matrices(1)%matrix => matrix_a 282 CALL setup_arnoldi_data(my_arnoldi, arnoldi_matrices, max_iter=max_iter_internal, & 283 threshold=threshold, selection_crit=1, nval_request=2, nrestarts=nrestarts, & 284 generalized_ev=.FALSE., iram=.TRUE.) 285 CALL arnoldi_ev(arnoldi_matrices, my_arnoldi) 286 converged = arnoldi_is_converged(my_arnoldi) 287 max_eV = REAL(get_selected_ritz_val(my_arnoldi, 2), dp) 288 min_eV = REAL(get_selected_ritz_val(my_arnoldi, 1), dp) 289 CALL deallocate_arnoldi_data(my_arnoldi) 290 DEALLOCATE (arnoldi_matrices) 291 292 CALL timestop(handle) 293 294 END SUBROUTINE arnoldi_extremal 295 296! ************************************************************************************************** 297!> \brief Wrapper for conjugated gradient algorithm for Ax=b 298!> \param matrix_a input mat 299!> \param vec_x input right hand side vector; output solution vector, fully replicated! 300!> \param matrix_p input preconditioner (optional) 301!> \param converged ... 302!> \param threshold target precision 303!> \param max_iter max allowed iterations 304! ************************************************************************************************** 305 SUBROUTINE arnoldi_conjugate_gradient(matrix_a, vec_x, matrix_p, converged, threshold, max_iter) 306 TYPE(dbcsr_type), INTENT(IN), TARGET :: matrix_a 307 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: vec_x 308 TYPE(dbcsr_type), INTENT(IN), OPTIONAL, TARGET :: matrix_p 309 LOGICAL, INTENT(OUT) :: converged 310 REAL(KIND=dp), INTENT(IN) :: threshold 311 INTEGER, INTENT(IN) :: max_iter 312 313 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_conjugate_gradient', & 314 routineP = moduleN//':'//routineN 315 316 INTEGER :: handle, i, j, mpgrp, nb, nloc, no 317 INTEGER, DIMENSION(:), POINTER :: rb_offset, rb_size 318 REAL(KIND=dp), DIMENSION(:, :), POINTER :: xvec 319 TYPE(arnoldi_control_type), POINTER :: control 320 TYPE(arnoldi_data_type) :: my_arnoldi 321 TYPE(dbcsr_iterator_type) :: dbcsr_iter 322 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: arnoldi_matrices 323 TYPE(dbcsr_type), TARGET :: x 324 TYPE(m_x_v_vectors_type) :: vectors 325 326 CALL timeset(routineN, handle) 327 328 !prepare the vector like matrices needed in the matrix vector products, 329 !they will be reused throughout the iterations 330 CALL create_col_vec_from_matrix(vectors%input_vec, matrix_a, 1) 331 CALL dbcsr_copy(vectors%result_vec, vectors%input_vec) 332 CALL create_replicated_col_vec_from_matrix(vectors%rep_col_vec, matrix_a, 1) 333 CALL create_replicated_row_vec_from_matrix(vectors%rep_row_vec, matrix_a, 1) 334 ! 335 CALL dbcsr_copy(x, vectors%input_vec) 336 ! 337 CALL dbcsr_get_info(x, nfullrows_local=nloc, row_blk_size=rb_size, row_blk_offset=rb_offset) 338 ! 339 CALL dbcsr_iterator_start(dbcsr_iter, x) 340 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 341 CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, xvec) 342 nb = rb_size(i) 343 no = rb_offset(i) 344 xvec(1:nb, 1) = vec_x(no:no + nb - 1) 345 END DO 346 CALL dbcsr_iterator_stop(dbcsr_iter) 347 348 ! Arnoldi interface (used just for the iterator information here 349 ALLOCATE (arnoldi_matrices(3)) 350 arnoldi_matrices(1)%matrix => matrix_a 351 IF (PRESENT(matrix_p)) THEN 352 arnoldi_matrices(2)%matrix => matrix_p 353 ELSE 354 NULLIFY (arnoldi_matrices(2)%matrix) 355 END IF 356 arnoldi_matrices(3)%matrix => x 357 CALL setup_arnoldi_data(my_arnoldi, arnoldi_matrices, max_iter=max_iter, & 358 threshold=threshold, selection_crit=1, nval_request=1, nrestarts=0, & 359 generalized_ev=.FALSE., iram=.FALSE.) 360 361 CALL conjugate_gradient(my_arnoldi, arnoldi_matrices, vectors) 362 363 converged = arnoldi_is_converged(my_arnoldi) 364 365 vec_x = 0.0_dp 366 CALL dbcsr_iterator_start(dbcsr_iter, x) 367 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 368 CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, xvec) 369 nb = rb_size(i) 370 no = rb_offset(i) 371 vec_x(no:no + nb - 1) = xvec(1:nb, 1) 372 END DO 373 CALL dbcsr_iterator_stop(dbcsr_iter) 374 control => get_control(my_arnoldi) 375 mpgrp = control%mp_group 376 CALL mp_sum(vec_x, mpgrp) 377 ! Deallocated the work vectors 378 CALL dbcsr_release(x) 379 CALL dbcsr_release(vectors%input_vec) 380 CALL dbcsr_release(vectors%result_vec) 381 CALL dbcsr_release(vectors%rep_col_vec) 382 CALL dbcsr_release(vectors%rep_row_vec) 383 384 CALL deallocate_arnoldi_data(my_arnoldi) 385 DEALLOCATE (arnoldi_matrices) 386 387 CALL timestop(handle) 388 389 END SUBROUTINE arnoldi_conjugate_gradient 390 391! ************************************************************************************************** 392!> \brief ... 393!> \param arnoldi_data ... 394!> \param arnoldi_matrices ... 395!> \param vectors ... 396! ************************************************************************************************** 397 SUBROUTINE conjugate_gradient(arnoldi_data, arnoldi_matrices, vectors) 398 TYPE(arnoldi_data_type) :: arnoldi_data 399 TYPE(dbcsr_p_type), DIMENSION(:) :: arnoldi_matrices 400 TYPE(m_x_v_vectors_type) :: vectors 401 402 CHARACTER(LEN=*), PARAMETER :: routineN = 'conjugate_gradient', & 403 routineP = moduleN//':'//routineN 404 405 INTEGER :: iter, mpgrp, pcgrp 406 REAL(KIND=dp) :: alpha, beta, pap, rsnew, rsold 407 TYPE(arnoldi_control_type), POINTER :: control 408 TYPE(dbcsr_type) :: apvec, pvec, rvec, zvec 409 TYPE(dbcsr_type), POINTER :: amat, pmat, xvec 410 411 control => get_control(arnoldi_data) 412 control%converged = .FALSE. 413 pcgrp = control%pcol_group 414 mpgrp = control%mp_group 415 416 NULLIFY (amat, pmat, xvec) 417 amat => arnoldi_matrices(1)%matrix 418 pmat => arnoldi_matrices(2)%matrix 419 xvec => arnoldi_matrices(3)%matrix 420 421 IF (ASSOCIATED(pmat)) THEN 422 ! Preconditioned conjugate gradients 423 CALL dbcsr_copy(vectors%input_vec, xvec) 424 CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, & 425 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec) 426 CALL dbcsr_copy(pvec, vectors%result_vec) 427 CALL dbcsr_copy(vectors%input_vec, pvec) 428 CALL dbcsr_matrix_colvec_multiply(amat, vectors%input_vec, vectors%result_vec, 1.0_dp, & 429 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec) 430 CALL dbcsr_copy(apvec, vectors%result_vec) 431 CALL dbcsr_copy(rvec, xvec) 432 CALL dbcsr_add(rvec, apvec, 1.0_dp, -1.0_dp) 433 CALL dbcsr_copy(xvec, pvec) 434 CALL dbcsr_copy(vectors%input_vec, rvec) 435 CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, & 436 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec) 437 CALL dbcsr_copy(zvec, vectors%result_vec) 438 CALL dbcsr_copy(pvec, zvec) 439 rsold = vec_dot_vec(rvec, zvec, mpgrp) 440 DO iter = 1, control%max_iter 441 CALL dbcsr_copy(vectors%input_vec, pvec) 442 CALL dbcsr_matrix_colvec_multiply(amat, vectors%input_vec, vectors%result_vec, 1.0_dp, & 443 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec) 444 CALL dbcsr_copy(apvec, vectors%result_vec) 445 446 pap = vec_dot_vec(pvec, apvec, mpgrp) 447 IF (ABS(pap) < 1.e-24_dp) THEN 448 alpha = 0.0_dp 449 ELSE 450 alpha = rsold/pap 451 END IF 452 453 CALL dbcsr_add(xvec, pvec, 1.0_dp, alpha) 454 CALL dbcsr_add(rvec, apvec, 1.0_dp, -alpha) 455 rsnew = vec_dot_vec(rvec, rvec, mpgrp) 456 IF (SQRT(rsnew) < control%threshold) EXIT 457 CPASSERT(alpha /= 0.0_dp) 458 459 CALL dbcsr_copy(vectors%input_vec, rvec) 460 CALL dbcsr_matrix_colvec_multiply(pmat, vectors%input_vec, vectors%result_vec, 1.0_dp, & 461 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec) 462 CALL dbcsr_copy(zvec, vectors%result_vec) 463 rsnew = vec_dot_vec(rvec, zvec, mpgrp) 464 beta = rsnew/rsold 465 CALL dbcsr_add(pvec, zvec, beta, 1.0_dp) 466 rsold = rsnew 467 END DO 468 IF (SQRT(rsnew) < control%threshold) control%converged = .TRUE. 469 CALL dbcsr_release(zvec) 470 CALL dbcsr_release(pvec) 471 CALL dbcsr_release(rvec) 472 CALL dbcsr_release(apvec) 473 474 ELSE 475 CALL dbcsr_copy(pvec, xvec) 476 CALL dbcsr_copy(rvec, xvec) 477 CALL dbcsr_set(xvec, 0.0_dp) 478 ! Conjugate gradients 479 rsold = vec_dot_vec(rvec, rvec, mpgrp) 480 DO iter = 1, control%max_iter 481 CALL dbcsr_copy(vectors%input_vec, pvec) 482 CALL dbcsr_matrix_colvec_multiply(amat, vectors%input_vec, vectors%result_vec, 1.0_dp, & 483 0.0_dp, vectors%rep_row_vec, vectors%rep_col_vec) 484 CALL dbcsr_copy(apvec, vectors%result_vec) 485 pap = vec_dot_vec(pvec, apvec, mpgrp) 486 IF (ABS(pap) < 1.e-24_dp) THEN 487 alpha = 0.0_dp 488 ELSE 489 alpha = rsold/pap 490 END IF 491 CALL dbcsr_add(xvec, pvec, 1.0_dp, alpha) 492 CALL dbcsr_add(rvec, apvec, 1.0_dp, -alpha) 493 rsnew = vec_dot_vec(rvec, rvec, mpgrp) 494 IF (SQRT(rsnew) < control%threshold) EXIT 495 CPASSERT(alpha /= 0.0_dp) 496 beta = rsnew/rsold 497 CALL dbcsr_add(pvec, rvec, beta, 1.0_dp) 498 rsold = rsnew 499 END DO 500 IF (SQRT(rsnew) < control%threshold) control%converged = .TRUE. 501 CALL dbcsr_release(pvec) 502 CALL dbcsr_release(rvec) 503 CALL dbcsr_release(apvec) 504 END IF 505 506 END SUBROUTINE conjugate_gradient 507 508! ************************************************************************************************** 509!> \brief ... 510!> \param avec ... 511!> \param bvec ... 512!> \param mpgrp ... 513!> \return ... 514! ************************************************************************************************** 515 FUNCTION vec_dot_vec(avec, bvec, mpgrp) RESULT(adotb) 516 TYPE(dbcsr_type) :: avec, bvec 517 INTEGER, INTENT(IN) :: mpgrp 518 REAL(KIND=dp) :: adotb 519 520 INTEGER :: i, j 521 LOGICAL :: found 522 REAL(KIND=dp), DIMENSION(:, :), POINTER :: av, bv 523 TYPE(dbcsr_iterator_type) :: dbcsr_iter 524 525 adotb = 0.0_dp 526 CALL dbcsr_iterator_start(dbcsr_iter, avec) 527 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 528 CALL dbcsr_iterator_next_block(dbcsr_iter, i, j, av) 529 CALL dbcsr_get_block_p(bvec, i, j, bv, found) 530 IF (found .AND. SIZE(bv) > 0) THEN 531 adotb = adotb + DOT_PRODUCT(av(:, 1), bv(:, 1)) 532 END IF 533 END DO 534 CALL dbcsr_iterator_stop(dbcsr_iter) 535 CALL mp_sum(adotb, mpgrp) 536 537 END FUNCTION vec_dot_vec 538 539END MODULE arnoldi_api 540