1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief methods for arnoldi iteration 8!> \par History 9!> 2014.09 created [Florian Schiffmann] 10!> \author Florian Schiffmann 11! ************************************************************************************************** 12 13#:include 'arnoldi.fypp' 14MODULE arnoldi_methods 15 USE arnoldi_geev, ONLY: arnoldi_general_local_diag,& 16 arnoldi_symm_local_diag,& 17 arnoldi_tridiag_local_diag 18 USE arnoldi_types, ONLY: & 19 arnoldi_control_type, arnoldi_data_c_type, arnoldi_data_d_type, arnoldi_data_s_type, & 20 arnoldi_data_type, arnoldi_data_z_type, get_control, get_data_c, get_data_d, get_data_s, & 21 get_data_z, has_d_cmplx, has_d_real, has_s_cmplx, has_s_real, m_x_v_vectors_type 22 USE dbcsr_api, ONLY: & 23 dbcsr_add, dbcsr_copy, dbcsr_get_data_p, dbcsr_get_info, dbcsr_iterator_blocks_left, & 24 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 25 dbcsr_p_type, dbcsr_scale, dbcsr_type 26 USE dbcsr_vector, ONLY: dbcsr_matrix_colvec_multiply 27 USE kinds, ONLY: real_4,& 28 real_8 29 USE message_passing, ONLY: mp_bcast,& 30 mp_sum 31#include "../base/base_uses.f90" 32 33 IMPLICIT NONE 34 35 PRIVATE 36 37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_methods' 38 39 PUBLIC :: arnoldi_init, build_subspace, compute_evals, arnoldi_iram, & 40 gev_arnoldi_init, gev_build_subspace, gev_update_data 41 42 INTERFACE convert_matrix 43 MODULE PROCEDURE convert_matrix_z_to_d, convert_matrix_s_to_c 44 MODULE PROCEDURE convert_matrix_d_to_z, convert_matrix_c_to_s 45 MODULE PROCEDURE convert_matrix_z_to_z, convert_matrix_c_to_c 46 END INTERFACE 47 48CONTAINS 49 50! ************************************************************************************************** 51!> \brief Interface for the routine calcualting the implicit restarts 52!> Currently all based on lapack 53!> \param arnoldi_data ... 54! ************************************************************************************************** 55 SUBROUTINE arnoldi_iram(arnoldi_data) 56 TYPE(arnoldi_data_type) :: arnoldi_data 57 58 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_iram', routineP = moduleN//':'//routineN 59 60 INTEGER :: handle 61 62 CALL timeset(routineN, handle) 63 64 IF (has_d_real(arnoldi_data)) CALL arnoldi_iram_d(arnoldi_data) 65 IF (has_s_real(arnoldi_data)) CALL arnoldi_iram_s(arnoldi_data) 66 IF (has_d_cmplx(arnoldi_data)) CALL arnoldi_iram_z(arnoldi_data) 67 IF (has_s_cmplx(arnoldi_data)) CALL arnoldi_iram_c(arnoldi_data) 68 69 CALL timestop(handle) 70 71 END SUBROUTINE arnoldi_iram 72 73! ************************************************************************************************** 74!> \brief Interface to compute the eigenvalues of a nonsymmetric matrix 75!> This is only the serial version 76!> \param arnoldi_data ... 77! ************************************************************************************************** 78 SUBROUTINE compute_evals(arnoldi_data) 79 TYPE(arnoldi_data_type) :: arnoldi_data 80 81 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_evals', routineP = moduleN//':'//routineN 82 83 INTEGER :: handle 84 85 CALL timeset(routineN, handle) 86 87 IF (has_d_real(arnoldi_data)) CALL compute_evals_d(arnoldi_data) 88 IF (has_s_real(arnoldi_data)) CALL compute_evals_s(arnoldi_data) 89 IF (has_d_cmplx(arnoldi_data)) CALL compute_evals_z(arnoldi_data) 90 IF (has_s_cmplx(arnoldi_data)) CALL compute_evals_c(arnoldi_data) 91 92 CALL timestop(handle) 93 94 END SUBROUTINE compute_evals 95 96! ************************************************************************************************** 97!> \brief Interface for the initialization of the arnoldi subspace creation 98!> currently it can only setup a random vector but can be improved to 99!> various types of restarts easily 100!> \param matrix pointer to the matrices as described in main interface 101!> \param vectors work vectors for the matrix vector multiplications 102!> \param arnoldi_data all data concerning the subspace 103! ************************************************************************************************** 104 SUBROUTINE arnoldi_init(matrix, vectors, arnoldi_data) 105 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 106 TYPE(m_x_v_vectors_type) :: vectors 107 TYPE(arnoldi_data_type) :: arnoldi_data 108 109 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_init', routineP = moduleN//':'//routineN 110 111 INTEGER :: handle 112 113 CALL timeset(routineN, handle) 114 115 IF (has_d_real(arnoldi_data)) CALL arnoldi_init_d(matrix, vectors, arnoldi_data) 116 IF (has_s_real(arnoldi_data)) CALL arnoldi_init_s(matrix, vectors, arnoldi_data) 117 IF (has_d_cmplx(arnoldi_data)) CALL arnoldi_init_z(matrix, vectors, arnoldi_data) 118 IF (has_s_cmplx(arnoldi_data)) CALL arnoldi_init_c(matrix, vectors, arnoldi_data) 119 120 CALL timestop(handle) 121 122 END SUBROUTINE arnoldi_init 123 124! ************************************************************************************************** 125!> \brief Interface for the initialization of the arnoldi subspace creation 126!> for the generalized eigenvalue problem 127!> \param matrix pointer to the matrices as described in main interface 128!> \param matrix_arnoldi ... 129!> \param vectors work vectors for the matrix vector multiplications 130!> \param arnoldi_data all data concerning the subspace 131! ************************************************************************************************** 132 SUBROUTINE gev_arnoldi_init(matrix, matrix_arnoldi, vectors, arnoldi_data) 133 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix, matrix_arnoldi 134 TYPE(m_x_v_vectors_type) :: vectors 135 TYPE(arnoldi_data_type) :: arnoldi_data 136 137 CHARACTER(LEN=*), PARAMETER :: routineN = 'gev_arnoldi_init', & 138 routineP = moduleN//':'//routineN 139 140 INTEGER :: handle 141 142 CALL timeset(routineN, handle) 143 144 IF (has_d_real(arnoldi_data)) CALL gev_arnoldi_init_d(matrix, matrix_arnoldi, vectors, arnoldi_data) 145 IF (has_s_real(arnoldi_data)) CALL gev_arnoldi_init_s(matrix, matrix_arnoldi, vectors, arnoldi_data) 146 IF (has_d_cmplx(arnoldi_data)) CALL gev_arnoldi_init_z(matrix, matrix_arnoldi, vectors, arnoldi_data) 147 IF (has_s_cmplx(arnoldi_data)) CALL gev_arnoldi_init_c(matrix, matrix_arnoldi, vectors, arnoldi_data) 148 149 CALL timestop(handle) 150 151 END SUBROUTINE gev_arnoldi_init 152 153! ************************************************************************************************** 154!> \brief here the iterations are performed and the krylov space is constructed 155!> \param matrix see above 156!> \param vectors see above 157!> \param arnoldi_data see above 158! ************************************************************************************************** 159 SUBROUTINE build_subspace(matrix, vectors, arnoldi_data) 160 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 161 TYPE(m_x_v_vectors_type) :: vectors 162 TYPE(arnoldi_data_type) :: arnoldi_data 163 164 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_subspace', routineP = moduleN//':'//routineN 165 166 INTEGER :: handle 167 168 CALL timeset(routineN, handle) 169 170 IF (has_d_real(arnoldi_data)) CALL build_subspace_d(matrix, vectors, arnoldi_data) 171 IF (has_s_real(arnoldi_data)) CALL build_subspace_s(matrix, vectors, arnoldi_data) 172 IF (has_d_cmplx(arnoldi_data)) CALL build_subspace_z(matrix, vectors, arnoldi_data) 173 IF (has_s_cmplx(arnoldi_data)) CALL build_subspace_c(matrix, vectors, arnoldi_data) 174 175 CALL timestop(handle) 176 177 END SUBROUTINE build_subspace 178 179! ************************************************************************************************** 180!> \brief here the iterations are performed and the krylov space for the generalized 181!> eigenvalue probelm is created 182!> \param matrix see above 183!> \param vectors see above 184!> \param arnoldi_data see above 185! ************************************************************************************************** 186 SUBROUTINE gev_build_subspace(matrix, vectors, arnoldi_data) 187 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 188 TYPE(m_x_v_vectors_type) :: vectors 189 TYPE(arnoldi_data_type) :: arnoldi_data 190 191 CHARACTER(LEN=*), PARAMETER :: routineN = 'gev_build_subspace', & 192 routineP = moduleN//':'//routineN 193 194 INTEGER :: handle 195 196 CALL timeset(routineN, handle) 197 198 IF (has_d_real(arnoldi_data)) CALL gev_build_subspace_d(matrix, vectors, arnoldi_data) 199 IF (has_s_real(arnoldi_data)) CALL gev_build_subspace_s(matrix, vectors, arnoldi_data) 200 IF (has_d_cmplx(arnoldi_data)) CALL gev_build_subspace_z(matrix, vectors, arnoldi_data) 201 IF (has_s_cmplx(arnoldi_data)) CALL gev_build_subspace_c(matrix, vectors, arnoldi_data) 202 203 CALL timestop(handle) 204 205 END SUBROUTINE gev_build_subspace 206 207! ************************************************************************************************** 208!> \brief in the generalized eigenvalue the matrix depende on the projection 209!> therefore the outer loop has to build a new set of matrices for the 210!> inner loop 211!> \param matrix see above 212!> \param matrix_arnoldi ... 213!> \param vectors ... 214!> \param arnoldi_data see above 215! ************************************************************************************************** 216 SUBROUTINE gev_update_data(matrix, matrix_arnoldi, vectors, arnoldi_data) 217 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix, matrix_arnoldi 218 TYPE(m_x_v_vectors_type) :: vectors 219 TYPE(arnoldi_data_type) :: arnoldi_data 220 221 CHARACTER(LEN=*), PARAMETER :: routineN = 'gev_update_data', & 222 routineP = moduleN//':'//routineN 223 224 INTEGER :: handle 225 226 CALL timeset(routineN, handle) 227 228 IF (has_d_real(arnoldi_data)) CALL gev_update_data_d(matrix, matrix_arnoldi, vectors, arnoldi_data) 229 IF (has_s_real(arnoldi_data)) CALL gev_update_data_s(matrix, matrix_arnoldi, vectors, arnoldi_data) 230 IF (has_d_cmplx(arnoldi_data)) CALL gev_update_data_z(matrix, matrix_arnoldi, vectors, arnoldi_data) 231 IF (has_s_cmplx(arnoldi_data)) CALL gev_update_data_c(matrix, matrix_arnoldi, vectors, arnoldi_data) 232 233 CALL timestop(handle) 234 235 END SUBROUTINE gev_update_data 236 237! ************************************************************************************************** 238!> \brief ... 239!> \param m_out ... 240!> \param m_in ... 241! ************************************************************************************************** 242 SUBROUTINE convert_matrix_z_to_d(m_out, m_in) 243 REAL(real_8), DIMENSION(:, :) :: m_out 244 COMPLEX(real_8), DIMENSION(:, :) :: m_in 245 246 m_out(:, :) = REAL(m_in(:, :), KIND=real_8) 247 END SUBROUTINE convert_matrix_z_to_d 248 249! ************************************************************************************************** 250!> \brief ... 251!> \param m_out ... 252!> \param m_in ... 253! ************************************************************************************************** 254 SUBROUTINE convert_matrix_c_to_s(m_out, m_in) 255 REAL(real_4), DIMENSION(:, :) :: m_out 256 COMPLEX(real_4), DIMENSION(:, :) :: m_in 257 258 m_out(:, :) = REAL(m_in, KIND=real_4) 259 END SUBROUTINE convert_matrix_c_to_s 260 261! ************************************************************************************************** 262!> \brief ... 263!> \param m_out ... 264!> \param m_in ... 265! ************************************************************************************************** 266 SUBROUTINE convert_matrix_d_to_z(m_out, m_in) 267 COMPLEX(real_8), DIMENSION(:, :) :: m_out 268 REAL(real_8), DIMENSION(:, :) :: m_in 269 270 m_out(:, :) = CMPLX(m_in, 0.0, KIND=real_8) 271 END SUBROUTINE convert_matrix_d_to_z 272 273! ************************************************************************************************** 274!> \brief ... 275!> \param m_out ... 276!> \param m_in ... 277! ************************************************************************************************** 278 SUBROUTINE convert_matrix_s_to_c(m_out, m_in) 279 COMPLEX(real_4), DIMENSION(:, :) :: m_out 280 REAL(real_4), DIMENSION(:, :) :: m_in 281 282 m_out(:, :) = CMPLX(m_in, 0.0, KIND=real_4) 283 END SUBROUTINE convert_matrix_s_to_c 284 285! I kno that one is stupid but like this it simplifies the template 286! ************************************************************************************************** 287!> \brief ... 288!> \param m_out ... 289!> \param m_in ... 290! ************************************************************************************************** 291 SUBROUTINE convert_matrix_z_to_z(m_out, m_in) 292 COMPLEX(real_8), DIMENSION(:, :) :: m_out, m_in 293 294 m_out(:, :) = m_in 295 END SUBROUTINE convert_matrix_z_to_z 296 297! ************************************************************************************************** 298!> \brief ... 299!> \param m_out ... 300!> \param m_in ... 301! ************************************************************************************************** 302 SUBROUTINE convert_matrix_c_to_c(m_out, m_in) 303 COMPLEX(real_4), DIMENSION(:, :) :: m_out, m_in 304 305 m_out(:, :) = m_in 306 END SUBROUTINE convert_matrix_c_to_c 307 308#:for nametype1, type_prec, type_nametype1, nametype_zero, nametype_one, nametype_negone, czero, cone, ctype, rnorm_to_norm, val_to_type in inst_params_2 309! ************************************************************************************************** 310!> \brief Call the correct eigensolver, in the arnoldi method only the right 311!> eigenvectors are used. Lefts are created here but dumped immediatly 312!> \param arnoldi_data ... 313! ************************************************************************************************** 314 SUBROUTINE compute_evals_${nametype1}$(arnoldi_data) 315 TYPE(arnoldi_data_type) :: arnoldi_data 316 317 CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_evals_${nametype1}$', & 318 routineP = moduleN//':'//routineN 319 320 COMPLEX(${type_prec}$), DIMENSION(:, :), ALLOCATABLE :: levec 321 TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data 322 INTEGER :: ndim 323 TYPE(arnoldi_control_type), POINTER :: control 324 325 ar_data=>get_data_${nametype1}$(arnoldi_data) 326 control=> get_control(arnoldi_data) 327 ndim=control%current_step 328 ALLOCATE(levec(ndim, ndim)) 329 330! Needs antoher interface as the calls to real and complex geev differ (sucks!) 331! only perform the diagonalization on processors which hold data 332 IF(control%generalized_ev)THEN 333 CALL arnoldi_symm_local_diag('V',ar_data%Hessenberg(1:ndim, 1:ndim), ndim,& 334 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim)) 335 ELSE 336 IF(control%symmetric)THEN 337 CALL arnoldi_tridiag_local_diag('N', 'V', ar_data%Hessenberg(1:ndim, 1:ndim), ndim, & 338 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim), levec) 339 ELSE 340 CALL arnoldi_general_local_diag('N', 'V', ar_data%Hessenberg(1:ndim, 1:ndim), ndim, & 341 ar_data%evals(1:ndim), ar_data%revec(1:ndim, 1:ndim), levec) 342 END IF 343 END IF 344 345 DEALLOCATE(levec) 346 347 END SUBROUTINE compute_evals_${nametype1}$ 348 349! ************************************************************************************************** 350!> \brief Initialization of the arnoldi vector. Here a random vector is used, 351!> might be generalized in the future 352!> \param matrix ... 353!> \param vectors ... 354!> \param arnoldi_data ... 355! ************************************************************************************************** 356 SUBROUTINE arnoldi_init_${nametype1}$ (matrix, vectors, arnoldi_data) 357 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 358 TYPE(m_x_v_vectors_type) :: vectors 359 TYPE(arnoldi_data_type) :: arnoldi_data 360 361 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_init_${nametype1}$', & 362 routineP = moduleN//':'//routineN 363 364 INTEGER :: i, iseed(4), row_size, col_size, & 365 nrow_local, ncol_local, pcol_group, & 366 row, col 367 REAL(${type_prec}$) :: rnorm 368 TYPE(dbcsr_iterator_type) :: iter 369 ${type_nametype1}$ :: norm 370 ${type_nametype1}$, DIMENSION(:), ALLOCATABLE :: & 371 v_vec, w_vec 372 ${type_nametype1}$, DIMENSION(:), POINTER :: data_vec 373 LOGICAL :: local_comp 374 TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data 375 TYPE(arnoldi_control_type), POINTER :: control 376 377 control=>get_control(arnoldi_data) 378 pcol_group=control%pcol_group 379 local_comp=control%local_comp 380 381 ar_data=>get_data_${nametype1}$(arnoldi_data) 382 383 ! create a local data copy to store the vectors and make Gram Schmidt a bit simpler 384 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local) 385 ALLOCATE(v_vec(nrow_local)) 386 ALLOCATE(w_vec(nrow_local)) 387 v_vec = ${nametype_zero}$ ; w_vec = ${nametype_zero}$ 388 ar_data%Hessenberg=${nametype_zero}$ 389 390 IF(control%has_initial_vector)THEN 391 ! after calling the set routine the initial vector is stored in f_vec 392 CALL transfer_local_array_to_dbcsr_${nametype1}$(vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp) 393 ELSE 394 ! Setup the initial normalized random vector (sufficient if it only happens on proc_col 0) 395 CALL dbcsr_iterator_start(iter, vectors%input_vec) 396 DO WHILE (dbcsr_iterator_blocks_left(iter)) 397 CALL dbcsr_iterator_next_block(iter, row, col, data_vec, row_size=row_size, col_size=col_size) 398 iseed(1)=2; iseed(2)=MOD(row, 4095); iseed(3)=MOD(col, 4095); iseed(4)=11 399 CALL ${nametype1}$larnv(2, iseed, row_size*col_size, data_vec) 400 END DO 401 CALL dbcsr_iterator_stop(iter) 402 END IF 403 404 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%input_vec, v_vec, nrow_local, control%local_comp) 405 406 ! compute the vector norm of the random vectorm, get it real valued as well (rnorm) 407 CALL compute_norms_${nametype1}$(v_vec, norm, rnorm, control%pcol_group) 408 409 IF (rnorm==0) rnorm=1 ! catch case where this rank has no actual data 410 CALL dbcsr_scale(vectors%input_vec, REAL(1.0, ${type_prec}$)/rnorm) 411 412 ! Everything prepared, initialize the Arnoldi iteration 413 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%input_vec, v_vec, nrow_local, control%local_comp) 414 415 ! This permits to compute the subspace of a matrix which is a product of multiple matrices 416 DO i=1, SIZE(matrix) 417 CALL dbcsr_matrix_colvec_multiply(matrix(i)%matrix, vectors%input_vec, vectors%result_vec, ${nametype_one}$, & 418 ${nametype_zero}$, vectors%rep_row_vec, vectors%rep_col_vec) 419 CALL dbcsr_copy(vectors%input_vec, vectors%result_vec) 420 END DO 421 422 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%result_vec, w_vec, nrow_local, control%local_comp) 423 424 ! Put the projection into the Hessenberg matrix, and make the vectors orthonormal 425 ar_data%Hessenberg(1, 1)=DOT_PRODUCT(v_vec, w_vec) 426 CALL mp_sum(ar_data%Hessenberg(1, 1), pcol_group) 427 ar_data%f_vec=w_vec-v_vec*ar_data%Hessenberg(1, 1) 428 429 ar_data%local_history(:, 1)=v_vec(:) 430 431 ! We did the first step in here so we should set the current step for the subspace generation accordingly 432 control%current_step=1 433 434 DEALLOCATE(v_vec, w_vec) 435 436 END SUBROUTINE arnoldi_init_${nametype1}$ 437 438! ************************************************************************************************** 439!> \brief Alogorithm for the implicit restarts in the arnoldi method 440!> this is an early implementaion which scales subspace size^4 441!> by replacing the lapack calls with direct math the 442!> QR and gemms can be made linear and a N^2 sacling will be acchieved 443!> however this already sets the framework but should be used with care 444!> \param arnoldi_data ... 445! ************************************************************************************************** 446 SUBROUTINE arnoldi_iram_${nametype1}$(arnoldi_data) 447 TYPE(arnoldi_data_type) :: arnoldi_data 448 449 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_iram_${nametype1}$', & 450 routineP = moduleN//':'//routineN 451 452 TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data 453 TYPE(arnoldi_control_type), POINTER :: control 454 COMPLEX(${type_prec}$), DIMENSION(:,:), ALLOCATABLE :: tmp_mat, safe_mat, Q, tmp_mat1 455 COMPLEX(${type_prec}$), DIMENSION(:), ALLOCATABLE :: work, tau, work_measure 456 INTEGER :: msize, lwork, i, j, info, nwant 457 ${type_nametype1}$ :: beta, sigma 458 ${type_nametype1}$,DIMENSION(:,:),ALLOCATABLE :: Qdata 459 460 461 ! This is just a terribly inefficient implementation but I hope it is correct and might serve as a reference 462 ar_data=>get_data_${nametype1}$(arnoldi_data) 463 control=>get_control(arnoldi_data) 464 msize=control%current_step 465 nwant=control%nval_out 466 ALLOCATE(tmp_mat(msize,msize)); ALLOCATE(safe_mat(msize,msize)) 467 ALLOCATE(Q(msize,msize)); ALLOCATE(tmp_mat1(msize,msize)) 468 ALLOCATE(work_measure(1)) 469 ALLOCATE(tau(msize)); ALLOCATE(Qdata(msize,msize)) 470 !make Q identity 471 Q=${czero}$ 472 DO i=1,msize 473 Q(i,i)=${cone}$ 474 END DO 475 476 ! Looks a bit odd, but safe_mat will contain the result in the end, while tmpmat gets violated by lapack 477 CALL convert_matrix(tmp_mat,ar_data%Hessenberg(1:msize,1:msize)) 478 safe_mat(:,:)=tmp_mat(:,:) 479 480 DO i=1,msize 481 ! A bit a strange check but in the end we only want to shift the unwanted evals 482 IF(ANY(control%selected_ind==i))CYCLE 483 ! Here we shift the matrix by subtracting unwanted evals from the diagonal 484 DO j=1,msize 485 tmp_mat(j,j)=tmp_mat(j,j)-ar_data%evals(i) 486 END DO 487 ! Now we repair the damage by QR factorizing 488 lwork=-1 489 CALL ${ctype}$geqrf(msize,msize,tmp_mat,msize,tau,work_measure,lwork,info) 490 lwork=INT(work_measure(1)) 491 IF (ALLOCATED(work)) THEN 492 IF (SIZE(work).LT.lwork) THEN 493 DEALLOCATE(work) 494 ENDIF 495 ENDIF 496 IF (.NOT.ALLOCATED(work)) ALLOCATE(work(lwork)) 497 CALL ${ctype}$geqrf(msize,msize,tmp_mat,msize,tau,work,lwork,info) 498 ! Ask Lapack to reconstruct Q from its own way of storing data (tmpmat will contain Q) 499 CALL ${ctype}$ungqr(msize,msize,msize,tmp_mat,msize,tau,work,lwork,info) 500 ! update Q=Q*Q_current 501 tmp_mat1(:,:)=Q(:,:) 502 CALL ${ctype}$gemm('N','N',msize,msize,msize,${cone}$,tmp_mat1,msize,tmp_mat,msize,${czero}$,& 503 Q,msize) 504 ! Update H=(Q*)HQ 505 CALL ${ctype}$gemm('C','N',msize,msize,msize,${cone}$,tmp_mat,msize,safe_mat,msize,${czero}$,& 506 tmp_mat1,msize) 507 CALL ${ctype}$gemm('N','N',msize,msize,msize,${cone}$,tmp_mat1,msize,tmp_mat,msize,${czero}$,& 508 safe_mat,msize) 509 510 ! this one is crucial for numerics not to accumulate noise in the subdiagonals 511 DO j=1,msize 512 safe_mat(j+2:msize,j)=${czero}$ 513 END DO 514 tmp_mat(:,:)=safe_mat(:,:) 515 END DO 516 517 ! Now we can compute our restart quantities 518 ar_data%Hessenberg=${nametype_zero}$ 519 CALL convert_matrix(ar_data%Hessenberg(1:msize,1:msize),safe_mat) 520 CALL convert_matrix(Qdata,Q) 521 522 beta=ar_data%Hessenberg(nwant+1,nwant); sigma=Qdata(msize,nwant) 523 524 !update the residuum and the basis vectors 525 IF(control%local_comp)THEN 526 ar_data%f_vec=MATMUL(ar_data%local_history(:,1:msize),Qdata(1:msize,nwant+1))*beta+ar_data%f_vec(:)*sigma 527 ar_data%local_history(:,1:nwant)=MATMUL(ar_data%local_history(:,1:msize),Qdata(1:msize,1:nwant)) 528 END IF 529 ! Set the current step to nwant so the subspace build knows where to start 530 control%current_step=nwant 531 532 DEALLOCATE(tmp_mat,safe_mat,Q,Qdata,tmp_mat1,work,tau,work_measure) 533 534 END SUBROUTINE arnoldi_iram_${nametype1}$ 535 536! ************************************************************************************************** 537!> \brief Here we create the Krylov subspace and fill the Hessenberg matrix 538!> convergence check is only performed on subspace convergence 539!> Gram Schidt is used to orthonogonalize. 540!> If this is numericall not sufficient a Daniel, Gragg, Kaufman and Steward 541!> correction is performed 542!> \param matrix ... 543!> \param vectors ... 544!> \param arnoldi_data ... 545! ************************************************************************************************** 546 SUBROUTINE build_subspace_${nametype1}$(matrix, vectors, arnoldi_data) 547 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 548 TYPE(m_x_v_vectors_type), TARGET :: vectors 549 TYPE(arnoldi_data_type) :: arnoldi_data 550 551 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_subspace_${nametype1}$', & 552 routineP = moduleN//':'//routineN 553 554 INTEGER :: i, j, ncol_local, nrow_local 555 REAL(${type_prec}$) :: rnorm 556 TYPE(arnoldi_control_type), POINTER :: control 557 TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data 558 ${type_nametype1}$ :: norm 559 ${type_nametype1}$, ALLOCATABLE, DIMENSION(:) :: h_vec, s_vec, v_vec, w_vec 560 TYPE(dbcsr_type), POINTER :: input_vec, result_vec, swap_vec 561 562 ar_data=>get_data_${nametype1}$(arnoldi_data) 563 control=>get_control(arnoldi_data) 564 control%converged=.FALSE. 565 566 ! create the vectors required during the iterations 567 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local) 568 ALLOCATE(v_vec(nrow_local)); ALLOCATE(w_vec(nrow_local)) 569 v_vec = ${nametype_zero}$ ; w_vec = ${nametype_zero}$ 570 ALLOCATE(s_vec(control%max_iter)); ALLOCATE(h_vec(control%max_iter)) 571 572 DO j=control%current_step, control%max_iter-1 573 574 ! compute the vector norm of the residuum, get it real valued as well (rnorm) 575 CALL compute_norms_${nametype1}$(ar_data%f_vec, norm, rnorm, control%pcol_group) 576 577 ! check convergence and inform everybody about it, a bit annoying to talk to everybody because of that 578 IF(control%myproc==0)control%converged=rnorm.LT.REAL(control%threshold, ${type_prec}$) 579 CALL mp_bcast(control%converged, 0, control%mp_group) 580 IF(control%converged)EXIT 581 582 ! transfer normalized residdum to history and its norm to the Hessenberg matrix 583 IF (rnorm==0) rnorm=1 ! catch case where this rank has no actual data 584 v_vec(:)=ar_data%f_vec(:)/rnorm; ar_data%local_history(:, j+1)=v_vec(:); ar_data%Hessenberg(j+1, j)=norm 585 586 input_vec=>vectors%input_vec 587 result_vec=>vectors%result_vec 588 CALL transfer_local_array_to_dbcsr_${nametype1}$(input_vec, v_vec, nrow_local, control%local_comp) 589 590 ! This permits to compute the subspace of a matrix which is a product of two matrices 591 DO i=1, SIZE(matrix) 592 CALL dbcsr_matrix_colvec_multiply(matrix(i)%matrix, input_vec, result_vec, ${nametype_one}$, & 593 ${nametype_zero}$, vectors%rep_row_vec, vectors%rep_col_vec) 594 swap_vec=>input_vec 595 input_vec=>result_vec 596 result_vec=>swap_vec 597 END DO 598 599 CALL transfer_dbcsr_to_local_array_${nametype1}$(input_vec, w_vec, nrow_local, control%local_comp) 600 601 ! Let's do the orthonormalization, to get the new f_vec. First try the Gram Schmidt scheme 602 CALL Gram_Schmidt_ortho_${nametype1}$(h_vec, ar_data%f_vec, s_vec, w_vec, nrow_local, j+1, & 603 ar_data%local_history, ar_data%local_history, control%local_comp, control%pcol_group) 604 605 ! A bit more expensive but simply always top up with a DGKS correction, otherwise numerics 606 ! can become a problem later on, there is probably a good check whether it's necessary, but we don't perform it 607 CALL DGKS_ortho_${nametype1}$(h_vec, ar_data%f_vec, s_vec, nrow_local, j+1, ar_data%local_history, & 608 ar_data%local_history, control%local_comp, control%pcol_group) 609 ! Finally we can put the projections into our Hessenberg matrix 610 ar_data%Hessenberg(1:j+1, j+1)= h_vec(1:j+1) 611 control%current_step=j+1 612 END DO 613 614 ! compute the vector norm of the final residuum and put it in to Hessenberg 615 CALL compute_norms_${nametype1}$(ar_data%f_vec, norm, rnorm, control%pcol_group) 616 ar_data%Hessenberg(control%current_step+1, control%current_step)=norm 617 618 ! broadcast the Hessenberg matrix so we don't need to care later on 619 CALL mp_bcast(ar_data%Hessenberg, 0, control%mp_group) 620 621 DEALLOCATE(v_vec, w_vec, h_vec, s_vec) 622 623 END SUBROUTINE build_subspace_${nametype1}$ 624 625! ************************************************************************************************** 626!> \brief Helper routine to transfer the all data of a dbcsr matrix to a local array 627!> \param vec ... 628!> \param array ... 629!> \param n ... 630!> \param is_local ... 631! ************************************************************************************************** 632 SUBROUTINE transfer_dbcsr_to_local_array_${nametype1}$(vec, array, n, is_local) 633 TYPE(dbcsr_type) :: vec 634 ${type_nametype1}$, DIMENSION(:) :: array 635 INTEGER :: n 636 LOGICAL :: is_local 637 ${type_nametype1}$, DIMENSION(:), POINTER :: data_vec 638 639 data_vec => dbcsr_get_data_p (vec, select_data_type=${nametype_zero}$) 640 IF(is_local) array(1:n)=data_vec(1:n) 641 642 END SUBROUTINE transfer_dbcsr_to_local_array_${nametype1}$ 643 644! ************************************************************************************************** 645!> \brief The inverse routine transfering data back from an array to a dbcsr 646!> \param vec ... 647!> \param array ... 648!> \param n ... 649!> \param is_local ... 650! ************************************************************************************************** 651 SUBROUTINE transfer_local_array_to_dbcsr_${nametype1}$(vec, array, n, is_local) 652 TYPE(dbcsr_type) :: vec 653 ${type_nametype1}$, DIMENSION(:) :: array 654 INTEGER :: n 655 LOGICAL :: is_local 656 ${type_nametype1}$, DIMENSION(:), POINTER :: data_vec 657 658 data_vec => dbcsr_get_data_p (vec, select_data_type=${nametype_zero}$) 659 IF(is_local) data_vec(1:n)=array(1:n) 660 661 END SUBROUTINE transfer_local_array_to_dbcsr_${nametype1}$ 662 663! ************************************************************************************************** 664!> \brief Gram-Schmidt in matrix vector form 665!> \param h_vec ... 666!> \param f_vec ... 667!> \param s_vec ... 668!> \param w_vec ... 669!> \param nrow_local ... 670!> \param j ... 671!> \param local_history ... 672!> \param reorth_mat ... 673!> \param local_comp ... 674!> \param pcol_group ... 675! ************************************************************************************************** 676 SUBROUTINE Gram_Schmidt_ortho_${nametype1}$(h_vec, f_vec, s_vec, w_vec, nrow_local,& 677 j, local_history, reorth_mat, local_comp, pcol_group) 678 ${type_nametype1}$, DIMENSION(:) :: h_vec, s_vec, f_vec, w_vec 679 ${type_nametype1}$, DIMENSION(:, :) :: local_history, reorth_mat 680 INTEGER :: nrow_local, j, pcol_group 681 LOGICAL :: local_comp 682 683 CHARACTER(LEN=*), PARAMETER :: routineN = 'Gram_Schmidt_ortho_${nametype1}$', & 684 routineP = moduleN//':'//routineN 685 INTEGER :: handle 686 687 CALL timeset(routineN, handle) 688 689 ! Let's do the orthonormalization, first try the Gram Schmidt scheme 690 h_vec=${nametype_zero}$; f_vec=${nametype_zero}$; s_vec=${nametype_zero}$ 691 IF(local_comp)CALL ${nametype1}$gemv('T', nrow_local, j, ${nametype_one}$, local_history, & 692 nrow_local, w_vec, 1, ${nametype_zero}$, h_vec, 1) 693 CALL mp_sum(h_vec(1:j), pcol_group) 694 695 IF(local_comp)CALL ${nametype1}$gemv('N', nrow_local, j, ${nametype_one}$, reorth_mat, & 696 nrow_local, h_vec, 1, ${nametype_zero}$, f_vec, 1) 697 f_vec(:)=w_vec(:)-f_vec(:) 698 699 CALL timestop(handle) 700 701 END SUBROUTINE Gram_Schmidt_ortho_${nametype1}$ 702 703! ************************************************************************************************** 704!> \brief Compute the Daniel, Gragg, Kaufman and Steward correction 705!> \param h_vec ... 706!> \param f_vec ... 707!> \param s_vec ... 708!> \param nrow_local ... 709!> \param j ... 710!> \param local_history ... 711!> \param reorth_mat ... 712!> \param local_comp ... 713!> \param pcol_group ... 714! ************************************************************************************************** 715 SUBROUTINE DGKS_ortho_${nametype1}$(h_vec, f_vec, s_vec, nrow_local, j, & 716 local_history, reorth_mat, local_comp, pcol_group) 717 ${type_nametype1}$, DIMENSION(:) :: h_vec, s_vec, f_vec 718 ${type_nametype1}$, DIMENSION(:, :) :: local_history, reorth_mat 719 INTEGER :: nrow_local, j, pcol_group 720 721 CHARACTER(LEN=*), PARAMETER :: routineN = 'DGKS_ortho_${nametype1}$', & 722 routineP = moduleN//':'//routineN 723 724 LOGICAL :: local_comp 725 INTEGER :: handle 726 727 CALL timeset(routineN, handle) 728 729 IF(local_comp)CALL ${nametype1}$gemv('T', nrow_local, j, ${nametype_one}$, local_history, & 730 nrow_local, f_vec, 1, ${nametype_zero}$, s_vec, 1) 731 CALL mp_sum(s_vec(1:j), pcol_group) 732 IF(local_comp)CALL ${nametype1}$gemv('N', nrow_local, j, ${nametype_negone}$, reorth_mat, & 733 nrow_local, s_vec, 1, ${nametype_one}$, f_vec, 1) 734 h_vec(1:j)=h_vec(1:j)+s_vec(1:j) 735 736 CALL timestop(handle) 737 738 END SUBROUTINE DGKS_ortho_${nametype1}$ 739 740! ************************************************************************************************** 741!> \brief Compute the norm of a vector distributed along proc_col 742!> as local arrays. Always return the real part next to the complex rep. 743!> \param vec ... 744!> \param norm ... 745!> \param rnorm ... 746!> \param pcol_group ... 747! ************************************************************************************************** 748 SUBROUTINE compute_norms_${nametype1}$(vec, norm, rnorm, pcol_group) 749 ${type_nametype1}$, DIMENSION(:) :: vec 750 REAL(${type_prec}$) :: rnorm 751 ${type_nametype1}$ :: norm 752 INTEGER :: pcol_group 753 754 ! the norm is computed along the processor column 755 norm=DOT_PRODUCT(vec, vec) 756 CALL mp_sum(norm, pcol_group) 757 rnorm=SQRT(REAL(norm, ${type_prec}$)) 758 ${rnorm_to_norm}$ 759 760 END SUBROUTINE compute_norms_${nametype1}$ 761 762! ************************************************************************************************** 763!> \brief Computes the initial guess for the solution of the generalized eigenvalue 764!> using the arnoldi method 765!> \param matrix ... 766!> \param matrix_arnoldi ... 767!> \param vectors ... 768!> \param arnoldi_data ... 769! ************************************************************************************************** 770 771 SUBROUTINE gev_arnoldi_init_${nametype1}$ (matrix, matrix_arnoldi, vectors, arnoldi_data) 772 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 773 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_arnoldi 774 TYPE(m_x_v_vectors_type) :: vectors 775 TYPE(arnoldi_data_type) :: arnoldi_data 776 777 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_init_${nametype1}$', & 778 routineP = moduleN//':'//routineN 779 780 INTEGER :: iseed(4), row_size, col_size, & 781 nrow_local, ncol_local, pcol_group, & 782 row, col 783 REAL(${type_prec}$) :: rnorm 784 TYPE(dbcsr_iterator_type) :: iter 785 ${type_nametype1}$ :: norm, denom 786 ${type_nametype1}$, DIMENSION(:), ALLOCATABLE :: & 787 v_vec, w_vec 788 ${type_nametype1}$, DIMENSION(:), POINTER :: data_vec 789 LOGICAL :: local_comp 790 TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data 791 TYPE(arnoldi_control_type), POINTER :: control 792 793 control=>get_control(arnoldi_data) 794 pcol_group=control%pcol_group 795 local_comp=control%local_comp 796 797 ar_data=>get_data_${nametype1}$(arnoldi_data) 798 799 ! create a local data copy to store the vectors and make Gram Schmidt a bit simpler 800 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local) 801 ALLOCATE(v_vec(nrow_local)) 802 ALLOCATE(w_vec(nrow_local)) 803 v_vec = ${nametype_zero}$ ; w_vec = ${nametype_zero}$ 804 ar_data%Hessenberg=${nametype_zero}$ 805 806 IF(control%has_initial_vector)THEN 807 ! after calling the set routine the initial vector is stored in f_vec 808 CALL transfer_local_array_to_dbcsr_${nametype1}$(vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp) 809 ELSE 810 ! Setup the initial normalized random vector (sufficient if it only happens on proc_col 0) 811 CALL dbcsr_iterator_start(iter, vectors%input_vec) 812 DO WHILE (dbcsr_iterator_blocks_left(iter)) 813 CALL dbcsr_iterator_next_block(iter, row, col, data_vec, row_size=row_size, col_size=col_size) 814 iseed(1)=2; iseed(2)=MOD(row, 4095); iseed(3)=MOD(col, 4095); iseed(4)=11 815 CALL ${nametype1}$larnv(2, iseed, row_size*col_size, data_vec) 816 END DO 817 CALL dbcsr_iterator_stop(iter) 818 END IF 819 820 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%input_vec, v_vec, nrow_local, control%local_comp) 821 822 ! compute the vector norm of the reandom vectorm, get it real valued as well (rnorm) 823 CALL compute_norms_${nametype1}$(v_vec, norm, rnorm, control%pcol_group) 824 825 IF (rnorm==0) rnorm=1 ! catch case where this rank has no actual data 826 CALL dbcsr_scale(vectors%input_vec, REAL(1.0, ${type_prec}$)/rnorm) 827 828 CALL dbcsr_matrix_colvec_multiply(matrix(1)%matrix, vectors%input_vec, vectors%result_vec, ${nametype_one}$, & 829 ${nametype_zero}$, vectors%rep_row_vec, vectors%rep_col_vec) 830 831 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%result_vec, w_vec, nrow_local, control%local_comp) 832 833 ar_data%rho_scale=${nametype_zero}$ 834 ar_data%rho_scale=DOT_PRODUCT(v_vec,w_vec) 835 CALL mp_sum(ar_data%rho_scale, pcol_group) 836 837 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, ${nametype_one}$, & 838 ${nametype_zero}$, vectors%rep_row_vec, vectors%rep_col_vec) 839 840 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%result_vec, w_vec, nrow_local, control%local_comp) 841 842 denom=${nametype_zero}$ 843 denom=DOT_PRODUCT(v_vec,w_vec) 844 CALL mp_sum(denom, pcol_group) 845 IF(control%myproc==0) ar_data%rho_scale=ar_data%rho_scale/denom 846 CALL mp_bcast(ar_data%rho_scale,0,control%mp_group) 847 848 ! if the maximum ev is requested we need to optimize with -A-rho*B 849 CALL dbcsr_copy(matrix_arnoldi(1)%matrix,matrix(1)%matrix) 850 CALL dbcsr_add(matrix_arnoldi(1)%matrix, matrix(2)%matrix, ${nametype_one}$, -ar_data%rho_scale) 851 852 ar_data%x_vec=v_vec 853 854 END SUBROUTINE gev_arnoldi_init_${nametype1}$ 855 856! ************************************************************************************************** 857!> \brief builds the basis rothogonal wrt. teh metric. 858!> The structure looks similar to normal arnoldi but norms, vectors and 859!> matrix_vector products are very differently defined. Therefore it is 860!> cleaner to put it in a separate subroutine to avoid confusion 861!> \param matrix ... 862!> \param vectors ... 863!> \param arnoldi_data ... 864! ************************************************************************************************** 865 SUBROUTINE gev_build_subspace_${nametype1}$(matrix, vectors, arnoldi_data) 866 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 867 TYPE(m_x_v_vectors_type) :: vectors 868 TYPE(arnoldi_data_type) :: arnoldi_data 869 870 CHARACTER(LEN=*), PARAMETER :: routineN = 'build_subspace_${nametype1}$', & 871 routineP = moduleN//':'//routineN 872 873 INTEGER :: j, ncol_local, nrow_local, pcol_group 874 TYPE(arnoldi_control_type), POINTER :: control 875 TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data 876 ${type_nametype1}$ :: norm 877 ${type_nametype1}$, ALLOCATABLE, DIMENSION(:) :: h_vec, s_vec, v_vec, w_vec 878 ${type_nametype1}$, ALLOCATABLE, DIMENSION(:,:) :: Zmat, CZmat , BZmat 879 880 ar_data=>get_data_${nametype1}$(arnoldi_data) 881 control=>get_control(arnoldi_data) 882 control%converged=.FALSE. 883 pcol_group=control%pcol_group 884 885 ! create the vectors required during the iterations 886 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local) 887 ALLOCATE(v_vec(nrow_local)); ALLOCATE(w_vec(nrow_local)) 888 v_vec = ${nametype_zero}$ ; w_vec = ${nametype_zero}$ 889 ALLOCATE(s_vec(control%max_iter)); ALLOCATE(h_vec(control%max_iter)) 890 ALLOCATE(Zmat(nrow_local,control%max_iter)); ALLOCATE(CZmat(nrow_local,control%max_iter)) 891 ALLOCATE(BZmat(nrow_local,control%max_iter)) 892 893 CALL transfer_local_array_to_dbcsr_${nametype1}$(vectors%input_vec, ar_data%x_vec, nrow_local, control%local_comp) 894 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, ${nametype_one}$, & 895 ${nametype_zero}$, vectors%rep_row_vec, vectors%rep_col_vec) 896 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%result_vec, BZmat(:,1), nrow_local, control%local_comp) 897 898 norm=${nametype_zero}$ 899 norm=DOT_PRODUCT(ar_data%x_vec,BZmat(:,1)) 900 CALL mp_sum(norm, pcol_group) 901 IF(control%local_comp)THEN 902 Zmat(:,1)=ar_data%x_vec/SQRT(norm); BZmat(:,1)= BZmat(:,1)/SQRT(norm) 903 END IF 904 905 DO j=1,control%max_iter 906 control%current_step=j 907 CALL transfer_local_array_to_dbcsr_${nametype1}$(vectors%input_vec, Zmat(:,j), nrow_local, control%local_comp) 908 CALL dbcsr_matrix_colvec_multiply(matrix(1)%matrix, vectors%input_vec, vectors%result_vec, ${nametype_one}$, & 909 ${nametype_zero}$, vectors%rep_row_vec, vectors%rep_col_vec) 910 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%result_vec, CZmat(:,j), nrow_local, control%local_comp) 911 w_vec(:)=CZmat(:,j) 912 913 ! Let's do the orthonormalization, to get the new f_vec. First try the Gram Schmidt scheme 914 CALL Gram_Schmidt_ortho_${nametype1}$(h_vec, ar_data%f_vec, s_vec, w_vec, nrow_local, j, & 915 BZmat, Zmat, control%local_comp, control%pcol_group) 916 917 ! A bit more expensive but simpliy always top up with a DGKS correction, otherwise numerics 918 ! can becom a problem later on, there is probably a good check, but we don't perform it 919 CALL DGKS_ortho_${nametype1}$(h_vec, ar_data%f_vec, s_vec, nrow_local, j, BZmat, & 920 Zmat, control%local_comp, control%pcol_group) 921 922 CALL transfer_local_array_to_dbcsr_${nametype1}$(vectors%input_vec, ar_data%f_vec, nrow_local, control%local_comp) 923 CALL dbcsr_matrix_colvec_multiply(matrix(2)%matrix, vectors%input_vec, vectors%result_vec, ${nametype_one}$, & 924 ${nametype_zero}$, vectors%rep_row_vec, vectors%rep_col_vec) 925 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%result_vec, v_vec, nrow_local, control%local_comp) 926 norm=${nametype_zero}$ 927 norm=DOT_PRODUCT(ar_data%f_vec,v_vec) 928 CALL mp_sum(norm, pcol_group) 929 930 IF(control%myproc==0)control%converged=REAL(norm,${type_prec}$).LT.EPSILON(REAL(1.0,${type_prec}$)) 931 CALL mp_bcast(control%converged, 0, control%mp_group) 932 IF(control%converged)EXIT 933 IF(j==control%max_iter-1)EXIT 934 935 IF(control%local_comp)THEN 936 Zmat(:,j+1)=ar_data%f_vec/SQRT(norm); BZmat(:,j+1)= v_vec(:)/SQRT(norm) 937 END IF 938 END DO 939 940! getting a bit more complicated here as the final matrix is again a product which has to be computed with the 941! ditributed vectors, therefore a sum along the first proc_col is necessary. As we want that matrix everywhere, 942! we set it to zero before and compute the distributed product only on the first col and then sum over the full grid 943 ar_data%Hessenberg=${nametype_zero}$ 944 IF(control%local_comp)THEN 945 ar_data%Hessenberg(1:control%current_step,1:control%current_step)=& 946 MATMUL(TRANSPOSE(CZmat(:,1:control%current_step)),Zmat(:,1:control%current_step)) 947 END IF 948 CALL mp_sum(ar_data%Hessenberg,control%mp_group) 949 950 ar_data%local_history=Zmat 951 ! broadcast the Hessenberg matrix so we don't need to care later on 952 953 DEALLOCATE(v_vec); DEALLOCATE(w_vec); DEALLOCATE(s_vec); DEALLOCATE(h_vec); DEALLOCATE(CZmat); 954 DEALLOCATE(Zmat); DEALLOCATE(BZmat) 955 956 END SUBROUTINE gev_build_subspace_${nametype1}$ 957 958! ************************************************************************************************** 959!> \brief Updates all data after an inner loop of the generalized ev arnoldi. 960!> Updates rho and C=A-rho*B accordingly. 961!> As an update scheme is used for he ev, the output ev has to be replaced 962!> with the updated one. 963!> Furthermore a convergence check is performed. The mv product could 964!> be skiiped by making clever use of the precomputed data, However, 965!> it is most likely not worth the effort. 966!> \param matrix ... 967!> \param matrix_arnoldi ... 968!> \param vectors ... 969!> \param arnoldi_data ... 970! ************************************************************************************************** 971 SUBROUTINE gev_update_data_${nametype1}$(matrix, matrix_arnoldi, vectors, arnoldi_data) 972 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix 973 TYPE(dbcsr_p_type), DIMENSION(:) :: matrix_arnoldi 974 TYPE(m_x_v_vectors_type) :: vectors 975 TYPE(arnoldi_data_type) :: arnoldi_data 976 977 CHARACTER(LEN=*), PARAMETER :: routineN = 'gev_update_data_${nametype1}$', & 978 routineP = moduleN//':'//routineN 979 980 TYPE(arnoldi_control_type), POINTER :: control 981 TYPE(arnoldi_data_${nametype1}$_type), POINTER :: ar_data 982 INTEGER :: ncol_local, nrow_local, ind, i 983 ${type_nametype1}$, ALLOCATABLE, DIMENSION(:) :: v_vec 984 REAL(${type_prec}$) :: rnorm 985 ${type_nametype1}$ :: norm 986 COMPLEX(${type_prec}$) :: val 987 988 control=>get_control(arnoldi_data) 989 990 ar_data=>get_data_${nametype1}$(arnoldi_data) 991 992! compute the new shift, hack around the problem templating the conversion 993 val=ar_data%evals(control%selected_ind(1)) 994 ar_data%rho_scale=ar_data%rho_scale+${val_to_type}$ 995! compute the new eigenvector / initial guess for the next arnoldi loop 996 ar_data%x_vec=${nametype_zero}$ 997 DO i=1,control%current_step 998 val=ar_data%revec(i,control%selected_ind(1)) 999 ar_data%x_vec(:)=ar_data%x_vec(:)+ar_data%local_history(:,i)*${val_to_type}$ 1000 END DO 1001! ar_data%x_vec(:)=MATMUL(ar_data%local_history(:,1:control%current_step),& 1002! ar_data%revec(1:control%current_step,control%selected_ind(1))) 1003 1004! update the C-matrix (A-rho*B), if teh maximum value is requested we have to use -A-rho*B 1005 CALL dbcsr_copy(matrix_arnoldi(1)%matrix,matrix(1)%matrix) 1006 CALL dbcsr_add(matrix_arnoldi(1)%matrix, matrix(2)%matrix, ${nametype_one}$, -ar_data%rho_scale) 1007 1008! compute convergence and set the correct eigenvalue and eigenvector 1009 CALL dbcsr_get_info(matrix=vectors%input_vec, nfullrows_local=nrow_local, nfullcols_local=ncol_local) 1010 IF (ncol_local>0) THEN 1011 ALLOCATE(v_vec(nrow_local)) 1012 CALL compute_norms_${nametype1}$(ar_data%x_vec, norm, rnorm, control%pcol_group) 1013 v_vec(:)=ar_data%x_vec(:)/rnorm 1014 ENDIF 1015 CALL transfer_local_array_to_dbcsr_${nametype1}$(vectors%input_vec, v_vec, nrow_local, control%local_comp) 1016 CALL dbcsr_matrix_colvec_multiply(matrix_arnoldi(1)%matrix, vectors%input_vec, vectors%result_vec, ${nametype_one}$, & 1017 ${nametype_zero}$, vectors%rep_row_vec, vectors%rep_col_vec) 1018 CALL transfer_dbcsr_to_local_array_${nametype1}$(vectors%result_vec, v_vec, nrow_local, control%local_comp) 1019 IF (ncol_local>0) THEN 1020 CALL compute_norms_${nametype1}$(v_vec, norm, rnorm, control%pcol_group) 1021 ! check convergence 1022 control%converged=rnorm.LT.control%threshold 1023 DEALLOCATE(v_vec) 1024 ENDIF 1025 ! and broadcast the real eigenvalue 1026 CALL mp_bcast(control%converged,0,control%mp_group) 1027 ind=control%selected_ind(1) 1028 CALL mp_bcast(ar_data%rho_scale,0,control%mp_group) 1029 1030! Again the maximum value request is done on -A therefore the eigenvalue needs the opposite sign 1031 ar_data%evals(ind)=ar_data%rho_scale 1032 1033 1034 END SUBROUTINE gev_update_data_${nametype1}$ 1035#:endfor 1036 1037END MODULE arnoldi_methods 1038