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