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