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