1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief provides a unified interface to lapack geev routines 8!> \par History 9!> 2014.09 created [Florian Schiffmann] 10!> \author Florian Schiffmann 11! ************************************************************************************************** 12 13MODULE arnoldi_geev 14 USE kinds, ONLY: real_4,& 15 real_8 16#include "../base/base_uses.f90" 17 18 IMPLICIT NONE 19 20 PRIVATE 21 22 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'arnoldi_geev' 23 24 PUBLIC :: arnoldi_general_local_diag, arnoldi_tridiag_local_diag, arnoldi_symm_local_diag 25 26 INTERFACE arnoldi_general_local_diag 27 MODULE PROCEDURE arnoldi_sgeev, arnoldi_dgeev, arnoldi_zgeev, arnoldi_cgeev 28 END INTERFACE 29 30 ! currently only specialzed for real matrices 31 INTERFACE arnoldi_tridiag_local_diag 32 MODULE PROCEDURE arnoldi_sstev, arnoldi_dstev, arnoldi_zgeev, arnoldi_cgeev 33 END INTERFACE 34 35 ! currently only specialzed for real matrices 36 INTERFACE arnoldi_symm_local_diag 37 MODULE PROCEDURE arnoldi_dsyevd, arnoldi_ssyevd, arnoldi_cheevd, arnoldi_zheevd 38 END INTERFACE 39 40CONTAINS 41 42! ************************************************************************************************** 43!> \brief ... 44!> \param jobvr ... 45!> \param matrix ... 46!> \param ndim ... 47!> \param evals ... 48!> \param revec ... 49! ************************************************************************************************** 50 SUBROUTINE arnoldi_zheevd(jobvr, matrix, ndim, evals, revec) 51 CHARACTER(1) :: jobvr 52 COMPLEX(real_8), DIMENSION(:, :) :: matrix 53 INTEGER :: ndim 54 COMPLEX(real_8), DIMENSION(:) :: evals 55 COMPLEX(real_8), DIMENSION(:, :) :: revec 56 57 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_zheevd', routineP = moduleN//':'//routineN 58 59 INTEGER :: i, info, liwork, lrwork, lwork, & 60 iwork(3 + 5*ndim) 61 COMPLEX(real_8) :: work(2*ndim + ndim**2), & 62 tmp_array(ndim, ndim) 63 REAL(real_8) :: rwork(1 + 5*ndim + 2*ndim**2) 64 65 tmp_array(:, :) = matrix(:, :) 66 lwork = 2*ndim + ndim**2 67 lrwork = 1 + 5*ndim + 2*ndim**2 68 liwork = 3 + 5*ndim 69 70 CALL zheevd(jobvr, 'U', ndim, tmp_array, evals, ndim, work, lwork, rwork, lrwork, iwork, liwork, info) 71 72 DO i = 1, ndim 73 revec(:, i) = tmp_array(:, i) 74 END DO 75 76 END SUBROUTINE arnoldi_zheevd 77 78! ************************************************************************************************** 79!> \brief ... 80!> \param jobvr ... 81!> \param matrix ... 82!> \param ndim ... 83!> \param evals ... 84!> \param revec ... 85! ************************************************************************************************** 86 SUBROUTINE arnoldi_cheevd(jobvr, matrix, ndim, evals, revec) 87 CHARACTER(1) :: jobvr 88 COMPLEX(real_4), DIMENSION(:, :) :: matrix 89 INTEGER :: ndim 90 COMPLEX(real_4), DIMENSION(:) :: evals 91 COMPLEX(real_4), DIMENSION(:, :) :: revec 92 93 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_cheevd', routineP = moduleN//':'//routineN 94 95 INTEGER :: i, info, liwork, lrwork, lwork, & 96 iwork(3 + 5*ndim) 97 COMPLEX(real_4) :: work(2*ndim + ndim**2), & 98 tmp_array(ndim, ndim) 99 REAL(real_4) :: rwork(1 + 5*ndim + 2*ndim**2) 100 101 tmp_array(:, :) = matrix(:, :) 102 lwork = 2*ndim + ndim**2 103 lrwork = 1 + 5*ndim + 2*ndim**2 104 liwork = 3 + 5*ndim 105 106 CALL zheevd(jobvr, 'U', ndim, tmp_array, evals, ndim, work, lwork, rwork, lrwork, iwork, liwork, info) 107 108 DO i = 1, ndim 109 revec(:, i) = tmp_array(:, i) 110 END DO 111 112 END SUBROUTINE arnoldi_cheevd 113 114! ************************************************************************************************** 115!> \brief ... 116!> \param jobvr ... 117!> \param matrix ... 118!> \param ndim ... 119!> \param evals ... 120!> \param revec ... 121! ************************************************************************************************** 122 SUBROUTINE arnoldi_dsyevd(jobvr, matrix, ndim, evals, revec) 123 CHARACTER(1) :: jobvr 124 REAL(real_8), DIMENSION(:, :) :: matrix 125 INTEGER :: ndim 126 COMPLEX(real_8), DIMENSION(:) :: evals 127 COMPLEX(real_8), DIMENSION(:, :) :: revec 128 129 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_dsyevd', routineP = moduleN//':'//routineN 130 131 INTEGER :: i, info, liwork, lwork, iwork(3 + 5*ndim) 132 REAL(real_8) :: tmp_array(ndim, ndim), & 133 work(1 + 6*ndim + 2*ndim**2) 134 REAL(real_8), DIMENSION(ndim) :: eval 135 136 lwork = 1 + 6*ndim + 2*ndim**2 137 liwork = 3 + 5*ndim 138 139 tmp_array(:, :) = matrix(:, :) 140 CALL dsyevd(jobvr, "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info) 141 142 DO i = 1, ndim 143 revec(:, i) = CMPLX(tmp_array(:, i), REAL(0.0, real_8), real_8) 144 evals(i) = CMPLX(eval(i), 0.0, real_8) 145 END DO 146 147 END SUBROUTINE arnoldi_dsyevd 148 149! ************************************************************************************************** 150!> \brief ... 151!> \param jobvr ... 152!> \param matrix ... 153!> \param ndim ... 154!> \param evals ... 155!> \param revec ... 156! ************************************************************************************************** 157 SUBROUTINE arnoldi_ssyevd(jobvr, matrix, ndim, evals, revec) 158 CHARACTER(1) :: jobvr 159 REAL(real_4), DIMENSION(:, :) :: matrix 160 INTEGER :: ndim 161 COMPLEX(real_4), DIMENSION(:) :: evals 162 COMPLEX(real_4), DIMENSION(:, :) :: revec 163 164 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_ssyevd', routineP = moduleN//':'//routineN 165 166 INTEGER :: i, info, liwork, lwork, iwork(3 + 5*ndim) 167 REAL(real_4) :: tmp_array(ndim, ndim), & 168 work(1 + 6*ndim + 2*ndim**2) 169 REAL(real_4), DIMENSION(ndim) :: eval 170 171 MARK_USED(jobvr) !the argument has to be here for the template to work 172 lwork = 1 + 6*ndim + 2*ndim**2 173 liwork = 3 + 5*ndim 174 175 tmp_array(:, :) = matrix(:, :) 176 CALL ssyevd("V", "U", ndim, tmp_array, ndim, eval, work, lwork, iwork, liwork, info) 177 178 DO i = 1, ndim 179 revec(:, i) = CMPLX(tmp_array(:, i), REAL(0.0, real_4), real_4) 180 evals(i) = CMPLX(eval(i), 0.0, real_4) 181 END DO 182 183 END SUBROUTINE arnoldi_ssyevd 184 185! ************************************************************************************************** 186!> \brief ... 187!> \param jobvl ... 188!> \param jobvr ... 189!> \param matrix ... 190!> \param ndim ... 191!> \param evals ... 192!> \param revec ... 193!> \param levec ... 194! ************************************************************************************************** 195 SUBROUTINE arnoldi_sstev(jobvl, jobvr, matrix, ndim, evals, revec, levec) 196 CHARACTER(1) :: jobvl, jobvr 197 REAL(real_4), DIMENSION(:, :) :: matrix 198 INTEGER :: ndim 199 COMPLEX(real_4), DIMENSION(:) :: evals 200 COMPLEX(real_4), DIMENSION(:, :) :: revec, levec 201 202 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_sstev', routineP = moduleN//':'//routineN 203 204 INTEGER :: i, info 205 REAL(real_4) :: work(20*ndim) 206 REAL(real_4), DIMENSION(ndim) :: diag, offdiag 207 REAL(real_4), DIMENSION(ndim, ndim) :: evec_r 208 209 MARK_USED(jobvl) !the argument has to be here for the template to work 210 211 levec(1, 1) = CMPLX(0.0, 0.0, real_4) 212 info = 0 213 diag(ndim) = matrix(ndim, ndim) 214 DO i = 1, ndim - 1 215 diag(i) = matrix(i, i) 216 offdiag(i) = matrix(i + 1, i) 217 END DO 218 219 CALL sstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info) 220 221 DO i = 1, ndim 222 revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, real_4), real_4) 223 evals(i) = CMPLX(diag(i), 0.0, real_4) 224 END DO 225 END SUBROUTINE arnoldi_sstev 226 227! ************************************************************************************************** 228!> \brief ... 229!> \param jobvl ... 230!> \param jobvr ... 231!> \param matrix ... 232!> \param ndim ... 233!> \param evals ... 234!> \param revec ... 235!> \param levec ... 236! ************************************************************************************************** 237 SUBROUTINE arnoldi_dstev(jobvl, jobvr, matrix, ndim, evals, revec, levec) 238 CHARACTER(1) :: jobvl, jobvr 239 REAL(real_8), DIMENSION(:, :) :: matrix 240 INTEGER :: ndim 241 COMPLEX(real_8), DIMENSION(:) :: evals 242 COMPLEX(real_8), DIMENSION(:, :) :: revec, levec 243 244 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_dstev', routineP = moduleN//':'//routineN 245 246 INTEGER :: i, info 247 REAL(real_8) :: work(20*ndim) 248 REAL(real_8), DIMENSION(ndim) :: diag, offdiag 249 REAL(real_8), DIMENSION(ndim, ndim) :: evec_r 250 251 MARK_USED(jobvl) !the argument has to be here for the template to work 252 253 levec(1, 1) = CMPLX(0.0, 0.0, real_8) 254 info = 0 255 diag(ndim) = matrix(ndim, ndim) 256 DO i = 1, ndim - 1 257 diag(i) = matrix(i, i) 258 offdiag(i) = matrix(i + 1, i) 259 260 END DO 261 262 CALL dstev(jobvr, ndim, diag, offdiag, evec_r, ndim, work, info) 263 264 DO i = 1, ndim 265 revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, real_8), real_8) 266 evals(i) = CMPLX(diag(i), 0.0, real_8) 267 END DO 268 END SUBROUTINE arnoldi_dstev 269! ************************************************************************************************** 270!> \brief ... 271!> \param jobvl ... 272!> \param jobvr ... 273!> \param matrix ... 274!> \param ndim ... 275!> \param evals ... 276!> \param revec ... 277!> \param levec ... 278! ************************************************************************************************** 279 SUBROUTINE arnoldi_sgeev(jobvl, jobvr, matrix, ndim, evals, revec, levec) 280 CHARACTER(1) :: jobvl, jobvr 281 REAL(real_4), DIMENSION(:, :) :: matrix 282 INTEGER :: ndim 283 COMPLEX(real_4), DIMENSION(:) :: evals 284 COMPLEX(real_4), DIMENSION(:, :) :: revec, levec 285 286 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_sgeev', routineP = moduleN//':'//routineN 287 288 INTEGER :: i, info, lwork 289 LOGICAL :: selects(ndim) 290 REAL(real_4) :: norm, tmp_array(ndim, ndim), & 291 work(20*ndim) 292 REAL(real_4), DIMENSION(ndim) :: eval1, eval2 293 REAL(real_4), DIMENSION(ndim, ndim) :: evec_l, evec_r 294 295 MARK_USED(jobvr) !the argument has to be here for the template to work 296 MARK_USED(jobvl) !the argument has to be here for the template to work 297 298 eval1 = REAL(0.0, real_4); eval2 = REAL(0.0, real_4) 299 tmp_array(:, :) = matrix(:, :) 300 ! ask lapack how much space it would like in the work vector, don't ask me why 301 lwork = -1 302 CALL shseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info) 303 304 lwork = MIN(20*ndim, INT(work(1))) 305 CALL shseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info) 306 CALL strevc('R', 'B', selects, ndim, tmp_array, ndim, evec_l, ndim, evec_r, ndim, ndim, ndim, work, info) 307 308 ! compose the eigenvectors, lapacks way of storing them is a pain 309 ! if eval is complex, then the complex conj pair of evec can be constructed from the i and i+1st evec 310 ! Unfortunately dtrevc computes the ev such that the largest is set to one and not normalized 311 i = 1 312 DO WHILE (i .LE. ndim) 313 IF (ABS(eval2(i)) .LT. EPSILON(REAL(0.0, real_4))) THEN 314 evec_r(:, i) = evec_r(:, i)/SQRT(DOT_PRODUCT(evec_r(:, i), evec_r(:, i))) 315 revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, real_4), real_4) 316 levec(:, i) = CMPLX(evec_l(:, i), REAL(0.0, real_4), real_4) 317 i = i + 1 318 ELSE IF (eval2(i) .GT. EPSILON(REAL(0.0, real_4))) THEN 319 norm = SQRT(SUM(evec_r(:, i)**2.0_real_4) + SUM(evec_r(:, i + 1)**2.0_real_4)) 320 revec(:, i) = CMPLX(evec_r(:, i), evec_r(:, i + 1), real_4)/norm 321 revec(:, i + 1) = CMPLX(evec_r(:, i), -evec_r(:, i + 1), real_4)/norm 322 levec(:, i) = CMPLX(evec_l(:, i), evec_l(:, i + 1), real_4) 323 levec(:, i + 1) = CMPLX(evec_l(:, i), -evec_l(:, i + 1), real_4) 324 i = i + 2 325 ELSE 326 CPABORT('something went wrong while sorting the EV in arnoldi_geev') 327 END IF 328 END DO 329 330 ! this is to keep the interface consistent with complex geev 331 DO i = 1, ndim 332 evals(i) = CMPLX(eval1(i), eval2(i), real_4) 333 END DO 334 335 END SUBROUTINE arnoldi_sgeev 336 337! ************************************************************************************************** 338!> \brief ... 339!> \param jobvl ... 340!> \param jobvr ... 341!> \param matrix ... 342!> \param ndim ... 343!> \param evals ... 344!> \param revec ... 345!> \param levec ... 346! ************************************************************************************************** 347 SUBROUTINE arnoldi_dgeev(jobvl, jobvr, matrix, ndim, evals, revec, levec) 348 CHARACTER(1) :: jobvl, jobvr 349 REAL(real_8), DIMENSION(:, :) :: matrix 350 INTEGER :: ndim 351 COMPLEX(real_8), DIMENSION(:) :: evals 352 COMPLEX(real_8), DIMENSION(:, :) :: revec, levec 353 354 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_dgeev', routineP = moduleN//':'//routineN 355 356 INTEGER :: i, info, lwork 357 LOGICAL :: selects(ndim) 358 REAL(real_8) :: norm, tmp_array(ndim, ndim), & 359 work(20*ndim) 360 REAL(real_8), DIMENSION(ndim) :: eval1, eval2 361 REAL(real_8), DIMENSION(ndim, ndim) :: evec_l, evec_r 362 363 MARK_USED(jobvr) !the argument has to be here for the template to work 364 MARK_USED(jobvl) !the argument has to be here for the template to work 365 366 eval1 = REAL(0.0, real_8); eval2 = REAL(0.0, real_8) 367 tmp_array(:, :) = matrix(:, :) 368 ! ask lapack how much space it would like in the work vector, don't ask me why 369 lwork = -1 370 CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info) 371 372 lwork = MIN(20*ndim, INT(work(1))) 373 CALL dhseqr('S', 'I', ndim, 1, ndim, tmp_array, ndim, eval1, eval2, evec_r, ndim, work, lwork, info) 374 CALL dtrevc('R', 'B', selects, ndim, tmp_array, ndim, evec_l, ndim, evec_r, ndim, ndim, ndim, work, info) 375 376 ! compose the eigenvectors, lapacks way of storing them is a pain 377 ! if eval is complex, then the complex conj pair of evec can be constructed from the i and i+1st evec 378 ! Unfortunately dtrevc computes the ev such that the largest is set to one and not normalized 379 i = 1 380 DO WHILE (i .LE. ndim) 381 IF (ABS(eval2(i)) .LT. EPSILON(REAL(0.0, real_8))) THEN 382 evec_r(:, i) = evec_r(:, i)/SQRT(DOT_PRODUCT(evec_r(:, i), evec_r(:, i))) 383 revec(:, i) = CMPLX(evec_r(:, i), REAL(0.0, real_8), real_8) 384 levec(:, i) = CMPLX(evec_l(:, i), REAL(0.0, real_8), real_8) 385 i = i + 1 386 ELSE IF (eval2(i) .GT. EPSILON(REAL(0.0, real_8))) THEN 387 norm = SQRT(SUM(evec_r(:, i)**2.0_real_8) + SUM(evec_r(:, i + 1)**2.0_real_8)) 388 revec(:, i) = CMPLX(evec_r(:, i), evec_r(:, i + 1), real_8)/norm 389 revec(:, i + 1) = CMPLX(evec_r(:, i), -evec_r(:, i + 1), real_8)/norm 390 levec(:, i) = CMPLX(evec_l(:, i), evec_l(:, i + 1), real_8) 391 levec(:, i + 1) = CMPLX(evec_l(:, i), -evec_l(:, i + 1), real_8) 392 i = i + 2 393 ELSE 394 CPABORT('something went wrong while sorting the EV in arnoldi_geev') 395 END IF 396 END DO 397 398 ! this is to keep the interface consistent with complex geev 399 DO i = 1, ndim 400 evals(i) = CMPLX(eval1(i), eval2(i), real_8) 401 END DO 402 403 END SUBROUTINE arnoldi_dgeev 404 405! ************************************************************************************************** 406!> \brief ... 407!> \param jobvl ... 408!> \param jobvr ... 409!> \param matrix ... 410!> \param ndim ... 411!> \param evals ... 412!> \param revec ... 413!> \param levec ... 414! ************************************************************************************************** 415 SUBROUTINE arnoldi_zgeev(jobvl, jobvr, matrix, ndim, evals, revec, levec) 416 CHARACTER(1) :: jobvl, jobvr 417 COMPLEX(real_8), DIMENSION(:, :) :: matrix 418 INTEGER :: ndim 419 COMPLEX(real_8), DIMENSION(:) :: evals 420 COMPLEX(real_8), DIMENSION(:, :) :: revec, levec 421 422 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_zgeev', routineP = moduleN//':'//routineN 423 424 INTEGER :: info, lwork 425 COMPLEX(real_8) :: work(20*ndim), tmp_array(ndim, ndim) 426 REAL(real_8) :: work2(2*ndim) 427 428 evals = CMPLX(0.0, 0.0, real_8) 429 ! ask lapack how much space it would like in the work vector, don't ask me why 430 lwork = -1 431 CALL ZGEEV(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info) 432 lwork = MIN(20*ndim, INT(work(1))) 433 434 tmp_array(:, :) = matrix(:, :) 435 CALL ZGEEV(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info) 436 437 END SUBROUTINE arnoldi_zgeev 438 439! ************************************************************************************************** 440!> \brief ... 441!> \param jobvl ... 442!> \param jobvr ... 443!> \param matrix ... 444!> \param ndim ... 445!> \param evals ... 446!> \param revec ... 447!> \param levec ... 448! ************************************************************************************************** 449 SUBROUTINE arnoldi_cgeev(jobvl, jobvr, matrix, ndim, evals, revec, levec) 450 CHARACTER(1) :: jobvl, jobvr 451 COMPLEX(real_4), DIMENSION(:, :) :: matrix 452 INTEGER :: ndim 453 COMPLEX(real_4), DIMENSION(:) :: evals 454 COMPLEX(real_4), DIMENSION(:, :) :: revec, levec 455 456 CHARACTER(LEN=*), PARAMETER :: routineN = 'arnoldi_cgeev', routineP = moduleN//':'//routineN 457 458 INTEGER :: info, lwork 459 COMPLEX(real_4) :: work(20*ndim), tmp_array(ndim, ndim) 460 REAL(real_4) :: work2(2*ndim) 461 462 evals = CMPLX(0.0, 0.0, real_4) 463 ! ask lapack how much space it would like in the work vector, don't ask me why 464 lwork = -1 465 CALL CGEEV(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info) 466 lwork = MIN(20*ndim, INT(work(1))) 467 468 tmp_array(:, :) = matrix(:, :) 469 CALL CGEEV(jobvl, jobvr, ndim, tmp_array, ndim, evals, levec, ndim, revec, ndim, work, lwork, work2, info) 470 471 END SUBROUTINE arnoldi_cgeev 472 473END MODULE arnoldi_geev 474