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