Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Collection of simple mathematical functions and subroutines
10 : !> \par History
11 : !> FUNCTION angle updated and FUNCTION dihedral angle added; cleaned
12 : !> (13.03.2004,MK)
13 : !> \author MK (15.11.1998)
14 : ! **************************************************************************************************
15 : MODULE mathlib
16 :
17 : USE kinds, ONLY: default_string_length, &
18 : dp
19 : USE mathconstants, ONLY: euler, &
20 : fac, &
21 : oorootpi, &
22 : z_one, &
23 : z_zero
24 : #include "../base/base_uses.f90"
25 : #:include 'expint.fypp'
26 :
27 : IMPLICIT NONE
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mathlib'
31 : REAL(KIND=dp), PARAMETER :: eps_geo = 1.0E-6_dp
32 :
33 : ! Public subroutines
34 :
35 : PUBLIC :: build_rotmat, &
36 : jacobi, &
37 : diamat_all, &
38 : invmat, &
39 : invmat_symm, &
40 : invert_matrix, &
41 : get_pseudo_inverse_svd, &
42 : get_pseudo_inverse_diag, &
43 : symmetrize_matrix, &
44 : unit_matrix, diag, &
45 : erfc_cutoff, &
46 : diag_antisym, &
47 : diag_complex, &
48 : gemm_square, &
49 : geeig_right
50 :
51 : ! Public functions
52 :
53 : PUBLIC :: angle, &
54 : binomial, &
55 : binomial_gen, &
56 : multinomial, &
57 : det_3x3, &
58 : dihedral_angle, &
59 : digamma, &
60 : gcd, &
61 : inv_3x3, &
62 : lcm, &
63 : vector_product, &
64 : pswitch, &
65 : rotate_vector, &
66 : reflect_vector, &
67 : expint, abnormal_value, &
68 : get_diag, &
69 : set_diag
70 :
71 : INTERFACE det_3x3
72 : MODULE PROCEDURE det_3x3_1, det_3x3_2
73 : END INTERFACE
74 :
75 : INTERFACE invert_matrix
76 : MODULE PROCEDURE invert_matrix_d, invert_matrix_z
77 : END INTERFACE
78 :
79 : INTERFACE set_diag
80 : MODULE PROCEDURE set_diag_scalar_d, set_diag_scalar_z
81 : END INTERFACE
82 :
83 : INTERFACE swap
84 : MODULE PROCEDURE swap_scalar, swap_vector
85 : END INTERFACE
86 :
87 : INTERFACE unit_matrix
88 : MODULE PROCEDURE unit_matrix_d, unit_matrix_z
89 : END INTERFACE
90 :
91 : INTERFACE gemm_square
92 : MODULE PROCEDURE zgemm_square_2, zgemm_square_3, dgemm_square_2, dgemm_square_3
93 : END INTERFACE
94 :
95 : CONTAINS
96 :
97 : ! **************************************************************************************************
98 : !> \brief Polynomial (5th degree) switching function
99 : !> f(a) = 1 .... f(b) = 0 with f'(a) = f"(a) = f'(b) = f"(b) = 0
100 : !> \param x ...
101 : !> \param a ...
102 : !> \param b ...
103 : !> \param order ...
104 : !> \return =0 : f(x)
105 : !> \return =1 : f'(x)
106 : !> \return =2 : f"(x)
107 : ! **************************************************************************************************
108 29870 : FUNCTION pswitch(x, a, b, order) RESULT(fx)
109 : REAL(KIND=dp) :: x, a, b
110 : INTEGER :: order
111 : REAL(KIND=dp) :: fx
112 :
113 : REAL(KIND=dp) :: u, u2, u3
114 :
115 29870 : CPASSERT(b > a)
116 29870 : IF (x < a .OR. x > b) THEN
117 : ! outside switching intervall
118 27934 : IF (order > 0) THEN
119 : ! derivatives are 0
120 : fx = 0.0_dp
121 : ELSE
122 13967 : IF (x < a) THEN
123 : ! x < a => f(x) = 1
124 : fx = 1.0_dp
125 : ELSE
126 : ! x > b => f(x) = 0
127 13640 : fx = 0.0_dp
128 : END IF
129 : END IF
130 : ELSE
131 : ! renormalized coordinate
132 1936 : u = (x - a)/(b - a)
133 2904 : SELECT CASE (order)
134 : CASE (0)
135 968 : u2 = u*u
136 968 : u3 = u2*u
137 968 : fx = 1._dp - 10._dp*u3 + 15._dp*u2*u2 - 6._dp*u2*u3
138 : CASE (1)
139 968 : u2 = u*u
140 968 : fx = -30._dp*u2 + 60._dp*u*u2 - 30._dp*u2*u2
141 968 : fx = fx/(b - a)
142 : CASE (2)
143 0 : u2 = u*u
144 0 : fx = -60._dp*u + 180._dp*u2 - 120._dp*u*u2
145 0 : fx = fx/(b - a)**2
146 : CASE DEFAULT
147 1936 : CPABORT('order not defined')
148 : END SELECT
149 : END IF
150 :
151 29870 : END FUNCTION pswitch
152 :
153 : ! **************************************************************************************************
154 : !> \brief determines if a value is not normal (e.g. for Inf and Nan)
155 : !> based on IO to work also under optimization.
156 : !> \param a input value
157 : !> \return TRUE for NaN and Inf
158 : ! **************************************************************************************************
159 4035837 : LOGICAL FUNCTION abnormal_value(a)
160 : REAL(KIND=dp) :: a
161 :
162 : CHARACTER(LEN=32) :: buffer
163 :
164 4035837 : abnormal_value = .FALSE.
165 : ! the function should work when compiled with -ffast-math and similar
166 : ! unfortunately, that option asserts that all numbers are normals,
167 : ! which the compiler uses to optimize the function to .FALSE. if based on the IEEE module
168 : ! therefore, pass this to the Fortran runtime/printf, if things are NaN or Inf, error out.
169 4035837 : WRITE (buffer, *) a
170 4035837 : IF (INDEX(buffer, "N") /= 0 .OR. INDEX(buffer, "n") /= 0) abnormal_value = .TRUE.
171 :
172 4035837 : END FUNCTION abnormal_value
173 :
174 : ! **************************************************************************************************
175 : !> \brief Calculation of the angle between the vectors a and b.
176 : !> The angle is returned in radians.
177 : !> \param a ...
178 : !> \param b ...
179 : !> \return ...
180 : !> \date 14.10.1998
181 : !> \author MK
182 : !> \version 1.0
183 : ! **************************************************************************************************
184 751660 : PURE FUNCTION angle(a, b) RESULT(angle_ab)
185 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: a, b
186 : REAL(KIND=dp) :: angle_ab
187 :
188 : REAL(KIND=dp) :: length_of_a, length_of_b
189 751660 : REAL(KIND=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
190 :
191 3006640 : length_of_a = NORM2(a)
192 3006640 : length_of_b = NORM2(b)
193 :
194 751660 : IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
195 3006640 : a_norm(:) = a(:)/length_of_a
196 3006640 : b_norm(:) = b(:)/length_of_b
197 3006640 : angle_ab = ACOS(MIN(MAX(DOT_PRODUCT(a_norm, b_norm), -1.0_dp), 1.0_dp))
198 : ELSE
199 : angle_ab = 0.0_dp
200 : END IF
201 :
202 751660 : END FUNCTION angle
203 :
204 : ! **************************************************************************************************
205 : !> \brief The binomial coefficient n over k for 0 <= k <= n is calculated,
206 : !> otherwise zero is returned.
207 : !> \param n ...
208 : !> \param k ...
209 : !> \return ...
210 : !> \date 08.03.1999
211 : !> \author MK
212 : !> \version 1.0
213 : ! **************************************************************************************************
214 562419623 : ELEMENTAL FUNCTION binomial(n, k) RESULT(n_over_k)
215 : INTEGER, INTENT(IN) :: n, k
216 : REAL(KIND=dp) :: n_over_k
217 :
218 562419623 : IF ((k >= 0) .AND. (k <= n)) THEN
219 559041211 : n_over_k = fac(n)/(fac(n - k)*fac(k))
220 : ELSE
221 : n_over_k = 0.0_dp
222 : END IF
223 :
224 562419623 : END FUNCTION binomial
225 :
226 : ! **************************************************************************************************
227 : !> \brief The generalized binomial coefficient z over k for 0 <= k <= n is calculated.
228 : !> (z) z*(z-1)*...*(z-k+2)*(z-k+1)
229 : !> ( ) = ---------------------------
230 : !> (k) k!
231 : !> \param z ...
232 : !> \param k ...
233 : !> \return ...
234 : !> \date 11.11.2019
235 : !> \author FS
236 : !> \version 1.0
237 : ! **************************************************************************************************
238 21826063 : ELEMENTAL FUNCTION binomial_gen(z, k) RESULT(z_over_k)
239 : REAL(KIND=dp), INTENT(IN) :: z
240 : INTEGER, INTENT(IN) :: k
241 : REAL(KIND=dp) :: z_over_k
242 :
243 : INTEGER :: i
244 :
245 21826063 : IF (k >= 0) THEN
246 : z_over_k = 1.0_dp
247 35391981 : DO i = 1, k
248 35391981 : z_over_k = z_over_k*(z - i + 1)/REAL(i, dp)
249 : END DO
250 : ELSE
251 : z_over_k = 0.0_dp
252 : END IF
253 :
254 21826063 : END FUNCTION binomial_gen
255 :
256 : ! **************************************************************************************************
257 : !> \brief Calculates the multinomial coefficients
258 : !> \param n ...
259 : !> \param k ...
260 : !> \return ...
261 : !> \author Ole Schuett
262 : ! **************************************************************************************************
263 6150 : PURE FUNCTION multinomial(n, k) RESULT(res)
264 : INTEGER, INTENT(IN) :: n
265 : INTEGER, DIMENSION(:), INTENT(IN) :: k
266 : REAL(KIND=dp) :: res
267 :
268 : INTEGER :: i
269 : REAL(KIND=dp) :: denom
270 :
271 49200 : IF (ALL(k >= 0) .AND. SUM(k) == n) THEN
272 3156 : denom = 1.0_dp
273 12624 : DO i = 1, SIZE(k)
274 12624 : denom = denom*fac(k(i))
275 : END DO
276 3156 : res = fac(n)/denom
277 : ELSE
278 : res = 0.0_dp
279 : END IF
280 :
281 6150 : END FUNCTION multinomial
282 :
283 : ! **************************************************************************************************
284 : !> \brief The rotation matrix rotmat which rotates a vector about a
285 : !> rotation axis defined by the vector a is build up.
286 : !> The rotation angle is phi (radians).
287 : !> \param phi ...
288 : !> \param a ...
289 : !> \param rotmat ...
290 : !> \date 16.10.1998
291 : !> \author MK
292 : !> \version 1.0
293 : ! **************************************************************************************************
294 15759 : PURE SUBROUTINE build_rotmat(phi, a, rotmat)
295 : REAL(KIND=dp), INTENT(IN) :: phi
296 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a
297 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: rotmat
298 :
299 : REAL(KIND=dp) :: cosp, cost, length_of_a, sinp
300 : REAL(KIND=dp), DIMENSION(3) :: d
301 :
302 15759 : length_of_a = SQRT(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
303 : ! Check the length of the vector a
304 15759 : IF (length_of_a > eps_geo) THEN
305 :
306 63036 : d(:) = a(:)/length_of_a
307 :
308 15759 : cosp = COS(phi)
309 15759 : sinp = SIN(phi)
310 15759 : cost = 1.0_dp - cosp
311 :
312 15759 : rotmat(1, 1) = d(1)*d(1)*cost + cosp
313 15759 : rotmat(1, 2) = d(1)*d(2)*cost - d(3)*sinp
314 15759 : rotmat(1, 3) = d(1)*d(3)*cost + d(2)*sinp
315 15759 : rotmat(2, 1) = d(2)*d(1)*cost + d(3)*sinp
316 15759 : rotmat(2, 2) = d(2)*d(2)*cost + cosp
317 15759 : rotmat(2, 3) = d(2)*d(3)*cost - d(1)*sinp
318 15759 : rotmat(3, 1) = d(3)*d(1)*cost - d(2)*sinp
319 15759 : rotmat(3, 2) = d(3)*d(2)*cost + d(1)*sinp
320 15759 : rotmat(3, 3) = d(3)*d(3)*cost + cosp
321 : ELSE
322 0 : CALL unit_matrix(rotmat)
323 : END IF
324 :
325 15759 : END SUBROUTINE build_rotmat
326 :
327 : ! **************************************************************************************************
328 : !> \brief Returns the determinante of the 3x3 matrix a.
329 : !> \param a ...
330 : !> \return ...
331 : !> \date 13.03.2004
332 : !> \author MK
333 : !> \version 1.0
334 : ! **************************************************************************************************
335 16588056 : PURE FUNCTION det_3x3_1(a) RESULT(det_a)
336 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
337 : REAL(KIND=dp) :: det_a
338 :
339 : det_a = a(1, 1)*(a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2)) + &
340 : a(1, 2)*(a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3)) + &
341 16588056 : a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
342 :
343 16588056 : END FUNCTION det_3x3_1
344 :
345 : ! **************************************************************************************************
346 : !> \brief Returns the determinante of the 3x3 matrix a given by its columns.
347 : !> \param a1 ...
348 : !> \param a2 ...
349 : !> \param a3 ...
350 : !> \return ...
351 : !> \date 13.03.2004
352 : !> \author MK
353 : !> \version 1.0
354 : ! **************************************************************************************************
355 5665 : PURE FUNCTION det_3x3_2(a1, a2, a3) RESULT(det_a)
356 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3
357 : REAL(KIND=dp) :: det_a
358 :
359 : det_a = a1(1)*(a2(2)*a3(3) - a3(2)*a2(3)) + &
360 : a2(1)*(a3(2)*a1(3) - a1(2)*a3(3)) + &
361 5665 : a3(1)*(a1(2)*a2(3) - a2(2)*a1(3))
362 :
363 5665 : END FUNCTION det_3x3_2
364 :
365 : ! **************************************************************************************************
366 : !> \brief Diagonalize the symmetric n by n matrix a using the LAPACK
367 : !> library. Only the upper triangle of matrix a is used.
368 : !> Externals (LAPACK 3.0)
369 : !> \param a ...
370 : !> \param eigval ...
371 : !> \param dac ...
372 : !> \date 29.03.1999
373 : !> \par Variables
374 : !> - a : Symmetric matrix to be diagonalized (input; upper triangle) ->
375 : !> - eigenvectors of the matrix a (output).
376 : !> - dac : If true, then the divide-and-conquer algorithm is applied.
377 : !> - eigval : Eigenvalues of the matrix a (output).
378 : !> \author MK
379 : !> \version 1.0
380 : ! **************************************************************************************************
381 77855 : SUBROUTINE diamat_all(a, eigval, dac)
382 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
383 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigval
384 : LOGICAL, INTENT(IN), OPTIONAL :: dac
385 :
386 : CHARACTER(len=*), PARAMETER :: routineN = 'diamat_all'
387 :
388 : INTEGER :: handle, info, liwork, lwork, n, nb
389 77855 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
390 : INTEGER, EXTERNAL :: ilaenv
391 : LOGICAL :: divide_and_conquer
392 77855 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
393 :
394 : EXTERNAL dsyev, dsyevd
395 :
396 77855 : CALL timeset(routineN, handle)
397 :
398 : ! Get the size of the matrix a
399 77855 : n = SIZE(a, 1)
400 :
401 : ! Check the size of matrix a
402 77855 : IF (SIZE(a, 2) /= n) THEN
403 0 : CPABORT("Check the size of matrix a (parameter #1)")
404 : END IF
405 :
406 : ! Check the size of vector eigval
407 77855 : IF (SIZE(eigval) /= n) THEN
408 0 : CPABORT("The dimension of vector eigval is too small")
409 : END IF
410 :
411 : ! Check, if the divide-and-conquer algorithm is requested
412 :
413 77855 : IF (PRESENT(dac)) THEN
414 275 : divide_and_conquer = dac
415 : ELSE
416 : divide_and_conquer = .FALSE.
417 : END IF
418 :
419 : ! Get the optimal work storage size
420 :
421 275 : IF (divide_and_conquer) THEN
422 275 : lwork = 2*n**2 + 6*n + 1
423 275 : liwork = 5*n + 3
424 : ELSE
425 77580 : nb = ilaenv(1, "DSYTRD", "U", n, -1, -1, -1)
426 77580 : lwork = (nb + 2)*n
427 : END IF
428 :
429 : ! Allocate work storage
430 :
431 233565 : ALLOCATE (work(lwork))
432 77855 : IF (divide_and_conquer) THEN
433 825 : ALLOCATE (iwork(liwork))
434 : END IF
435 :
436 : ! Diagonalize the matrix a
437 :
438 77855 : info = 0
439 : IF (divide_and_conquer) THEN
440 275 : CALL dsyevd("V", "U", n, a, n, eigval, work, lwork, iwork, liwork, info)
441 : ELSE
442 79404 : CALL dsyev("V", "U", n, a, n, eigval, work, lwork, info)
443 : END IF
444 :
445 77855 : IF (info /= 0) THEN
446 0 : IF (divide_and_conquer) THEN
447 0 : CPABORT("The matrix diagonalization with dsyevd failed")
448 : ELSE
449 0 : CPABORT("The matrix diagonalization with dsyev failed")
450 : END IF
451 : END IF
452 :
453 : ! Release work storage
454 77855 : DEALLOCATE (work)
455 :
456 77855 : IF (divide_and_conquer) THEN
457 275 : DEALLOCATE (iwork)
458 : END IF
459 :
460 77855 : CALL timestop(handle)
461 :
462 155710 : END SUBROUTINE diamat_all
463 :
464 : ! **************************************************************************************************
465 : !> \brief Returns the dihedral angle, i.e. the angle between the planes
466 : !> defined by the vectors (-ab,bc) and (cd,-bc).
467 : !> The dihedral angle is returned in radians.
468 : !> \param ab ...
469 : !> \param bc ...
470 : !> \param cd ...
471 : !> \return ...
472 : !> \date 13.03.2004
473 : !> \author MK
474 : !> \version 1.0
475 : ! **************************************************************************************************
476 117 : PURE FUNCTION dihedral_angle(ab, bc, cd) RESULT(dihedral_angle_abcd)
477 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ab, bc, cd
478 : REAL(KIND=dp) :: dihedral_angle_abcd
479 :
480 : REAL(KIND=dp) :: det_abcd
481 : REAL(KIND=dp), DIMENSION(3) :: abc, bcd
482 :
483 468 : abc = vector_product(bc, -ab)
484 468 : bcd = vector_product(cd, -bc)
485 : ! Calculate the normal vectors of the planes
486 : ! defined by the points a,b,c and b,c,d
487 :
488 468 : det_abcd = det_3x3(abc, bcd, -bc)
489 117 : dihedral_angle_abcd = SIGN(1.0_dp, det_abcd)*angle(abc, bcd)
490 :
491 117 : END FUNCTION dihedral_angle
492 :
493 : ! **************************************************************************************************
494 : !> \brief Return the diagonal elements of matrix a as a vector.
495 : !> \param a ...
496 : !> \return ...
497 : !> \date 20.11.1998
498 : !> \author MK
499 : !> \version 1.0
500 : ! **************************************************************************************************
501 63502 : PURE FUNCTION get_diag(a) RESULT(a_diag)
502 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
503 : REAL(KIND=dp), &
504 : DIMENSION(MIN(SIZE(a, 1), SIZE(a, 2))) :: a_diag
505 :
506 : INTEGER :: i, n
507 :
508 63502 : n = MIN(SIZE(a, 1), SIZE(a, 2))
509 :
510 6369514 : DO i = 1, n
511 6369514 : a_diag(i) = a(i, i)
512 : END DO
513 :
514 63502 : END FUNCTION get_diag
515 :
516 : ! **************************************************************************************************
517 : !> \brief Returns the inverse of the 3 x 3 matrix a.
518 : !> \param a ...
519 : !> \return ...
520 : !> \date 13.03.2004
521 : !> \author MK
522 : !> \version 1.0
523 : ! **************************************************************************************************
524 15784313 : PURE FUNCTION inv_3x3(a) RESULT(a_inv)
525 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: a
526 : REAL(KIND=dp), DIMENSION(3, 3) :: a_inv
527 :
528 : REAL(KIND=dp) :: det_a
529 :
530 15784313 : det_a = 1.0_dp/det_3x3(a)
531 :
532 15784313 : a_inv(1, 1) = (a(2, 2)*a(3, 3) - a(3, 2)*a(2, 3))*det_a
533 15784313 : a_inv(2, 1) = (a(2, 3)*a(3, 1) - a(3, 3)*a(2, 1))*det_a
534 15784313 : a_inv(3, 1) = (a(2, 1)*a(3, 2) - a(3, 1)*a(2, 2))*det_a
535 :
536 15784313 : a_inv(1, 2) = (a(1, 3)*a(3, 2) - a(3, 3)*a(1, 2))*det_a
537 15784313 : a_inv(2, 2) = (a(1, 1)*a(3, 3) - a(3, 1)*a(1, 3))*det_a
538 15784313 : a_inv(3, 2) = (a(1, 2)*a(3, 1) - a(3, 2)*a(1, 1))*det_a
539 :
540 15784313 : a_inv(1, 3) = (a(1, 2)*a(2, 3) - a(2, 2)*a(1, 3))*det_a
541 15784313 : a_inv(2, 3) = (a(1, 3)*a(2, 1) - a(2, 3)*a(1, 1))*det_a
542 15784313 : a_inv(3, 3) = (a(1, 1)*a(2, 2) - a(2, 1)*a(1, 2))*det_a
543 :
544 15784313 : END FUNCTION inv_3x3
545 :
546 : ! **************************************************************************************************
547 : !> \brief returns inverse of matrix using the lapack routines DGETRF and DGETRI
548 : !> \param a ...
549 : !> \param info ...
550 : ! **************************************************************************************************
551 9269 : SUBROUTINE invmat(a, info)
552 : REAL(KIND=dp), INTENT(INOUT) :: a(:, :)
553 : INTEGER, INTENT(OUT) :: info
554 :
555 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invmat'
556 :
557 : INTEGER :: handle, lwork, n
558 9269 : INTEGER, ALLOCATABLE :: ipiv(:)
559 9269 : REAL(KIND=dp), ALLOCATABLE :: work(:)
560 :
561 9269 : CALL timeset(routineN, handle)
562 :
563 9269 : n = SIZE(a, 1)
564 9269 : lwork = 20*n
565 27807 : ALLOCATE (ipiv(n))
566 27807 : ALLOCATE (work(lwork))
567 9269 : ipiv = 0
568 9269 : work = 0._dp
569 9269 : info = 0
570 487485 : CALL dgetrf(n, n, a, n, ipiv, info)
571 9269 : IF (info == 0) THEN
572 486335 : CALL dgetri(n, a, n, ipiv, work, lwork, info)
573 : END IF
574 9269 : DEALLOCATE (ipiv, work)
575 :
576 9269 : CALL timestop(handle)
577 :
578 9269 : END SUBROUTINE invmat
579 :
580 : ! **************************************************************************************************
581 : !> \brief returns inverse of real symmetric, positive definite matrix
582 : !> \param a matrix
583 : !> \param potrf if cholesky decomposition of a was already done using dpotrf.
584 : !> If not given, cholesky decomposition of a will be done before inversion.
585 : !> \param uplo indicating if the upper or lower triangle of a is stored.
586 : !> \author Dorothea Golze [02.2015]
587 : ! **************************************************************************************************
588 3075 : SUBROUTINE invmat_symm(a, potrf, uplo)
589 : REAL(KIND=dp), INTENT(INOUT) :: a(:, :)
590 : LOGICAL, INTENT(IN), OPTIONAL :: potrf
591 : CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: uplo
592 :
593 : CHARACTER(LEN=*), PARAMETER :: routineN = 'invmat_symm'
594 :
595 : CHARACTER(LEN=1) :: myuplo
596 : INTEGER :: handle, info, n
597 : LOGICAL :: do_potrf
598 :
599 3075 : CALL timeset(routineN, handle)
600 :
601 3075 : do_potrf = .TRUE.
602 3075 : IF (PRESENT(potrf)) do_potrf = potrf
603 :
604 3075 : myuplo = 'U'
605 3075 : IF (PRESENT(uplo)) myuplo = uplo
606 :
607 3075 : n = SIZE(a, 1)
608 3075 : info = 0
609 :
610 : ! do cholesky decomposition
611 3075 : IF (do_potrf) THEN
612 4058 : CALL dpotrf(myuplo, n, a, n, info)
613 2858 : IF (info /= 0) CPABORT("DPOTRF failed")
614 : END IF
615 :
616 : ! do inversion using the cholesky decomposition
617 4275 : CALL dpotri(myuplo, n, a, n, info)
618 3075 : IF (info /= 0) CPABORT("Matrix inversion failed")
619 :
620 : ! complete the matrix
621 3075 : IF ((myuplo == "U") .OR. (myuplo == "u")) THEN
622 3075 : CALL symmetrize_matrix(a, "upper_to_lower")
623 : ELSE
624 0 : CALL symmetrize_matrix(a, "lower_to_upper")
625 : END IF
626 :
627 3075 : CALL timestop(handle)
628 :
629 3075 : END SUBROUTINE invmat_symm
630 :
631 : ! **************************************************************************************************
632 : !> \brief Compute the inverse of the n by n real matrix a using the LAPACK
633 : !> library
634 : !> \param a ...
635 : !> \param a_inverse ...
636 : !> \param eval_error ...
637 : !> \param option ...
638 : !> \param improve ...
639 : !> \date 23.03.1999
640 : !> \par Variables
641 : !> - a : Real matrix to be inverted (input).
642 : !> - a_inverse: Inverse of the matrix a (output).
643 : !> - a_lu : LU factorization of matrix a.
644 : !> - a_norm : Norm of matrix a.
645 : !> - error : Estimated error of the inversion.
646 : !> - r_cond : Reciprocal condition number of the matrix a.
647 : !> - trans : "N" => invert a
648 : !> - "T" => invert transpose(a)
649 : !> \author MK
650 : !> \version 1.0
651 : !> \note NB add improve argument, used to disable call to dgerfs
652 : ! **************************************************************************************************
653 11266 : SUBROUTINE invert_matrix_d(a, a_inverse, eval_error, option, improve)
654 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
655 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: a_inverse
656 : REAL(KIND=dp), INTENT(OUT) :: eval_error
657 : CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: option
658 : LOGICAL, INTENT(IN), OPTIONAL :: improve
659 :
660 : CHARACTER(LEN=1) :: norm, trans
661 : CHARACTER(LEN=default_string_length) :: message
662 : INTEGER :: info, iter, n
663 11266 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv, iwork
664 : LOGICAL :: do_improve
665 : REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
666 11266 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, work
667 11266 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
668 : REAL(KIND=dp), EXTERNAL :: dlange
669 :
670 : EXTERNAL dgecon, dgerfs, dgetrf, dgetrs
671 :
672 : ! Check for optional parameter
673 11266 : IF (PRESENT(option)) THEN
674 580 : trans = option
675 : ELSE
676 10686 : trans = "N"
677 : END IF
678 :
679 11266 : IF (PRESENT(improve)) THEN
680 1088 : do_improve = improve
681 : ELSE
682 : do_improve = .TRUE.
683 : END IF
684 :
685 : ! Get the dimension of matrix a
686 11266 : n = SIZE(a, 1)
687 :
688 : ! Check array dimensions
689 11266 : IF (n == 0) THEN
690 0 : CPABORT("Matrix to be inverted of zero size")
691 : END IF
692 :
693 11266 : IF (n /= SIZE(a, 2)) THEN
694 0 : CPABORT("Check the array bounds of parameter #1")
695 : END IF
696 :
697 11266 : IF ((n /= SIZE(a_inverse, 1)) .OR. &
698 : (n /= SIZE(a_inverse, 2))) THEN
699 0 : CPABORT("Check the array bounds of parameter #2")
700 : END IF
701 :
702 : ! Allocate work storage
703 45064 : ALLOCATE (a_lu(n, n))
704 33798 : ALLOCATE (b(n, n))
705 33798 : ALLOCATE (berr(n))
706 22532 : ALLOCATE (ferr(n))
707 33798 : ALLOCATE (ipiv(n))
708 22532 : ALLOCATE (iwork(n))
709 33798 : ALLOCATE (work(4*n))
710 :
711 2900102 : a_lu(1:n, 1:n) = a(1:n, 1:n)
712 :
713 : ! Compute the LU factorization of the matrix a
714 11266 : CALL dgetrf(n, n, a_lu, n, ipiv, info)
715 :
716 11266 : IF (info /= 0) THEN
717 0 : CPABORT("The LU factorization in dgetrf failed")
718 : END IF
719 :
720 : ! Compute the norm of the matrix a
721 :
722 11266 : IF (trans == "N") THEN
723 10988 : norm = '1'
724 : ELSE
725 278 : norm = 'I'
726 : END IF
727 :
728 11266 : a_norm = dlange(norm, n, n, a, n, work)
729 :
730 : ! Compute the reciprocal of the condition number of a
731 :
732 11266 : CALL dgecon(norm, n, a_lu, n, a_norm, r_cond, work, iwork, info)
733 :
734 11266 : IF (info /= 0) THEN
735 0 : CPABORT("The computation of the condition number in dgecon failed")
736 : END IF
737 :
738 11266 : IF (r_cond < EPSILON(0.0_dp)) THEN
739 0 : WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
740 : CALL cp_abort(__LOCATION__, &
741 : "Bad condition number "//TRIM(message)//" (smaller than the machine "// &
742 0 : "working precision)")
743 : END IF
744 :
745 : ! Solve a system of linear equations using the LU factorization computed by dgetrf
746 :
747 11266 : CALL unit_matrix(a_inverse)
748 :
749 11266 : CALL dgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
750 :
751 11266 : IF (info /= 0) THEN
752 0 : CPABORT("Solving the system of linear equations in dgetrs failed")
753 : END IF
754 :
755 : ! Improve the computed solution iteratively
756 11266 : CALL unit_matrix(b) ! Initialize right-hand sides
757 :
758 11266 : eval_error = 0.0_dp
759 :
760 11266 : IF (do_improve) THEN
761 25526 : DO iter = 1, 10
762 :
763 : CALL dgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
764 25266 : work, iwork, info)
765 :
766 25266 : IF (info /= 0) THEN
767 0 : CPABORT("Improving the computed solution in dgerfs failed")
768 : END IF
769 :
770 25266 : old_eval_error = eval_error
771 280060 : eval_error = MAXVAL(ferr)
772 :
773 25526 : IF (ABS(eval_error - old_eval_error) <= EPSILON(1.0_dp)) EXIT
774 :
775 : END DO
776 : END IF
777 :
778 : ! Release work storage
779 11266 : DEALLOCATE (work)
780 11266 : DEALLOCATE (iwork)
781 11266 : DEALLOCATE (ipiv)
782 11266 : DEALLOCATE (ferr)
783 11266 : DEALLOCATE (berr)
784 11266 : DEALLOCATE (b)
785 11266 : DEALLOCATE (a_lu)
786 :
787 11266 : END SUBROUTINE invert_matrix_d
788 :
789 : ! **************************************************************************************************
790 : !> \brief Compute the inverse of the n by n complex matrix a using the LAPACK
791 : !> library
792 : !> \param a ...
793 : !> \param a_inverse ...
794 : !> \param eval_error ...
795 : !> \param option ...
796 : !> \date 08.06.2009
797 : !> \par Variables
798 : !> - a : Complex matrix to be inverted (input).
799 : !> - a_inverse: Inverse of the matrix a (output).
800 : !> - a_lu : LU factorization of matrix a.
801 : !> - a_norm : Norm of matrix a.
802 : !> - error : Estimated error of the inversion.
803 : !> - r_cond : Reciprocal condition number of the matrix a.
804 : !> - trans : "N" => invert a
805 : !> - "T" => invert transpose(a)
806 : !> \author MK
807 : !> \version 1.0
808 : ! **************************************************************************************************
809 0 : SUBROUTINE invert_matrix_z(a, a_inverse, eval_error, option)
810 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
811 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: a_inverse
812 : REAL(KIND=dp), INTENT(OUT) :: eval_error
813 : CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: option
814 :
815 : CHARACTER(LEN=1) :: norm, trans
816 : CHARACTER(LEN=default_string_length) :: message
817 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
818 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
819 : INTEGER :: info, iter, n
820 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
821 : REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
822 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, rwork
823 : REAL(KIND=dp), EXTERNAL :: zlange
824 :
825 : EXTERNAL zgecon, zgerfs, zgetrf, zgetrs
826 :
827 : ! Check for optional parameter
828 0 : IF (PRESENT(option)) THEN
829 0 : trans = option
830 : ELSE
831 0 : trans = "N"
832 : END IF
833 :
834 : ! Get the dimension of matrix a
835 0 : n = SIZE(a, 1)
836 :
837 : ! Check array dimensions
838 0 : IF (n == 0) THEN
839 0 : CPABORT("Matrix to be inverted of zero size")
840 : END IF
841 :
842 0 : IF (n /= SIZE(a, 2)) THEN
843 0 : CPABORT("Check the array bounds of parameter #1")
844 : END IF
845 :
846 0 : IF ((n /= SIZE(a_inverse, 1)) .OR. &
847 : (n /= SIZE(a_inverse, 2))) THEN
848 0 : CPABORT("Check the array bounds of parameter #2")
849 : END IF
850 :
851 : ! Allocate work storage
852 0 : ALLOCATE (a_lu(n, n))
853 0 : ALLOCATE (b(n, n))
854 0 : ALLOCATE (berr(n))
855 0 : ALLOCATE (ferr(n))
856 0 : ALLOCATE (ipiv(n))
857 0 : ALLOCATE (rwork(2*n))
858 0 : ALLOCATE (work(2*n))
859 :
860 0 : a_lu(1:n, 1:n) = a(1:n, 1:n)
861 :
862 : ! Compute the LU factorization of the matrix a
863 0 : CALL zgetrf(n, n, a_lu, n, ipiv, info)
864 :
865 0 : IF (info /= 0) THEN
866 0 : CPABORT("The LU factorization in dgetrf failed")
867 : END IF
868 :
869 : ! Compute the norm of the matrix a
870 :
871 0 : IF (trans == "N") THEN
872 0 : norm = '1'
873 : ELSE
874 0 : norm = 'I'
875 : END IF
876 :
877 0 : a_norm = zlange(norm, n, n, a, n, work)
878 :
879 : ! Compute the reciprocal of the condition number of a
880 :
881 0 : CALL zgecon(norm, n, a_lu, n, a_norm, r_cond, work, rwork, info)
882 :
883 0 : IF (info /= 0) THEN
884 0 : CPABORT("The computation of the condition number in dgecon failed")
885 : END IF
886 :
887 0 : IF (r_cond < EPSILON(0.0_dp)) THEN
888 0 : WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
889 : CALL cp_abort(__LOCATION__, &
890 : "Bad condition number "//TRIM(message)//" (smaller than the machine "// &
891 0 : "working precision)")
892 : END IF
893 :
894 : ! Solve a system of linear equations using the LU factorization computed by dgetrf
895 :
896 0 : CALL unit_matrix(a_inverse)
897 :
898 0 : CALL zgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
899 :
900 0 : IF (info /= 0) THEN
901 0 : CPABORT("Solving the system of linear equations in dgetrs failed")
902 : END IF
903 :
904 : ! Improve the computed solution iteratively
905 0 : CALL unit_matrix(b) ! Initialize right-hand sides
906 :
907 0 : eval_error = 0.0_dp
908 :
909 0 : DO iter = 1, 10
910 :
911 : CALL zgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
912 0 : work, rwork, info)
913 :
914 0 : IF (info /= 0) THEN
915 0 : CPABORT("Improving the computed solution in dgerfs failed")
916 : END IF
917 :
918 0 : old_eval_error = eval_error
919 0 : eval_error = MAXVAL(ferr)
920 :
921 0 : IF (ABS(eval_error - old_eval_error) <= EPSILON(1.0_dp)) EXIT
922 :
923 : END DO
924 :
925 : ! Release work storage
926 0 : DEALLOCATE (work)
927 0 : DEALLOCATE (rwork)
928 0 : DEALLOCATE (ipiv)
929 0 : DEALLOCATE (ferr)
930 0 : DEALLOCATE (berr)
931 0 : DEALLOCATE (b)
932 0 : DEALLOCATE (a_lu)
933 :
934 0 : END SUBROUTINE invert_matrix_z
935 :
936 : ! **************************************************************************************************
937 : !> \brief returns the pseudoinverse of a real, square matrix using singular
938 : !> value decomposition
939 : !> \param a matrix a
940 : !> \param a_pinverse pseudoinverse of matrix a
941 : !> \param rskip parameter for setting small singular values to zero
942 : !> \param determinant determinant of matrix a (optional output)
943 : !> \param sval array holding singular values of matrix a (optional output)
944 : !> \author Dorothea Golze [02.2015]
945 : ! **************************************************************************************************
946 1505 : SUBROUTINE get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
947 : REAL(KIND=dp), DIMENSION(:, :) :: a, a_pinverse
948 : REAL(KIND=dp), INTENT(IN) :: rskip
949 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: determinant
950 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT), &
951 : OPTIONAL, POINTER :: sval
952 :
953 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_pseudo_inverse_svd'
954 :
955 : INTEGER :: handle, i, info, lwork, n
956 1505 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
957 1505 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sig, work
958 1505 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sig_plus, temp_mat, u, vt
959 :
960 1505 : CALL timeset(routineN, handle)
961 :
962 1505 : n = SIZE(a, 1)
963 21070 : ALLOCATE (u(n, n), vt(n, n), sig(n), sig_plus(n, n), iwork(8*n), work(1), temp_mat(n, n))
964 1505 : u(:, :) = 0.0_dp
965 1505 : vt(:, :) = 0.0_dp
966 1505 : sig(:) = 0.0_dp
967 1505 : sig_plus = 0.0_dp
968 1505 : work = 0.0_dp
969 1505 : iwork = 0
970 1505 : IF (PRESENT(determinant)) determinant = 1.0_dp
971 :
972 : ! work size query
973 1505 : lwork = -1
974 : CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
975 1505 : lwork, iwork(1), info)
976 :
977 1505 : IF (info /= 0) THEN
978 0 : CPABORT("ERROR in DGESDD: Could not retrieve work array sizes")
979 : END IF
980 1505 : lwork = INT(work(1))
981 1505 : DEALLOCATE (work)
982 4515 : ALLOCATE (work(lwork))
983 :
984 : ! do SVD
985 : CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
986 1505 : lwork, iwork(1), info)
987 :
988 1505 : IF (info /= 0) THEN
989 0 : CPABORT("SVD failed")
990 : END IF
991 :
992 1505 : IF (PRESENT(sval)) THEN
993 0 : CPASSERT(.NOT. ASSOCIATED(sval))
994 0 : ALLOCATE (sval(n))
995 0 : sval(:) = sig
996 : END IF
997 :
998 : ! set singular values that are too small to zero
999 9567 : DO i = 1, n
1000 92401 : IF (sig(i) > rskip*MAXVAL(sig)) THEN
1001 7942 : IF (PRESENT(determinant)) THEN
1002 0 : determinant = determinant*sig(i)
1003 : END IF
1004 7942 : sig_plus(i, i) = 1._dp/sig(i)
1005 : ELSE
1006 120 : sig_plus(i, i) = 0.0_dp
1007 : END IF
1008 : END DO
1009 :
1010 : ! build pseudoinverse: V*sig_plus*UT
1011 1505 : CALL dgemm("N", "T", n, n, n, 1._dp, sig_plus, n, u, n, 0._dp, temp_mat, n)
1012 1505 : CALL dgemm("T", "N", n, n, n, 1._dp, vt, n, temp_mat, n, 0._dp, a_pinverse, n)
1013 :
1014 1505 : DEALLOCATE (u, vt, sig, iwork, work, sig_plus, temp_mat)
1015 :
1016 1505 : CALL timestop(handle)
1017 :
1018 1505 : END SUBROUTINE get_pseudo_inverse_svd
1019 :
1020 : ! **************************************************************************************************
1021 : !> \brief returns the pseudoinverse of a real, symmetric and positive definite
1022 : !> matrix using diagonalization.
1023 : !> \param a matrix a
1024 : !> \param a_pinverse pseudoinverse of matrix a
1025 : !> \param rskip parameter for setting small eigenvalues to zero
1026 : !> \author Dorothea Golze [02.2015]
1027 : ! **************************************************************************************************
1028 2287 : SUBROUTINE get_pseudo_inverse_diag(a, a_pinverse, rskip)
1029 : REAL(KIND=dp), DIMENSION(:, :) :: a, a_pinverse
1030 : REAL(KIND=dp), INTENT(IN) :: rskip
1031 :
1032 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_pseudo_inverse_diag'
1033 :
1034 : INTEGER :: handle, i, info, lwork, n
1035 2287 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
1036 2287 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dinv, temp_mat
1037 :
1038 2287 : CALL timeset(routineN, handle)
1039 :
1040 2287 : info = 0
1041 2287 : n = SIZE(a, 1)
1042 18296 : ALLOCATE (dinv(n, n), eig(n), work(1), temp_mat(n, n))
1043 2287 : dinv(:, :) = 0.0_dp
1044 2287 : eig(:) = 0.0_dp
1045 2287 : work(:) = 0.0_dp
1046 2287 : temp_mat = 0.0_dp
1047 :
1048 : ! work size query
1049 2287 : lwork = -1
1050 2287 : CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1051 2287 : IF (info /= 0) THEN
1052 0 : CPABORT("ERROR in DSYEV: Could not retrieve work array sizes")
1053 : END IF
1054 2287 : lwork = INT(work(1))
1055 2287 : DEALLOCATE (work)
1056 6861 : ALLOCATE (work(lwork))
1057 2287 : work = 0.0_dp
1058 :
1059 : ! get eigenvalues and eigenvectors
1060 2287 : CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1061 :
1062 2287 : IF (info /= 0) THEN
1063 0 : CPABORT("Matrix diagonalization failed")
1064 : END IF
1065 :
1066 : ! set eigenvalues that are too small to zero
1067 29641 : DO i = 1, n
1068 361789 : IF (eig(i) > rskip*MAXVAL(eig)) THEN
1069 24682 : dinv(i, i) = 1.0_dp/eig(i)
1070 : ELSE
1071 2672 : dinv(i, i) = 0._dp
1072 : END IF
1073 : END DO
1074 :
1075 : ! build pseudoinverse: U*dinv*UT
1076 2287 : CALL dgemm("N", "T", n, n, n, 1._dp, dinv, n, a, n, 0._dp, temp_mat, n)
1077 2287 : CALL dgemm("N", "N", n, n, n, 1._dp, a, n, temp_mat, n, 0._dp, a_pinverse, n)
1078 :
1079 2287 : DEALLOCATE (eig, work, dinv, temp_mat)
1080 :
1081 2287 : CALL timestop(handle)
1082 :
1083 2287 : END SUBROUTINE get_pseudo_inverse_diag
1084 :
1085 : ! **************************************************************************************************
1086 : !> \brief Reflection of the vector a through a mirror plane defined by the
1087 : !> normal vector b. The reflected vector a is stored in a_mirror.
1088 : !> \param a ...
1089 : !> \param b ...
1090 : !> \return ...
1091 : !> \date 16.10.1998
1092 : !> \author MK
1093 : !> \version 1.0
1094 : ! **************************************************************************************************
1095 19487 : PURE FUNCTION reflect_vector(a, b) RESULT(a_mirror)
1096 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a, b
1097 : REAL(KIND=dp), DIMENSION(3) :: a_mirror
1098 :
1099 : REAL(KIND=dp) :: length_of_b, scapro
1100 : REAL(KIND=dp), DIMENSION(3) :: d
1101 :
1102 19487 : length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1103 :
1104 19487 : IF (length_of_b > eps_geo) THEN
1105 :
1106 77948 : d(:) = b(:)/length_of_b
1107 :
1108 : ! Calculate the mirror image a_mirror of the vector a
1109 19487 : scapro = a(1)*d(1) + a(2)*d(2) + a(3)*d(3)
1110 :
1111 77948 : a_mirror(:) = a(:) - 2.0_dp*scapro*d(:)
1112 :
1113 : ELSE
1114 :
1115 0 : a_mirror(:) = 0.0_dp
1116 :
1117 : END IF
1118 :
1119 19487 : END FUNCTION reflect_vector
1120 :
1121 : ! **************************************************************************************************
1122 : !> \brief Rotation of the vector a about an rotation axis defined by the
1123 : !> vector b. The rotation angle is phi (radians). The rotated vector
1124 : !> a is stored in a_rot.
1125 : !> \param a ...
1126 : !> \param phi ...
1127 : !> \param b ...
1128 : !> \return ...
1129 : !> \date 16.10.1998
1130 : !> \author MK
1131 : !> \version 1.0
1132 : ! **************************************************************************************************
1133 4098 : PURE FUNCTION rotate_vector(a, phi, b) RESULT(a_rot)
1134 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a
1135 : REAL(KIND=dp), INTENT(IN) :: phi
1136 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: b
1137 : REAL(KIND=dp), DIMENSION(3) :: a_rot
1138 :
1139 : REAL(KIND=dp) :: length_of_b
1140 : REAL(KIND=dp), DIMENSION(3, 3) :: rotmat
1141 :
1142 4098 : length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1143 4098 : IF (length_of_b > eps_geo) THEN
1144 :
1145 : ! Build up the rotation matrix rotmat
1146 3964 : CALL build_rotmat(phi, b, rotmat)
1147 :
1148 : ! Rotate the vector a by phi about the axis defined by vector b
1149 63424 : a_rot(:) = MATMUL(rotmat, a)
1150 :
1151 : ELSE
1152 :
1153 536 : a_rot(:) = 0.0_dp
1154 :
1155 : END IF
1156 :
1157 4098 : END FUNCTION rotate_vector
1158 :
1159 : ! **************************************************************************************************
1160 : !> \brief Set the diagonal elements of matrix a to b.
1161 : !> \param a ...
1162 : !> \param b ...
1163 : !> \date 20.11.1998
1164 : !> \author MK
1165 : !> \version 1.0
1166 : ! **************************************************************************************************
1167 59838 : PURE SUBROUTINE set_diag_scalar_d(a, b)
1168 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1169 : REAL(KIND=dp), INTENT(IN) :: b
1170 :
1171 : INTEGER :: i, n
1172 :
1173 59838 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1174 312193 : DO i = 1, n
1175 312193 : a(i, i) = b
1176 : END DO
1177 :
1178 59838 : END SUBROUTINE set_diag_scalar_d
1179 :
1180 : ! **************************************************************************************************
1181 : !> \brief ...
1182 : !> \param a ...
1183 : !> \param b ...
1184 : ! **************************************************************************************************
1185 0 : PURE SUBROUTINE set_diag_scalar_z(a, b)
1186 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1187 : COMPLEX(KIND=dp), INTENT(IN) :: b
1188 :
1189 : INTEGER :: i, n
1190 :
1191 0 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1192 0 : DO i = 1, n
1193 0 : a(i, i) = b
1194 : END DO
1195 :
1196 0 : END SUBROUTINE set_diag_scalar_z
1197 :
1198 : ! **************************************************************************************************
1199 : !> \brief Symmetrize the matrix a.
1200 : !> \param a ...
1201 : !> \param option ...
1202 : !> \date 16.10.1998
1203 : !> \author MK
1204 : !> \version 1.0
1205 : ! **************************************************************************************************
1206 20149 : SUBROUTINE symmetrize_matrix(a, option)
1207 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1208 : CHARACTER(LEN=*), INTENT(IN) :: option
1209 :
1210 : INTEGER :: i, n
1211 :
1212 20149 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1213 :
1214 20149 : IF (option == "lower_to_upper") THEN
1215 0 : DO i = 1, n - 1
1216 0 : a(i, i + 1:n) = a(i + 1:n, i)
1217 : END DO
1218 20149 : ELSE IF (option == "upper_to_lower") THEN
1219 523491 : DO i = 1, n - 1
1220 29745547 : a(i + 1:n, i) = a(i, i + 1:n)
1221 : END DO
1222 42 : ELSE IF (option == "anti_lower_to_upper") THEN
1223 0 : DO i = 1, n - 1
1224 0 : a(i, i + 1:n) = -a(i + 1:n, i)
1225 : END DO
1226 42 : ELSE IF (option == "anti_upper_to_lower") THEN
1227 660 : DO i = 1, n - 1
1228 5532 : a(i + 1:n, i) = -a(i, i + 1:n)
1229 : END DO
1230 : ELSE
1231 0 : CPABORT("Invalid option <"//TRIM(option)//"> was specified for parameter #2")
1232 : END IF
1233 :
1234 20149 : END SUBROUTINE symmetrize_matrix
1235 :
1236 : ! **************************************************************************************************
1237 : !> \brief Set the matrix a to be a unit matrix.
1238 : !> \param a ...
1239 : !> \date 16.10.1998
1240 : !> \author MK
1241 : !> \version 1.0
1242 : ! **************************************************************************************************
1243 59838 : PURE SUBROUTINE unit_matrix_d(a)
1244 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1245 :
1246 7694482 : a(:, :) = 0.0_dp
1247 59838 : CALL set_diag(a, 1.0_dp)
1248 :
1249 59838 : END SUBROUTINE unit_matrix_d
1250 :
1251 : ! **************************************************************************************************
1252 : !> \brief ...
1253 : !> \param a ...
1254 : ! **************************************************************************************************
1255 0 : PURE SUBROUTINE unit_matrix_z(a)
1256 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1257 :
1258 0 : a(:, :) = (0.0_dp, 0.0_dp)
1259 0 : CALL set_diag(a, (1.0_dp, 0.0_dp))
1260 :
1261 0 : END SUBROUTINE unit_matrix_z
1262 :
1263 : ! **************************************************************************************************
1264 : !> \brief Calculation of the vector product c = a x b.
1265 : !> \param a ...
1266 : !> \param b ...
1267 : !> \return ...
1268 : !> \date 16.10.1998
1269 : !> \author MK
1270 : !> \version 1.0
1271 : ! **************************************************************************************************
1272 8037 : PURE FUNCTION vector_product(a, b) RESULT(c)
1273 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: a, b
1274 : REAL(KIND=dp), DIMENSION(3) :: c
1275 :
1276 8037 : c(1) = a(2)*b(3) - a(3)*b(2)
1277 8037 : c(2) = a(3)*b(1) - a(1)*b(3)
1278 8037 : c(3) = a(1)*b(2) - a(2)*b(1)
1279 :
1280 8037 : END FUNCTION vector_product
1281 :
1282 : ! **************************************************************************************************
1283 : !> \brief computes the greatest common divisor of two number
1284 : !> \param a ...
1285 : !> \param b ...
1286 : !> \return ...
1287 : !> \author Joost VandeVondele
1288 : ! **************************************************************************************************
1289 1186489 : ELEMENTAL FUNCTION gcd(a, b)
1290 : INTEGER, INTENT(IN) :: a, b
1291 : INTEGER :: gcd
1292 :
1293 : INTEGER :: aa, ab, l, rem, s
1294 :
1295 1186489 : aa = ABS(a)
1296 1186489 : ab = ABS(b)
1297 1186489 : IF (aa < ab) THEN
1298 : s = aa
1299 : l = ab
1300 : ELSE
1301 1140551 : s = ab
1302 1140551 : l = aa
1303 : END IF
1304 1186489 : IF (s /= 0) THEN
1305 : DO
1306 1186691 : rem = MOD(l, s)
1307 1186691 : IF (rem == 0) EXIT
1308 : l = s
1309 1186489 : s = rem
1310 : END DO
1311 : GCD = s
1312 : ELSE
1313 : GCD = l
1314 : END IF
1315 1186489 : END FUNCTION gcd
1316 :
1317 : ! **************************************************************************************************
1318 : !> \brief computes the least common multiplier of two numbers
1319 : !> \param a ...
1320 : !> \param b ...
1321 : !> \return ...
1322 : !> \author Joost VandeVondele
1323 : ! **************************************************************************************************
1324 546924 : ELEMENTAL FUNCTION lcm(a, b)
1325 : INTEGER, INTENT(IN) :: a, b
1326 : INTEGER :: lcm
1327 :
1328 : INTEGER :: tmp
1329 :
1330 546924 : tmp = gcd(a, b)
1331 546924 : IF (tmp == 0) THEN
1332 : lcm = 0
1333 : ELSE
1334 : ! could still overflow if the true lcm is larger than maxint
1335 546924 : lcm = ABS((a/tmp)*b)
1336 : END IF
1337 546924 : END FUNCTION lcm
1338 :
1339 : ! **************************************************************************************************
1340 : !> \brief computes the exponential integral
1341 : !> Ei(x) = Int(exp(-x*t)/t,t=1..infinity) x>0
1342 : !> \param x ...
1343 : !> \return ...
1344 : !> \author JGH (adapted from Numerical recipies)
1345 : ! **************************************************************************************************
1346 0 : FUNCTION ei(x)
1347 : REAL(dp) :: x, ei
1348 :
1349 : INTEGER, PARAMETER :: maxit = 100
1350 : REAL(dp), PARAMETER :: eps = EPSILON(0.0_dp), &
1351 : fpmin = TINY(0.0_dp)
1352 :
1353 : INTEGER :: k
1354 : REAL(dp) :: fact, prev, sum1, term
1355 :
1356 0 : IF (x <= 0._dp) THEN
1357 0 : CPABORT("Invalid argument")
1358 : END IF
1359 :
1360 0 : IF (x < fpmin) THEN
1361 0 : ei = LOG(x) + euler
1362 0 : ELSE IF (x <= -LOG(EPS)) THEN
1363 : sum1 = 0._dp
1364 : fact = 1._dp
1365 0 : DO k = 1, maxit
1366 0 : fact = fact*x/REAL(k, dp)
1367 0 : term = fact/REAL(k, dp)
1368 0 : sum1 = sum1 + term
1369 0 : IF (term < eps*sum1) EXIT
1370 : END DO
1371 0 : ei = sum1 + LOG(x) + euler
1372 : ELSE
1373 : sum1 = 0._dp
1374 : term = 1._dp
1375 0 : DO k = 1, maxit
1376 0 : prev = term
1377 0 : term = term*REAL(k, dp)/x
1378 0 : IF (term < eps) EXIT
1379 0 : IF (term < prev) THEN
1380 0 : sum1 = sum1 + term
1381 : ELSE
1382 0 : sum1 = sum1 - prev
1383 0 : EXIT
1384 : END IF
1385 : END DO
1386 0 : ei = EXP(x)*(1._dp + sum1)/x
1387 : END IF
1388 :
1389 0 : END FUNCTION ei
1390 :
1391 : ! **************************************************************************************************
1392 : !> \brief Computes the digamma function, the logarithmic derivative of the gamma function.
1393 : !> \param xx Argument of the digamma function.
1394 : !> \return The digamma function at xx. Returns zero at poles and for arguments outside the
1395 : !> supported numerical range.
1396 : !> \par History
1397 : !> Adapted from the FUNPACK routine PSI, written at Argonne National Laboratory and modified
1398 : !> by A. H. Morris (NSWC).
1399 : !> \note The implementation uses the rational Chebyshev approximations of Cody, Strecok, and
1400 : !> Thacher, Math. Comp. 27, 123-127 (1973), together with the reflection formula.
1401 : ! **************************************************************************************************
1402 7564 : PURE FUNCTION digamma(xx) RESULT(fn_val)
1403 : REAL(dp), INTENT(IN) :: xx
1404 : REAL(dp) :: fn_val
1405 :
1406 : REAL(dp), PARAMETER :: dx0 = 1.461632144968362341262659542325721325e0_dp, p1(7) = [ &
1407 : .895385022981970e-02_dp, .477762828042627e+01_dp, .142441585084029e+03_dp, &
1408 : .118645200713425e+04_dp, .363351846806499e+04_dp, .413810161269013e+04_dp, &
1409 : .130560269827897e+04_dp], p2(4) = [-.212940445131011e+01_dp, -.701677227766759e+01_dp, &
1410 : -.448616543918019e+01_dp, -.648157123766197e+00_dp], piov4 = .785398163397448e0_dp, q1(6) &
1411 : = [.448452573429826e+02_dp, .520752771467162e+03_dp, .221000799247830e+04_dp, &
1412 : .364127349079381e+04_dp, .190831076596300e+04_dp, .691091682714533e-05_dp]
1413 : REAL(dp), PARAMETER :: q2(4) = [.322703493791143e+02_dp, .892920700481861e+02_dp, &
1414 : .546117738103215e+02_dp, .777788548522962e+01_dp]
1415 :
1416 : INTEGER :: i, m, n, nq
1417 : REAL(dp) :: aug, den, sgn, upper, w, x, xmax1, xmx0, &
1418 : xsmall, z
1419 :
1420 : ! piov4 is pi/4 and dx0 is the positive zero of the digamma function. The coefficient sets
1421 : ! p1/q1 approximate digamma(x)/(x - dx0) for 0.5 <= x <= 3, while p2/q2 approximate
1422 : ! digamma(x) - log(x) + 1/(2*x) for x > 3.
1423 :
1424 : ! xmax1 limits integer argument reduction for large negative values and marks the point beyond
1425 : ! which the asymptotic value log(x) is sufficient. For very small x, pi*cot(pi*x) is
1426 : ! approximated by 1/x.
1427 7564 : xmax1 = MIN(REAL(HUGE(0), KIND=dp), 1.0e0_dp/EPSILON(1.0e0_dp))
1428 7564 : xsmall = 1.e-9_dp
1429 :
1430 7564 : x = xx
1431 7564 : aug = 0.0e0_dp
1432 7564 : IF (x < 0.5e0_dp) THEN
1433 : ! Use digamma(x) = digamma(1 - x) - pi*cot(pi*x).
1434 0 : IF (ABS(x) <= xsmall) THEN
1435 0 : IF (x == 0.0e0_dp) THEN
1436 : ! The digamma function has a pole at zero.
1437 7564 : fn_val = 0.0e0_dp
1438 : RETURN
1439 : END IF
1440 : ! For small nonzero x, use -1/x for -pi*cot(pi*x).
1441 0 : aug = -1.0e0_dp/x
1442 0 : x = 1.0e0_dp - x
1443 : ELSE
1444 : ! Reduce the cotangent argument to the first quadrant and keep track of its sign.
1445 0 : w = -x
1446 0 : sgn = piov4
1447 0 : IF (w <= 0.0e0_dp) THEN
1448 0 : w = -w
1449 0 : sgn = -sgn
1450 : END IF
1451 0 : IF (w >= xmax1) THEN
1452 : ! Argument reduction is no longer reliable.
1453 7564 : fn_val = 0.0e0_dp
1454 : RETURN
1455 : END IF
1456 0 : nq = INT(w)
1457 0 : w = w - nq
1458 0 : nq = INT(w*4.0e0_dp)
1459 0 : w = 4.0e0_dp*(w - nq*.25e0_dp)
1460 :
1461 0 : n = nq/2
1462 0 : IF ((n + n) /= nq) w = 1.0e0_dp - w
1463 0 : z = piov4*w
1464 0 : m = n/2
1465 0 : IF ((m + m) /= n) sgn = -sgn
1466 :
1467 : ! Evaluate -pi*cot(pi*x), using either tan(z) or cot(z) after quadrant reduction.
1468 0 : n = (nq + 1)/2
1469 0 : m = n/2
1470 0 : m = m + m
1471 0 : IF (m /= n) THEN
1472 0 : aug = sgn*((SIN(z)/COS(z))*4.0e0_dp)
1473 : ELSE
1474 0 : IF (z == 0.0e0_dp) THEN
1475 : ! Negative integers are poles of the digamma function.
1476 7564 : fn_val = 0.0e0_dp
1477 : RETURN
1478 : END IF
1479 0 : aug = sgn*((COS(z)/SIN(z))*4.0e0_dp)
1480 : END IF
1481 0 : x = 1.0e0_dp - x
1482 : END IF
1483 : END IF
1484 :
1485 7564 : IF (x <= 3.0e0_dp) THEN
1486 : ! Rational approximation on 0.5 <= x <= 3.
1487 0 : den = x
1488 0 : upper = p1(1)*x
1489 :
1490 0 : DO i = 1, 5
1491 0 : den = (den + q1(i))*x
1492 0 : upper = (upper + p1(i + 1))*x
1493 : END DO
1494 :
1495 0 : den = (upper + p1(7))/(den + q1(6))
1496 0 : xmx0 = x - dx0
1497 0 : fn_val = den*xmx0 + aug
1498 0 : RETURN
1499 : END IF
1500 :
1501 7564 : IF (x < xmax1) THEN
1502 : ! Asymptotic rational correction for 3 < x < xmax1.
1503 7564 : w = 1.0e0_dp/(x*x)
1504 7564 : den = w
1505 7564 : upper = p2(1)*w
1506 :
1507 30256 : DO i = 1, 3
1508 22692 : den = (den + q2(i))*w
1509 30256 : upper = (upper + p2(i + 1))*w
1510 : END DO
1511 :
1512 7564 : aug = upper/(den + q2(4)) - 0.5e0_dp/x + aug
1513 : END IF
1514 : ! For x >= xmax1, the correction is negligible and this reduces to log(x).
1515 7564 : fn_val = aug + LOG(x)
1516 :
1517 7564 : END FUNCTION digamma
1518 :
1519 193589066 : @:define_expint(impure=1)
1520 :
1521 : ! **************************************************************************************************
1522 : !> \brief Jacobi matrix diagonalization. The eigenvalues are returned in
1523 : !> vector d and the eigenvectors are returned in matrix v in ascending
1524 : !> order.
1525 : !>
1526 : !> \param a ...
1527 : !> \param d ...
1528 : !> \param v ...
1529 : !> \par History
1530 : !> - Creation (20.11.98, Matthias Krack)
1531 : ! **************************************************************************************************
1532 37248 : SUBROUTINE jacobi(a, d, v)
1533 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1534 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: d
1535 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: v
1536 :
1537 : INTEGER :: n
1538 :
1539 37248 : n = SIZE(d(:))
1540 :
1541 : ! Diagonalize matrix a
1542 37248 : CALL diag(n, a, d, v)
1543 :
1544 : ! Sort eigenvalues and eigenvector in ascending order
1545 37248 : CALL eigsrt(n, d, v)
1546 :
1547 37248 : END SUBROUTINE jacobi
1548 :
1549 : ! **************************************************************************************************
1550 : !> \brief Diagonalize matrix a. The eigenvalues are returned in vector d
1551 : !> and the eigenvectors are returned in matrix v.
1552 : !>
1553 : !> \param n matrix/vector extent (problem size)
1554 : !> \param a matrix to be diagonalised
1555 : !> \param d vector of eigenvalues
1556 : !> \param v matrix of eigenvectors
1557 : !> \par History
1558 : !> - Creation (20.11.98, Matthias Krack)
1559 : ! **************************************************************************************************
1560 37248 : SUBROUTINE diag(n, a, d, v)
1561 : INTEGER, INTENT(IN) :: n
1562 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1563 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: d
1564 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: v
1565 :
1566 : CHARACTER(len=*), PARAMETER :: routineN = 'diag'
1567 : INTEGER, PARAMETER :: maxit = 100
1568 : REAL(KIND=dp), PARAMETER :: a_eps = 1.0E-10_dp, d_eps = 1.0E-3_dp
1569 :
1570 : CHARACTER(LEN=default_string_length) :: max_iter
1571 : INTEGER :: handle, i, ip, iq
1572 : REAL(KIND=dp) :: a_max, apq, c, d_min, dip, diq, g, h, s, &
1573 : t, tau, theta, tresh
1574 37248 : REAL(KIND=dp), DIMENSION(n) :: b, z
1575 :
1576 37248 : CALL timeset(routineN, handle)
1577 :
1578 37248 : a_max = 0.0_dp
1579 127917 : DO ip = 1, n - 1
1580 890890 : a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip + 1:n))))
1581 127917 : b(ip) = a(ip, ip) ! get_diag(a)
1582 : END DO
1583 37248 : b(n) = a(n, n)
1584 :
1585 37248 : CALL unit_matrix(v)
1586 :
1587 171090 : DO i = 1, maxit
1588 1019528 : d = b
1589 1019528 : d_min = MAX(d_eps, MINVAL(ABS(b)))
1590 171090 : IF (a_max < a_eps*d_min) THEN
1591 37248 : CALL timestop(handle)
1592 : RETURN
1593 : END IF
1594 133842 : tresh = MERGE(a_max, 0.0_dp, (i < 4))
1595 854363 : z = 0.0_dp
1596 720521 : DO ip = 1, n - 1
1597 6480283 : DO iq = ip + 1, n
1598 5759762 : dip = d(ip)
1599 5759762 : diq = d(iq)
1600 5759762 : apq = a(ip, iq)
1601 5759762 : g = 100.0_dp*ABS(apq)
1602 6346441 : IF (tresh < ABS(apq)) THEN
1603 3364536 : h = diq - dip
1604 3364536 : IF ((ABS(h) + g) /= ABS(h)) THEN
1605 2548697 : theta = 0.5_dp*h/apq
1606 2548697 : t = 1.0_dp/(ABS(theta) + SQRT(1.0_dp + theta**2))
1607 2548697 : IF (theta < 0.0_dp) t = -t
1608 : ELSE
1609 815839 : t = apq/h
1610 : END IF
1611 3364536 : c = 1.0_dp/SQRT(1.0_dp + t**2)
1612 3364536 : s = t*c
1613 3364536 : tau = s/(1.0_dp + c)
1614 3364536 : h = t*apq
1615 3364536 : z(ip) = z(ip) - h
1616 3364536 : z(iq) = z(iq) + h
1617 3364536 : d(ip) = dip - h
1618 3364536 : d(iq) = diq + h
1619 3364536 : a(ip, iq) = 0.0_dp
1620 3364536 : CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1621 3364536 : CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1622 3364536 : CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1623 3364536 : CALL jrotate(v(:, ip), v(:, iq), s, tau)
1624 : ELSE IF ((4 < i) .AND. &
1625 0 : ((ABS(dip) + g) == ABS(dip)) .AND. &
1626 2395226 : ((ABS(diq) + g) == ABS(diq))) THEN
1627 0 : a(ip, iq) = 0.0_dp
1628 : END IF
1629 : END DO
1630 : END DO
1631 854363 : b = b + z
1632 : a_max = 0.0_dp
1633 720521 : DO ip = 1, n - 1
1634 6480283 : a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip + 1:n))))
1635 : END DO
1636 : END DO
1637 0 : WRITE (max_iter, *) maxit
1638 : CALL cp_warn(__LOCATION__, &
1639 : "Jacobi matrix diagonalization does not converge after "// &
1640 0 : TRIM(ADJUSTL(max_iter))//" iterations; stop here anyway")
1641 :
1642 0 : CALL timestop(handle)
1643 :
1644 : END SUBROUTINE diag
1645 :
1646 : ! **************************************************************************************************
1647 : !> \brief Perform a Jacobi rotation of the vectors a and b.
1648 : !>
1649 : !> \param a ...
1650 : !> \param b ...
1651 : !> \param ss ...
1652 : !> \param tt ...
1653 : !> \par History
1654 : !> - Creation (20.11.98, Matthias Krack)
1655 : ! **************************************************************************************************
1656 13458144 : PURE SUBROUTINE jrotate(a, b, ss, tt)
1657 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1658 : REAL(KIND=dp), INTENT(IN) :: ss, tt
1659 :
1660 : REAL(KIND=dp) :: u, v
1661 :
1662 13458144 : u = 1.0_dp - ss*tt
1663 13458144 : v = ss/u
1664 :
1665 229055204 : a = a*u - b*ss
1666 229055204 : b = b*(u + ss*v) + a*v
1667 :
1668 13458144 : END SUBROUTINE jrotate
1669 :
1670 : ! **************************************************************************************************
1671 : !> \brief Sort the values in vector d in ascending order and swap the
1672 : !> corresponding columns of matrix v.
1673 : !>
1674 : !> \param n ...
1675 : !> \param d ...
1676 : !> \param v ...
1677 : !> \par History
1678 : !> - Creation (20.11.98, Matthias Krack)
1679 : ! **************************************************************************************************
1680 37248 : SUBROUTINE eigsrt(n, d, v)
1681 : INTEGER, INTENT(IN) :: n
1682 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: d
1683 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: v
1684 :
1685 : INTEGER :: i, j
1686 :
1687 127917 : DO i = 1, n - 1
1688 1072228 : j = SUM(MINLOC(d(i:n))) + i - 1
1689 127917 : IF (j /= i) THEN
1690 43920 : CALL swap(d(i), d(j))
1691 43920 : CALL swap(v(:, i), v(:, j))
1692 : END IF
1693 : END DO
1694 :
1695 37248 : END SUBROUTINE eigsrt
1696 :
1697 : ! **************************************************************************
1698 : !> \brief Swap two scalars
1699 : !>
1700 : !> \param a ...
1701 : !> \param b ...
1702 : !> \par History
1703 : !> - Creation (20.11.98, Matthias Krack)
1704 : ! **************************************************************************************************
1705 43920 : ELEMENTAL SUBROUTINE swap_scalar(a, b)
1706 : REAL(KIND=dp), INTENT(INOUT) :: a, b
1707 :
1708 : REAL(KIND=dp) :: c
1709 :
1710 43920 : c = a
1711 43920 : a = b
1712 43920 : b = c
1713 :
1714 43920 : END SUBROUTINE swap_scalar
1715 :
1716 : ! **************************************************************************
1717 : !> \brief Swap two vectors
1718 : !>
1719 : !> \param a ...
1720 : !> \param b ...
1721 : !> \par History
1722 : !> - Creation (20.11.98, Matthias Krack)
1723 : ! **************************************************************************************************
1724 43920 : SUBROUTINE swap_vector(a, b)
1725 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1726 :
1727 : INTEGER :: i, n
1728 : REAL(KIND=dp) :: c
1729 :
1730 43920 : n = SIZE(a)
1731 :
1732 43920 : IF (n /= SIZE(b)) THEN
1733 0 : CPABORT("Check the array bounds of the parameters")
1734 : END IF
1735 :
1736 962981 : DO i = 1, n
1737 919061 : c = a(i)
1738 919061 : a(i) = b(i)
1739 962981 : b(i) = c
1740 : END DO
1741 :
1742 43920 : END SUBROUTINE swap_vector
1743 :
1744 : ! **************************************************************************************************
1745 : !> \brief - compute a truncation radius for the shortrange operator
1746 : !> \param eps target accuracy!> \param omg screening parameter
1747 : !> \param omg ...
1748 : !> \param r_cutoff cutoff radius
1749 : !> \par History
1750 : !> 10.2012 created [Hossein Banihashemian]
1751 : !> 05.2019 moved here from hfx_types (A. Bussy)
1752 : !> \author Hossein Banihashemian
1753 : ! **************************************************************************************************
1754 66 : SUBROUTINE erfc_cutoff(eps, omg, r_cutoff)
1755 : REAL(dp), INTENT(in) :: eps, omg
1756 : REAL(dp), INTENT(out) :: r_cutoff
1757 :
1758 : CHARACTER(LEN=*), PARAMETER :: routineN = 'erfc_cutoff'
1759 :
1760 : REAL(dp), PARAMETER :: abstol = 1E-10_dp, soltol = 1E-16_dp
1761 : REAL(dp) :: r0, f0, fprime0, delta_r
1762 : INTEGER :: iter, handle
1763 : INTEGER, PARAMETER :: iterMAX = 1000
1764 :
1765 66 : CALL timeset(routineN, handle)
1766 :
1767 : ! initial guess assuming that we are in the asymptotic regime of the erf, and the solution is about 10.
1768 66 : r0 = SQRT(-LOG(eps*omg*10**2))/omg
1769 66 : CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1770 :
1771 638 : DO iter = 1, iterMAX
1772 638 : delta_r = f0/fprime0
1773 638 : r0 = r0 - delta_r
1774 638 : CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1775 638 : IF (ABS(delta_r) < abstol .OR. ABS(f0) < soltol) EXIT
1776 : END DO
1777 66 : CPASSERT(iter <= itermax)
1778 66 : r_cutoff = r0
1779 :
1780 66 : CALL timestop(handle)
1781 : CONTAINS
1782 : ! **************************************************************************************************
1783 : !> \brief ...
1784 : !> \param r ...
1785 : !> \param eps ...
1786 : !> \param omega ...
1787 : !> \param fn ...
1788 : !> \param df ...
1789 : ! **************************************************************************************************
1790 704 : ELEMENTAL SUBROUTINE eval_transc_func(r, eps, omega, fn, df)
1791 : REAL(dp), INTENT(in) :: r, eps, omega
1792 : REAL(dp), INTENT(out) :: fn, df
1793 :
1794 : REAL(dp) :: qr
1795 :
1796 704 : qr = omega*r
1797 704 : fn = ERFC(qr) - r*eps
1798 704 : df = -2.0_dp*oorootpi*omega*EXP(-qr**2) - eps
1799 704 : END SUBROUTINE eval_transc_func
1800 : END SUBROUTINE erfc_cutoff
1801 :
1802 : ! **************************************************************************************************
1803 : !> \brief Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd
1804 : !> \param matrix Hermitian matrix is preserved
1805 : !> \param eigenvectors ...
1806 : !> \param eigenvalues ...
1807 : !> \author A. Bussy
1808 : ! **************************************************************************************************
1809 101013 : SUBROUTINE diag_complex(matrix, eigenvectors, eigenvalues)
1810 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
1811 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
1812 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1813 :
1814 : CHARACTER(len=*), PARAMETER :: routineN = 'diag_complex'
1815 :
1816 101013 : COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1817 : INTEGER :: handle, info, liwork, lrwork, lwork, n
1818 101013 : INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1819 101013 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: rwork
1820 :
1821 101013 : CALL timeset(routineN, handle)
1822 :
1823 101013 : IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) CPABORT("Expected square matrix")
1824 : ! IF (MAXVAL(ABS(matrix - CONJG(TRANSPOSE(matrix)))) > 1e-14_dp) CPABORT("Expected hermitian matrix")
1825 :
1826 101013 : n = SIZE(matrix, 1)
1827 101013 : ALLOCATE (iwork(1), rwork(1), work(1))
1828 :
1829 : ! work space query
1830 101013 : lwork = -1
1831 101013 : lrwork = -1
1832 101013 : liwork = -1
1833 :
1834 101013 : CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1835 :
1836 101013 : lwork = CEILING(REAL(work(1), KIND=dp))
1837 101013 : lrwork = CEILING(rwork(1))
1838 101013 : liwork = iwork(1)
1839 :
1840 101013 : DEALLOCATE (iwork, rwork, work)
1841 707091 : ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
1842 3928883 : eigenvectors(:, :) = matrix(:, :)
1843 :
1844 : ! final diagonalization
1845 101013 : CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1846 :
1847 101013 : DEALLOCATE (iwork, rwork, work)
1848 :
1849 101013 : IF (info /= 0) CPABORT("Diagonalisation of a complex matrix failed")
1850 :
1851 101013 : CALL timestop(handle)
1852 :
1853 101013 : END SUBROUTINE diag_complex
1854 :
1855 : ! **************************************************************************************************
1856 : !> \brief Helper routine for diagonalizing anti symmetric matrices
1857 : !> \param matrix ...
1858 : !> \param evecs ...
1859 : !> \param evals ...
1860 : ! **************************************************************************************************
1861 1049 : SUBROUTINE diag_antisym(matrix, evecs, evals)
1862 : REAL(dp), DIMENSION(:, :) :: matrix
1863 : COMPLEX(dp), DIMENSION(:, :) :: evecs
1864 : COMPLEX(dp), DIMENSION(:) :: evals
1865 :
1866 : COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: matrix_c
1867 : INTEGER :: n
1868 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1869 :
1870 1049 : IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) CPABORT("Expected square matrix")
1871 : ! IF (MAXVAL(ABS(matrix + TRANSPOSE(matrix))) > 1e-14_dp) CPABORT("Expected anti-symmetric matrix")
1872 :
1873 1049 : n = SIZE(matrix, 1)
1874 6294 : ALLOCATE (matrix_c(n, n), eigenvalues(n))
1875 :
1876 32519 : matrix_c(:, :) = CMPLX(0.0_dp, -matrix, kind=dp)
1877 1049 : CALL diag_complex(matrix_c, evecs, eigenvalues)
1878 6294 : evals = CMPLX(0.0_dp, eigenvalues, kind=dp)
1879 :
1880 1049 : DEALLOCATE (matrix_c, eigenvalues)
1881 1049 : END SUBROUTINE diag_antisym
1882 : ! **************************************************************************************************
1883 : !> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged
1884 : !> \param A_in Input matrix 1
1885 : !> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1886 : !> \param B_in Input matrix 2
1887 : !> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1888 : !> \param C_out Output matrix
1889 : !> \par History
1890 : !> 11.2025 created [Stepan Marek]
1891 : ! **************************************************************************************************
1892 88014 : SUBROUTINE zgemm_square_2(A_in, A_trans, B_in, B_trans, C_out)
1893 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1894 : CHARACTER, INTENT(IN) :: A_trans
1895 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1896 : CHARACTER, INTENT(IN) :: B_trans
1897 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: C_out
1898 :
1899 : CHARACTER(len=*), PARAMETER :: routineN = 'zgemm_square_2'
1900 :
1901 : INTEGER :: handle, n
1902 :
1903 88014 : n = SIZE(A_in, 1)
1904 88014 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1905 88014 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1906 88014 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1907 88014 : IF (n /= SIZE(C_out, 1)) CPABORT("Incompatible (rows) result array 3 (C).")
1908 88014 : IF (n /= SIZE(C_out, 2)) CPABORT("Incompatible (cols) result array 3 (C).")
1909 88014 : IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
1910 : A_trans == 'T' .OR. A_trans == 't' .OR. &
1911 : A_trans == 'C' .OR. A_trans == 'c')) THEN
1912 0 : CPABORT("Unknown transpose character for array 1 (A).")
1913 : END IF
1914 88014 : IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
1915 : B_trans == 'T' .OR. B_trans == 't' .OR. &
1916 : B_trans == 'C' .OR. B_trans == 'c')) THEN
1917 0 : CPABORT("Unknown transpose character for array 2 (B).")
1918 : END IF
1919 :
1920 88014 : CALL timeset(routineN, handle)
1921 :
1922 88014 : CALL ZGEMM(A_trans, B_trans, n, n, n, z_one, A_in, n, B_in, n, z_zero, C_out, n)
1923 :
1924 88014 : CALL timestop(handle)
1925 :
1926 88014 : END SUBROUTINE zgemm_square_2
1927 : ! **************************************************************************************************
1928 : !> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, real matrices
1929 : !> \param A_in Input matrix 1
1930 : !> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1931 : !> \param B_in Input matrix 2
1932 : !> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1933 : !> \param C_out Output matrix
1934 : !> \par History
1935 : !> 11.2025 created [Stepan Marek]
1936 : ! **************************************************************************************************
1937 2 : SUBROUTINE dgemm_square_2(A_in, A_trans, B_in, B_trans, C_out)
1938 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1939 : CHARACTER, INTENT(IN) :: A_trans
1940 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1941 : CHARACTER, INTENT(IN) :: B_trans
1942 : REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: C_out
1943 :
1944 : CHARACTER(len=*), PARAMETER :: routineN = 'dgemm_square_2'
1945 :
1946 : INTEGER :: handle, n
1947 :
1948 2 : n = SIZE(A_in, 1)
1949 2 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1950 2 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1951 2 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1952 2 : IF (n /= SIZE(C_out, 1)) CPABORT("Incompatible (rows) result array 3 (C).")
1953 2 : IF (n /= SIZE(C_out, 2)) CPABORT("Incompatible (cols) result array 3 (C).")
1954 2 : IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
1955 : A_trans == 'T' .OR. A_trans == 't' .OR. &
1956 : A_trans == 'C' .OR. A_trans == 'c')) THEN
1957 0 : CPABORT("Unknown transpose character for array 1 (A).")
1958 : END IF
1959 2 : IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
1960 : B_trans == 'T' .OR. B_trans == 't' .OR. &
1961 : B_trans == 'C' .OR. B_trans == 'c')) THEN
1962 0 : CPABORT("Unknown transpose character for array 2 (B).")
1963 : END IF
1964 :
1965 2 : CALL timeset(routineN, handle)
1966 :
1967 2 : CALL DGEMM(A_trans, B_trans, n, n, n, 1.0_dp, A_in, n, B_in, n, 0.0_dp, C_out, n)
1968 :
1969 2 : CALL timestop(handle)
1970 :
1971 2 : END SUBROUTINE dgemm_square_2
1972 : ! **************************************************************************************************
1973 : !> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, for 3 matrices
1974 : !> \param A_in Input matrix 1
1975 : !> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1976 : !> \param B_in Input matrix 2
1977 : !> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1978 : !> \param C_in Input matrix 3
1979 : !> \param C_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 3
1980 : !> \param D_out Output matrix
1981 : !> \par History
1982 : !> 11.2025 created [Stepan Marek]
1983 : ! **************************************************************************************************
1984 68137 : SUBROUTINE zgemm_square_3(A_in, A_trans, B_in, B_trans, C_in, C_trans, D_out)
1985 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1986 : CHARACTER, INTENT(IN) :: A_trans
1987 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1988 : CHARACTER, INTENT(IN) :: B_trans
1989 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: C_in
1990 : CHARACTER, INTENT(IN) :: C_trans
1991 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: D_out
1992 :
1993 : CHARACTER(len=*), PARAMETER :: routineN = 'zgemm_square_3'
1994 :
1995 68137 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1996 : INTEGER :: handle, n
1997 :
1998 68137 : n = SIZE(A_in, 1)
1999 68137 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
2000 68137 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
2001 68137 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
2002 68137 : IF (n /= SIZE(C_in, 1)) CPABORT("Incompatible (rows) array 3 (C).")
2003 68137 : IF (n /= SIZE(C_in, 2)) CPABORT("Non-square array 3 (C).")
2004 68137 : IF (n /= SIZE(D_out, 1)) CPABORT("Incompatible (rows) result array 4 (D).")
2005 68137 : IF (n /= SIZE(D_out, 2)) CPABORT("Incompatible (cols) result array 4 (D).")
2006 68137 : IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
2007 : A_trans == 'T' .OR. A_trans == 't' .OR. &
2008 : A_trans == 'C' .OR. A_trans == 'c')) THEN
2009 0 : CPABORT("Unknown transpose character for array 1 (A).")
2010 : END IF
2011 68137 : IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
2012 : B_trans == 'T' .OR. B_trans == 't' .OR. &
2013 : B_trans == 'C' .OR. B_trans == 'c')) THEN
2014 0 : CPABORT("Unknown transpose character for array 2 (B).")
2015 : END IF
2016 68137 : IF (.NOT. (C_trans == 'N' .OR. C_trans == 'n' .OR. &
2017 : C_trans == 'T' .OR. C_trans == 't' .OR. &
2018 : C_trans == 'C' .OR. C_trans == 'c')) THEN
2019 0 : CPABORT("Unknown transpose character for array 3 (C).")
2020 : END IF
2021 :
2022 68137 : CALL timeset(routineN, handle)
2023 :
2024 272548 : ALLOCATE (work(n, n), source=z_zero)
2025 :
2026 68137 : CALL ZGEMM(A_trans, B_trans, n, n, n, z_one, A_in, n, B_in, n, z_zero, work, n)
2027 68137 : CALL ZGEMM('N', C_trans, n, n, n, z_one, work, n, C_in, n, z_zero, D_out, n)
2028 :
2029 68137 : DEALLOCATE (work)
2030 :
2031 68137 : CALL timestop(handle)
2032 :
2033 68137 : END SUBROUTINE zgemm_square_3
2034 : ! **************************************************************************************************
2035 : !> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, for 3 matrices
2036 : !> \param A_in Input matrix 1
2037 : !> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
2038 : !> \param B_in Input matrix 2
2039 : !> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
2040 : !> \param C_in Input matrix 3
2041 : !> \param C_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 3
2042 : !> \param D_out Output matrix
2043 : !> \par History
2044 : !> 11.2025 created [Stepan Marek]
2045 : ! **************************************************************************************************
2046 4 : SUBROUTINE dgemm_square_3(A_in, A_trans, B_in, B_trans, C_in, C_trans, D_out)
2047 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
2048 : CHARACTER, INTENT(IN) :: A_trans
2049 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
2050 : CHARACTER, INTENT(IN) :: B_trans
2051 : REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: C_in
2052 : CHARACTER, INTENT(IN) :: C_trans
2053 : REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: D_out
2054 :
2055 : CHARACTER(len=*), PARAMETER :: routineN = 'dgemm_square_3'
2056 :
2057 : INTEGER :: handle, n
2058 4 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
2059 :
2060 4 : n = SIZE(A_in, 1)
2061 4 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
2062 4 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
2063 4 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
2064 4 : IF (n /= SIZE(C_in, 1)) CPABORT("Incompatible (rows) array 3 (C).")
2065 4 : IF (n /= SIZE(C_in, 2)) CPABORT("Non-square array 3 (C).")
2066 4 : IF (n /= SIZE(D_out, 1)) CPABORT("Incompatible (rows) result array 4 (D).")
2067 4 : IF (n /= SIZE(D_out, 2)) CPABORT("Incompatible (cols) result array 4 (D).")
2068 4 : IF (.NOT. (A_trans == 'N' .OR. A_trans == 'n' .OR. &
2069 : A_trans == 'T' .OR. A_trans == 't' .OR. &
2070 : A_trans == 'C' .OR. A_trans == 'c')) THEN
2071 0 : CPABORT("Unknown transpose character for array 1 (A).")
2072 : END IF
2073 4 : IF (.NOT. (B_trans == 'N' .OR. B_trans == 'n' .OR. &
2074 : B_trans == 'T' .OR. B_trans == 't' .OR. &
2075 : B_trans == 'C' .OR. B_trans == 'c')) THEN
2076 0 : CPABORT("Unknown transpose character for array 2 (B).")
2077 : END IF
2078 4 : IF (.NOT. (C_trans == 'N' .OR. C_trans == 'n' .OR. &
2079 : C_trans == 'T' .OR. C_trans == 't' .OR. &
2080 : C_trans == 'C' .OR. C_trans == 'c')) THEN
2081 0 : CPABORT("Unknown transpose character for array 3 (C).")
2082 : END IF
2083 :
2084 4 : CALL timeset(routineN, handle)
2085 :
2086 16 : ALLOCATE (work(n, n), source=0.0_dp)
2087 :
2088 4 : CALL DGEMM(A_trans, B_trans, n, n, n, 1.0_dp, A_in, n, B_in, n, 0.0_dp, work, n)
2089 4 : CALL DGEMM('N', C_trans, n, n, n, 1.0_dp, work, n, C_in, n, 0.0_dp, D_out, n)
2090 :
2091 4 : DEALLOCATE (work)
2092 :
2093 4 : CALL timestop(handle)
2094 :
2095 4 : END SUBROUTINE dgemm_square_3
2096 :
2097 : ! **************************************************************************************************
2098 : !> \brief Solve the generalized eigenvalue equation for complex matrices A*v = B*v*λ
2099 : !> \param A_in ...
2100 : !> \param B_in ...
2101 : !> \param eigenvalues ...
2102 : !> \param eigenvectors ...
2103 : !> \author Shridhar Shanbhag
2104 : ! **************************************************************************************************
2105 6 : SUBROUTINE geeig_right(A_in, B_in, eigenvalues, eigenvectors)
2106 : COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in, B_in
2107 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
2108 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
2109 :
2110 : CHARACTER(len=*), PARAMETER :: routineN = 'geeig_right'
2111 : COMPLEX(KIND=dp), PARAMETER :: cone = CMPLX(1.0_dp, 0.0_dp, KIND=dp), &
2112 : czero = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
2113 :
2114 6 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
2115 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: A, B, work
2116 : INTEGER :: handle, i, icol, irow, lda, ldb, ldc, &
2117 : nao, nc, ncol, nmo, nx
2118 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
2119 :
2120 6 : CALL timeset(routineN, handle)
2121 :
2122 : ! Test sizes
2123 6 : nao = SIZE(A_in, 1)
2124 6 : nmo = SIZE(eigenvalues)
2125 30 : ALLOCATE (evals(nao), cevals(nao))
2126 48 : ALLOCATE (work(nao, nao), B(nao, nao), A(nao, nao))
2127 22334 : A(:, :) = A_in(:, :)
2128 22334 : B(:, :) = B_in(:, :)
2129 :
2130 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
2131 22334 : B = -B
2132 6 : CALL diag_complex(B, work, evals)
2133 254 : evals(:) = -evals(:)
2134 6 : nc = nao
2135 254 : DO i = 1, nao
2136 254 : IF (evals(i) < -1.0_dp) THEN
2137 0 : nc = i - 1
2138 0 : EXIT
2139 : END IF
2140 : END DO
2141 6 : CPASSERT(nc /= 0)
2142 :
2143 6 : IF (nc /= nao) THEN
2144 0 : IF (nc < nmo) THEN
2145 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
2146 0 : ncol = nmo - nc
2147 0 : CALL zcopy(ncol*nao, work(1, nc + 1), 1, eigenvectors(1, nc + 1), 1)
2148 : END IF
2149 : ! Set NULL space in eigenvector matrix of S to zero
2150 0 : DO icol = nc + 1, nao
2151 0 : DO irow = 1, nao
2152 0 : work(irow, icol) = czero
2153 : END DO
2154 : END DO
2155 : ! Set small eigenvalues to a dummy save value
2156 0 : evals(nc + 1:nao) = 1.0_dp
2157 : END IF
2158 : ! Calculate U*s**(-1/2)
2159 254 : cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
2160 254 : DO i = 1, MIN(SIZE(work, 2), SIZE(cevals))
2161 22334 : work(1:MIN(SIZE(work, 2), SIZE(cevals)), i) = cevals(i)*work(1:MIN(SIZE(work, 2), SIZE(cevals)), i)
2162 : END DO
2163 : ! Reduce to get U^(-C) * H * U^(-1)
2164 6 : CALL gemm_square(work, 'C', A, 'N', B)
2165 6 : CALL gemm_square(B, 'N', work, 'N', A)
2166 6 : IF (nc /= nao) THEN
2167 : ! set diagonal values to save large value
2168 0 : DO icol = nc + 1, nao
2169 0 : A(icol, icol) = 10000*cone
2170 : END DO
2171 : END IF
2172 : ! Diagonalize
2173 6 : CALL diag_complex(A, B, evals)
2174 254 : eigenvalues(1:nmo) = evals(1:nmo)
2175 6 : nx = MIN(nc, nmo)
2176 : ! Restore vectors C = U^(-1) * C*
2177 6 : lda = SIZE(work, 1)
2178 6 : ldb = SIZE(B, 1)
2179 6 : ldc = SIZE(eigenvectors, 1)
2180 : CALL ZGEMM("N", "N", nao, nx, nc, cone, work, &
2181 6 : lda, B, ldb, czero, eigenvectors, ldc)
2182 :
2183 6 : DEALLOCATE (evals, cevals, work, B, A)
2184 :
2185 6 : CALL timestop(handle)
2186 6 : END SUBROUTINE geeig_right
2187 :
2188 : END MODULE mathlib
|
