(git:d18deda)
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-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! **************************************************************************************************
15MODULE 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
86CONTAINS
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 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 cpassert(b > a)
107 IF (x < a .OR. x > b) THEN
108 ! outside switching intervall
109 IF (order > 0) THEN
110 ! derivatives are 0
111 fx = 0.0_dp
112 ELSE
113 IF (x < a) THEN
114 ! x < a => f(x) = 1
115 fx = 1.0_dp
116 ELSE
117 ! x > b => f(x) = 0
118 fx = 0.0_dp
119 END IF
120 END IF
121 ELSE
122 ! renormalized coordinate
123 u = (x - a)/(b - a)
124 SELECT CASE (order)
125 CASE (0)
126 u2 = u*u
127 u3 = u2*u
128 fx = 1._dp - 10._dp*u3 + 15._dp*u2*u2 - 6._dp*u2*u3
129 CASE (1)
130 u2 = u*u
131 fx = -30._dp*u2 + 60._dp*u*u2 - 30._dp*u2*u2
132 fx = fx/(b - a)
133 CASE (2)
134 u2 = u*u
135 fx = -60._dp*u + 180._dp*u2 - 120._dp*u*u2
136 fx = fx/(b - a)**2
137 CASE DEFAULT
138 cpabort('order not defined')
139 END SELECT
140 END IF
141
142 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 LOGICAL FUNCTION abnormal_value(a)
151 REAL(kind=dp) :: a
152
153 CHARACTER(LEN=32) :: buffer
154
155 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 WRITE (buffer, *) a
161 IF (index(buffer, "N") .NE. 0 .OR. index(buffer, "n") .NE. 0) abnormal_value = .true.
162
163 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 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 REAL(kind=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
181
182 length_of_a = sqrt(dot_product(a, a))
183 length_of_b = sqrt(dot_product(b, b))
184
185 IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
186 a_norm(:) = a(:)/length_of_a
187 b_norm(:) = b(:)/length_of_b
188 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 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 ELEMENTAL FUNCTION binomial(n, k) RESULT(n_over_k)
206 INTEGER, INTENT(IN) :: n, k
207 REAL(kind=dp) :: n_over_k
208
209 IF ((k >= 0) .AND. (k <= n)) THEN
210 n_over_k = fac(n)/(fac(n - k)*fac(k))
211 ELSE
212 n_over_k = 0.0_dp
213 END IF
214
215 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 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 IF (k >= 0) THEN
237 z_over_k = 1.0_dp
238 DO i = 1, k
239 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 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 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 IF (all(k >= 0) .AND. sum(k) == n) THEN
263 denom = 1.0_dp
264 DO i = 1, SIZE(k)
265 denom = denom*fac(k(i))
266 END DO
267 res = fac(n)/denom
268 ELSE
269 res = 0.0_dp
270 END IF
271
272 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 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 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 IF (length_of_a > eps_geo) THEN
296
297 d(:) = a(:)/length_of_a
298
299 cosp = cos(phi)
300 sinp = sin(phi)
301 cost = 1.0_dp - cosp
302
303 rotmat(1, 1) = d(1)*d(1)*cost + cosp
304 rotmat(1, 2) = d(1)*d(2)*cost - d(3)*sinp
305 rotmat(1, 3) = d(1)*d(3)*cost + d(2)*sinp
306 rotmat(2, 1) = d(2)*d(1)*cost + d(3)*sinp
307 rotmat(2, 2) = d(2)*d(2)*cost + cosp
308 rotmat(2, 3) = d(2)*d(3)*cost - d(1)*sinp
309 rotmat(3, 1) = d(3)*d(1)*cost - d(2)*sinp
310 rotmat(3, 2) = d(3)*d(2)*cost + d(1)*sinp
311 rotmat(3, 3) = d(3)*d(3)*cost + cosp
312 ELSE
313 CALL unit_matrix(rotmat)
314 END IF
315
316 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 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 a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
333
334 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 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 a3(1)*(a1(2)*a2(3) - a2(2)*a1(3))
353
354 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
381 INTEGER, EXTERNAL :: ilaenv
382 LOGICAL :: divide_and_conquer
383 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
384
385 EXTERNAL dsyev, dsyevd
386
387 CALL timeset(routinen, handle)
388
389 ! Get the size of the matrix a
390 n = SIZE(a, 1)
391
392 ! Check the size of matrix a
393 IF (SIZE(a, 2) /= n) THEN
394 cpabort("Check the size of matrix a (parameter #1)")
395 END IF
396
397 ! Check the size of vector eigval
398 IF (SIZE(eigval) /= n) THEN
399 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 IF (PRESENT(dac)) THEN
405 divide_and_conquer = dac
406 ELSE
407 divide_and_conquer = .false.
408 END IF
409
410 ! Get the optimal work storage size
411
412 IF (divide_and_conquer) THEN
413 lwork = 2*n**2 + 6*n + 1
414 liwork = 5*n + 3
415 ELSE
416 nb = ilaenv(1, "DSYTRD", "U", n, -1, -1, -1)
417 lwork = (nb + 2)*n
418 END IF
419
420 ! Allocate work storage
421
422 ALLOCATE (work(lwork))
423 IF (divide_and_conquer) THEN
424 ALLOCATE (iwork(liwork))
425 END IF
426
427 ! Diagonalize the matrix a
428
429 info = 0
430 IF (divide_and_conquer) THEN
431 CALL dsyevd("V", "U", n, a, n, eigval, work, lwork, iwork, liwork, info)
432 ELSE
433 CALL dsyev("V", "U", n, a, n, eigval, work, lwork, info)
434 END IF
435
436 IF (info /= 0) THEN
437 IF (divide_and_conquer) THEN
438 cpabort("The matrix diagonalization with dsyevd failed")
439 ELSE
440 cpabort("The matrix diagonalization with dsyev failed")
441 END IF
442 END IF
443
444 ! Release work storage
445 DEALLOCATE (work)
446
447 IF (divide_and_conquer) THEN
448 DEALLOCATE (iwork)
449 END IF
450
451 CALL timestop(handle)
452
453 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 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 abc = vector_product(bc, -ab)
475 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 det_abcd = det_3x3(abc, bcd, -bc)
480 dihedral_angle_abcd = sign(1.0_dp, det_abcd)*angle(abc, bcd)
481
482 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 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 n = min(SIZE(a, 1), SIZE(a, 2))
500
501 DO i = 1, n
502 a_diag(i) = a(i, i)
503 END DO
504
505 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 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 det_a = 1.0_dp/det_3x3(a)
522
523 a_inv(1, 1) = (a(2, 2)*a(3, 3) - a(3, 2)*a(2, 3))*det_a
524 a_inv(2, 1) = (a(2, 3)*a(3, 1) - a(3, 3)*a(2, 1))*det_a
525 a_inv(3, 1) = (a(2, 1)*a(3, 2) - a(3, 1)*a(2, 2))*det_a
526
527 a_inv(1, 2) = (a(1, 3)*a(3, 2) - a(3, 3)*a(1, 2))*det_a
528 a_inv(2, 2) = (a(1, 1)*a(3, 3) - a(3, 1)*a(1, 3))*det_a
529 a_inv(3, 2) = (a(1, 2)*a(3, 1) - a(3, 2)*a(1, 1))*det_a
530
531 a_inv(1, 3) = (a(1, 2)*a(2, 3) - a(2, 2)*a(1, 3))*det_a
532 a_inv(2, 3) = (a(1, 3)*a(2, 1) - a(2, 3)*a(1, 1))*det_a
533 a_inv(3, 3) = (a(1, 1)*a(2, 2) - a(2, 1)*a(1, 2))*det_a
534
535 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 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 INTEGER, ALLOCATABLE :: ipiv(:)
550 REAL(kind=dp), ALLOCATABLE :: work(:)
551
552 CALL timeset(routinen, handle)
553
554 n = SIZE(a, 1)
555 lwork = 20*n
556 ALLOCATE (ipiv(n))
557 ALLOCATE (work(lwork))
558 ipiv = 0
559 work = 0._dp
560 info = 0
561 CALL dgetrf(n, n, a, n, ipiv, info)
562 IF (info == 0) THEN
563 CALL dgetri(n, a, n, ipiv, work, lwork, info)
564 END IF
565 DEALLOCATE (ipiv, work)
566
567 CALL timestop(handle)
568
569 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 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 CALL timeset(routinen, handle)
591
592 do_potrf = .true.
593 IF (PRESENT(potrf)) do_potrf = potrf
594
595 myuplo = 'U'
596 IF (PRESENT(uplo)) myuplo = uplo
597
598 n = SIZE(a, 1)
599 info = 0
600
601 ! do cholesky decomposition
602 IF (do_potrf) THEN
603 CALL dpotrf(myuplo, n, a, n, info)
604 IF (info /= 0) cpabort("DPOTRF failed")
605 END IF
606
607 ! do inversion using the cholesky decomposition
608 CALL dpotri(myuplo, n, a, n, info)
609 IF (info /= 0) cpabort("Matrix inversion failed")
610
611 ! complete the matrix
612 IF ((myuplo == "U") .OR. (myuplo == "u")) THEN
613 CALL symmetrize_matrix(a, "upper_to_lower")
614 ELSE
615 CALL symmetrize_matrix(a, "lower_to_upper")
616 END IF
617
618 CALL timestop(handle)
619
620 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv, iwork
655 LOGICAL :: do_improve
656 REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
657 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, work
658 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 IF (PRESENT(option)) THEN
665 trans = option
666 ELSE
667 trans = "N"
668 END IF
669
670 IF (PRESENT(improve)) THEN
671 do_improve = improve
672 ELSE
673 do_improve = .true.
674 END IF
675
676 ! Get the dimension of matrix a
677 n = SIZE(a, 1)
678
679 ! Check array dimensions
680 IF (n == 0) THEN
681 cpabort("Matrix to be inverted of zero size")
682 END IF
683
684 IF (n /= SIZE(a, 2)) THEN
685 cpabort("Check the array bounds of parameter #1")
686 END IF
687
688 IF ((n /= SIZE(a_inverse, 1)) .OR. &
689 (n /= SIZE(a_inverse, 2))) THEN
690 cpabort("Check the array bounds of parameter #2")
691 END IF
692
693 ! Allocate work storage
694 ALLOCATE (a_lu(n, n))
695 ALLOCATE (b(n, n))
696 ALLOCATE (berr(n))
697 ALLOCATE (ferr(n))
698 ALLOCATE (ipiv(n))
699 ALLOCATE (iwork(n))
700 ALLOCATE (work(4*n))
701
702 a_lu(1:n, 1:n) = a(1:n, 1:n)
703
704 ! Compute the LU factorization of the matrix a
705 CALL dgetrf(n, n, a_lu, n, ipiv, info)
706
707 IF (info /= 0) THEN
708 cpabort("The LU factorization in dgetrf failed")
709 END IF
710
711 ! Compute the norm of the matrix a
712
713 IF (trans == "N") THEN
714 norm = '1'
715 ELSE
716 norm = 'I'
717 END IF
718
719 a_norm = dlange(norm, n, n, a, n, work)
720
721 ! Compute the reciprocal of the condition number of a
722
723 CALL dgecon(norm, n, a_lu, n, a_norm, r_cond, work, iwork, info)
724
725 IF (info /= 0) THEN
726 cpabort("The computation of the condition number in dgecon failed")
727 END IF
728
729 IF (r_cond < epsilon(0.0_dp)) THEN
730 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 "working precision)")
734 END IF
735
736 ! Solve a system of linear equations using the LU factorization computed by dgetrf
737
738 CALL unit_matrix(a_inverse)
739
740 CALL dgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
741
742 IF (info /= 0) THEN
743 cpabort("Solving the system of linear equations in dgetrs failed")
744 END IF
745
746 ! Improve the computed solution iteratively
747 CALL unit_matrix(b) ! Initialize right-hand sides
748
749 eval_error = 0.0_dp
750
751 IF (do_improve) THEN
752 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 work, iwork, info)
756
757 IF (info /= 0) THEN
758 cpabort("Improving the computed solution in dgerfs failed")
759 END IF
760
761 old_eval_error = eval_error
762 eval_error = maxval(ferr)
763
764 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 DEALLOCATE (work)
771 DEALLOCATE (iwork)
772 DEALLOCATE (ipiv)
773 DEALLOCATE (ferr)
774 DEALLOCATE (berr)
775 DEALLOCATE (b)
776 DEALLOCATE (a_lu)
777
778 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 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 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
809 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
810 INTEGER :: info, iter, n
811 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
812 REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
813 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 IF (PRESENT(option)) THEN
820 trans = option
821 ELSE
822 trans = "N"
823 END IF
824
825 ! Get the dimension of matrix a
826 n = SIZE(a, 1)
827
828 ! Check array dimensions
829 IF (n == 0) THEN
830 cpabort("Matrix to be inverted of zero size")
831 END IF
832
833 IF (n /= SIZE(a, 2)) THEN
834 cpabort("Check the array bounds of parameter #1")
835 END IF
836
837 IF ((n /= SIZE(a_inverse, 1)) .OR. &
838 (n /= SIZE(a_inverse, 2))) THEN
839 cpabort("Check the array bounds of parameter #2")
840 END IF
841
842 ! Allocate work storage
843 ALLOCATE (a_lu(n, n))
844 ALLOCATE (b(n, n))
845 ALLOCATE (berr(n))
846 ALLOCATE (ferr(n))
847 ALLOCATE (ipiv(n))
848 ALLOCATE (rwork(2*n))
849 ALLOCATE (work(2*n))
850
851 a_lu(1:n, 1:n) = a(1:n, 1:n)
852
853 ! Compute the LU factorization of the matrix a
854 CALL zgetrf(n, n, a_lu, n, ipiv, info)
855
856 IF (info /= 0) THEN
857 cpabort("The LU factorization in dgetrf failed")
858 END IF
859
860 ! Compute the norm of the matrix a
861
862 IF (trans == "N") THEN
863 norm = '1'
864 ELSE
865 norm = 'I'
866 END IF
867
868 a_norm = zlange(norm, n, n, a, n, work)
869
870 ! Compute the reciprocal of the condition number of a
871
872 CALL zgecon(norm, n, a_lu, n, a_norm, r_cond, work, rwork, info)
873
874 IF (info /= 0) THEN
875 cpabort("The computation of the condition number in dgecon failed")
876 END IF
877
878 IF (r_cond < epsilon(0.0_dp)) THEN
879 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 "working precision)")
883 END IF
884
885 ! Solve a system of linear equations using the LU factorization computed by dgetrf
886
887 CALL unit_matrix(a_inverse)
888
889 CALL zgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
890
891 IF (info /= 0) THEN
892 cpabort("Solving the system of linear equations in dgetrs failed")
893 END IF
894
895 ! Improve the computed solution iteratively
896 CALL unit_matrix(b) ! Initialize right-hand sides
897
898 eval_error = 0.0_dp
899
900 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 work, rwork, info)
904
905 IF (info /= 0) THEN
906 cpabort("Improving the computed solution in dgerfs failed")
907 END IF
908
909 old_eval_error = eval_error
910 eval_error = maxval(ferr)
911
912 IF (abs(eval_error - old_eval_error) <= epsilon(1.0_dp)) EXIT
913
914 END DO
915
916 ! Release work storage
917 DEALLOCATE (work)
918 DEALLOCATE (rwork)
919 DEALLOCATE (ipiv)
920 DEALLOCATE (ferr)
921 DEALLOCATE (berr)
922 DEALLOCATE (b)
923 DEALLOCATE (a_lu)
924
925 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
948 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sig, work
949 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sig_plus, temp_mat, u, vt
950
951 CALL timeset(routinen, handle)
952
953 n = SIZE(a, 1)
954 ALLOCATE (u(n, n), vt(n, n), sig(n), sig_plus(n, n), iwork(8*n), work(1), temp_mat(n, n))
955 u(:, :) = 0.0_dp
956 vt(:, :) = 0.0_dp
957 sig(:) = 0.0_dp
958 sig_plus = 0.0_dp
959 work = 0.0_dp
960 iwork = 0
961 IF (PRESENT(determinant)) determinant = 1.0_dp
962
963 ! work size query
964 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 lwork, iwork(1), info)
967
968 IF (info /= 0) THEN
969 cpabort("ERROR in DGESDD: Could not retrieve work array sizes")
970 END IF
971 lwork = int(work(1))
972 DEALLOCATE (work)
973 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 lwork, iwork(1), info)
978
979 IF (info /= 0) THEN
980 cpabort("SVD failed")
981 END IF
982
983 IF (PRESENT(sval)) THEN
984 cpassert(.NOT. ASSOCIATED(sval))
985 ALLOCATE (sval(n))
986 sval(:) = sig
987 END IF
988
989 ! set singular values that are too small to zero
990 DO i = 1, n
991 IF (sig(i) > rskip*maxval(sig)) THEN
992 IF (PRESENT(determinant)) &
993 determinant = determinant*sig(i)
994 sig_plus(i, i) = 1._dp/sig(i)
995 ELSE
996 sig_plus(i, i) = 0.0_dp
997 END IF
998 END DO
999
1000 ! build pseudoinverse: V*sig_plus*UT
1001 CALL dgemm("N", "T", n, n, n, 1._dp, sig_plus, n, u, n, 0._dp, temp_mat, n)
1002 CALL dgemm("T", "N", n, n, n, 1._dp, vt, n, temp_mat, n, 0._dp, a_pinverse, n)
1003
1004 DEALLOCATE (u, vt, sig, iwork, work, sig_plus, temp_mat)
1005
1006 CALL timestop(handle)
1007
1008 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
1026 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dinv, temp_mat
1027
1028 CALL timeset(routinen, handle)
1029
1030 info = 0
1031 n = SIZE(a, 1)
1032 ALLOCATE (dinv(n, n), eig(n), work(1), temp_mat(n, n))
1033 dinv(:, :) = 0.0_dp
1034 eig(:) = 0.0_dp
1035 work(:) = 0.0_dp
1036 temp_mat = 0.0_dp
1037
1038 ! work size query
1039 lwork = -1
1040 CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1041 IF (info /= 0) THEN
1042 cpabort("ERROR in DSYEV: Could not retrieve work array sizes")
1043 END IF
1044 lwork = int(work(1))
1045 DEALLOCATE (work)
1046 ALLOCATE (work(lwork))
1047 work = 0.0_dp
1048
1049 ! get eigenvalues and eigenvectors
1050 CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1051
1052 IF (info /= 0) THEN
1053 cpabort("Matrix diagonalization failed")
1054 END IF
1055
1056 ! set eigenvalues that are too small to zero
1057 DO i = 1, n
1058 IF (eig(i) > rskip*maxval(eig)) THEN
1059 dinv(i, i) = 1.0_dp/eig(i)
1060 ELSE
1061 dinv(i, i) = 0._dp
1062 END IF
1063 END DO
1064
1065 ! build pseudoinverse: U*dinv*UT
1066 CALL dgemm("N", "T", n, n, n, 1._dp, dinv, n, a, n, 0._dp, temp_mat, n)
1067 CALL dgemm("N", "N", n, n, n, 1._dp, a, n, temp_mat, n, 0._dp, a_pinverse, n)
1068
1069 DEALLOCATE (eig, work, dinv, temp_mat)
1070
1071 CALL timestop(handle)
1072
1073 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 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 length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1093
1094 IF (length_of_b > eps_geo) THEN
1095
1096 d(:) = b(:)/length_of_b
1097
1098 ! Calculate the mirror image a_mirror of the vector a
1099 scapro = a(1)*d(1) + a(2)*d(2) + a(3)*d(3)
1100
1101 a_mirror(:) = a(:) - 2.0_dp*scapro*d(:)
1102
1103 ELSE
1104
1105 a_mirror(:) = 0.0_dp
1106
1107 END IF
1108
1109 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 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 length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1133 IF (length_of_b > eps_geo) THEN
1134
1135 ! Build up the rotation matrix rotmat
1136 CALL build_rotmat(phi, b, rotmat)
1137
1138 ! Rotate the vector a by phi about the axis defined by vector b
1139 a_rot(:) = matmul(rotmat, a)
1140
1141 ELSE
1142
1143 a_rot(:) = 0.0_dp
1144
1145 END IF
1146
1147 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 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 n = min(SIZE(a, 1), SIZE(a, 2))
1164 DO i = 1, n
1165 a(i, i) = b
1166 END DO
1167
1168 END SUBROUTINE set_diag_scalar_d
1169
1170! **************************************************************************************************
1171!> \brief ...
1172!> \param a ...
1173!> \param b ...
1174! **************************************************************************************************
1175 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 n = min(SIZE(a, 1), SIZE(a, 2))
1182 DO i = 1, n
1183 a(i, i) = b
1184 END DO
1185
1186 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 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 n = min(SIZE(a, 1), SIZE(a, 2))
1203
1204 IF (option == "lower_to_upper") THEN
1205 DO i = 1, n - 1
1206 a(i, i + 1:n) = a(i + 1:n, i)
1207 END DO
1208 ELSE IF (option == "upper_to_lower") THEN
1209 DO i = 1, n - 1
1210 a(i + 1:n, i) = a(i, i + 1:n)
1211 END DO
1212 ELSE IF (option == "anti_lower_to_upper") THEN
1213 DO i = 1, n - 1
1214 a(i, i + 1:n) = -a(i + 1:n, i)
1215 END DO
1216 ELSE IF (option == "anti_upper_to_lower") THEN
1217 DO i = 1, n - 1
1218 a(i + 1:n, i) = -a(i, i + 1:n)
1219 END DO
1220 ELSE
1221 cpabort("Invalid option <"//trim(option)//"> was specified for parameter #2")
1222 END IF
1223
1224 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 PURE SUBROUTINE unit_matrix_d(a)
1234 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1235
1236 a(:, :) = 0.0_dp
1237 CALL set_diag(a, 1.0_dp)
1238
1239 END SUBROUTINE unit_matrix_d
1240
1241! **************************************************************************************************
1242!> \brief ...
1243!> \param a ...
1244! **************************************************************************************************
1245 PURE SUBROUTINE unit_matrix_z(a)
1246 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1247
1248 a(:, :) = (0.0_dp, 0.0_dp)
1249 CALL set_diag(a, (1.0_dp, 0.0_dp))
1250
1251 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 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 c(1) = a(2)*b(3) - a(3)*b(2)
1267 c(2) = a(3)*b(1) - a(1)*b(3)
1268 c(3) = a(1)*b(2) - a(2)*b(1)
1269
1270 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 ELEMENTAL FUNCTION gcd(a, b)
1280 INTEGER, INTENT(IN) :: a, b
1281 INTEGER :: gcd
1282
1283 INTEGER :: aa, ab, l, rem, s
1284
1285 aa = abs(a)
1286 ab = abs(b)
1287 IF (aa < ab) THEN
1288 s = aa
1289 l = ab
1290 ELSE
1291 s = ab
1292 l = aa
1293 END IF
1294 IF (s .NE. 0) THEN
1295 DO
1296 rem = mod(l, s)
1297 IF (rem == 0) EXIT
1298 l = s
1299 s = rem
1300 END DO
1301 gcd = s
1302 ELSE
1303 gcd = l
1304 END IF
1305 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 ELEMENTAL FUNCTION lcm(a, b)
1315 INTEGER, INTENT(IN) :: a, b
1316 INTEGER :: lcm
1317
1318 INTEGER :: tmp
1319
1320 tmp = gcd(a, b)
1321 IF (tmp == 0) THEN
1322 lcm = 0
1323 ELSE
1324 ! could still overflow if the true lcm is larger than maxint
1325 lcm = abs((a/tmp)*b)
1326 END IF
1327 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 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 IF (x <= 0._dp) THEN
1347 cpabort("Invalid argument")
1348 END IF
1349
1350 IF (x < fpmin) THEN
1351 ei = log(x) + euler
1352 ELSE IF (x <= -log(eps)) THEN
1353 sum1 = 0._dp
1354 fact = 1._dp
1355 DO k = 1, maxit
1356 fact = fact*x/real(k, dp)
1357 term = fact/real(k, dp)
1358 sum1 = sum1 + term
1359 IF (term < eps*sum1) EXIT
1360 END DO
1361 ei = sum1 + log(x) + euler
1362 ELSE
1363 sum1 = 0._dp
1364 term = 1._dp
1365 DO k = 1, maxit
1366 prev = term
1367 term = term*real(k, dp)/x
1368 IF (term < eps) EXIT
1369 IF (term < prev) THEN
1370 sum1 = sum1 + term
1371 ELSE
1372 sum1 = sum1 - prev
1373 EXIT
1374 END IF
1375 END DO
1376 ei = exp(x)*(1._dp + sum1)/x
1377 END IF
1378
1379 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 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 nm1 = n - 1
1405
1406 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 cpabort("Invalid argument")
1408 ELSE IF (n .EQ. 0) THEN !Special case.
1409 expint = exp(-x)/x
1410 ELSE IF (x .EQ. 0.0_dp) THEN !Another special case.
1411 expint = 1.0_dp/nm1
1412 ELSE IF (x .GT. 1.0_dp) THEN !Lentz's algorithm (5.2).
1413 b = x + n
1414 c = 1.0_dp/fpmin
1415 d = 1.0_dp/b
1416 h = d
1417 DO i = 1, maxit
1418 a = -i*(nm1 + i)
1419 b = b + 2.0_dp
1420 d = 1.0_dp/(a*d + b)
1421 c = b + a/c
1422 del = c*d
1423 h = h*del
1424 IF (abs(del - 1.0_dp) .LT. eps) THEN
1425 expint = h*exp(-x)
1426 RETURN
1427 END IF
1428 END DO
1429 cpabort("continued fraction failed in expint")
1430 ELSE !Evaluate series.
1431 IF (nm1 .NE. 0) THEN !Set first term.
1432 expint = 1.0_dp/nm1
1433 ELSE
1434 expint = -log(x) - euler
1435 END IF
1436 fact = 1.0_dp
1437 DO i = 1, maxit
1438 fact = -fact*x/i
1439 IF (i .NE. nm1) THEN
1440 del = -fact/(i - nm1)
1441 ELSE
1442 psi = -euler !Compute I(n).
1443 DO ii = 1, nm1
1444 psi = psi + 1.0_dp/ii
1445 END DO
1446 del = fact*(-log(x) + psi)
1447 END IF
1448 expint = expint + del
1449 IF (abs(del) .LT. abs(expint)*eps) RETURN
1450 END DO
1451 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 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 n = SIZE(d(:))
1475
1476 ! Diagonalize matrix a
1477 CALL diag(n, a, d, v)
1478
1479 ! Sort eigenvalues and eigenvector in ascending order
1480 CALL eigsrt(n, d, v)
1481
1482 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 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 REAL(kind=dp), DIMENSION(n) :: b, z
1508
1509 CALL timeset(routinen, handle)
1510
1511 a_max = 0.0_dp
1512 DO ip = 1, n - 1
1513 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1514 b(ip) = a(ip, ip) ! get_diag(a)
1515 END DO
1516 b(n) = a(n, n)
1517
1518 CALL unit_matrix(v)
1519
1520 ! Go for 50 iterations
1521 DO i = 1, 50
1522 d = b
1523 d_min = max(d_eps, minval(abs(b)))
1524 IF (a_max < a_eps*d_min) THEN
1525 CALL timestop(handle)
1526 RETURN
1527 END IF
1528 tresh = merge(a_max, 0.0_dp, (i < 4))
1529 z = 0.0_dp
1530 DO ip = 1, n - 1
1531 DO iq = ip + 1, n
1532 dip = d(ip)
1533 diq = d(iq)
1534 apq = a(ip, iq)
1535 g = 100.0_dp*abs(apq)
1536 IF (tresh < abs(apq)) THEN
1537 h = diq - dip
1538 IF ((abs(h) + g) .NE. abs(h)) THEN
1539 theta = 0.5_dp*h/apq
1540 t = 1.0_dp/(abs(theta) + sqrt(1.0_dp + theta**2))
1541 IF (theta < 0.0_dp) t = -t
1542 ELSE
1543 t = apq/h
1544 END IF
1545 c = 1.0_dp/sqrt(1.0_dp + t**2)
1546 s = t*c
1547 tau = s/(1.0_dp + c)
1548 h = t*apq
1549 z(ip) = z(ip) - h
1550 z(iq) = z(iq) + h
1551 d(ip) = dip - h
1552 d(iq) = diq + h
1553 a(ip, iq) = 0.0_dp
1554 CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1555 CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1556 CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1557 CALL jrotate(v(:, ip), v(:, iq), s, tau)
1558 ELSE IF ((4 < i) .AND. &
1559 ((abs(dip) + g) == abs(dip)) .AND. &
1560 ((abs(diq) + g) == abs(diq))) THEN
1561 a(ip, iq) = 0.0_dp
1562 END IF
1563 END DO
1564 END DO
1565 b = b + z
1566 a_max = 0.0_dp
1567 DO ip = 1, n - 1
1568 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1569 END DO
1570 END DO
1571 WRITE (*, '(/,T2,A,/)') 'Too many iterations in jacobi'
1572
1573 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 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 u = 1.0_dp - ss*tt
1594 v = ss/u
1595
1596 a = a*u - b*ss
1597 b = b*(u + ss*v) + a*v
1598
1599 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 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 DO i = 1, n - 1
1619 j = sum(minloc(d(i:n))) + i - 1
1620 IF (j /= i) THEN
1621 CALL swap(d(i), d(j))
1622 CALL swap(v(:, i), v(:, j))
1623 END IF
1624 END DO
1625
1626 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 ELEMENTAL SUBROUTINE swap_scalar(a, b)
1637 REAL(kind=dp), INTENT(INOUT) :: a, b
1638
1639 REAL(kind=dp) :: c
1640
1641 c = a
1642 a = b
1643 b = c
1644
1645 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 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 n = SIZE(a)
1662
1663 IF (n /= SIZE(b)) THEN
1664 cpabort("Check the array bounds of the parameters")
1665 END IF
1666
1667 DO i = 1, n
1668 c = a(i)
1669 a(i) = b(i)
1670 b(i) = c
1671 END DO
1672
1673 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 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 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 r0 = sqrt(-log(eps*omg*10**2))/omg
1702 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1703
1704 DO iter = 1, itermax
1705 delta_r = f0/fprime0
1706 r0 = r0 - delta_r
1707 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1708 IF (abs(delta_r) .LT. abstol .OR. abs(f0) .LT. soltol) EXIT
1709 END DO
1710 cpassert(iter <= itermax)
1711 r_cutoff = r0
1712
1713 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 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 qr = omega*r
1730 fn = erfc(qr) - r*eps
1731 df = -2.0_dp*oorootpi*omega*exp(-qr**2) - eps
1732 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 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 COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1750 INTEGER :: handle, info, liwork, lrwork, lwork, n
1751 INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1752 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: rwork
1753
1754 CALL timeset(routinen, handle)
1755
1756 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 n = SIZE(matrix, 1)
1760 ALLOCATE (iwork(1), rwork(1), work(1))
1761
1762 ! work space query
1763 lwork = -1
1764 lrwork = -1
1765 liwork = -1
1766
1767 CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1768
1769 lwork = ceiling(real(work(1), kind=dp))
1770 lrwork = ceiling(rwork(1))
1771 liwork = iwork(1)
1772
1773 DEALLOCATE (iwork, rwork, work)
1774 ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
1775 eigenvectors(:, :) = matrix(:, :)
1776
1777 ! final diagonalization
1778 CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1779
1780 DEALLOCATE (iwork, rwork, work)
1781
1782 IF (info /= 0) cpabort("Diagonalisation of a complex matrix failed")
1783
1784 CALL timestop(handle)
1785
1786 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 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 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 n = SIZE(matrix, 1)
1807 ALLOCATE (matrix_c(n, n), eigenvalues(n))
1808
1809 matrix_c(:, :) = cmplx(0.0_dp, -matrix, kind=dp)
1810 CALL diag_complex(matrix_c, evecs, eigenvalues)
1811 evals = cmplx(0.0_dp, eigenvalues, kind=dp)
1812
1813 DEALLOCATE (matrix_c, eigenvalues)
1814 END SUBROUTINE diag_antisym
1815
1816END 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:206
mathlib::det_3x3
Definition mathlib.F:66
mathlib::invert_matrix
Definition mathlib.F:70
mathlib::set_diag
Definition mathlib.F:74
mathlib::unit_matrix
Definition mathlib.F:82
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::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
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:938
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:1468
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:206
mathlib::lcm
elemental integer function, public lcm(a, b)
computes the least common multiplier of two numbers
Definition mathlib.F:1315
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:1019
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:230
mathlib::invmat_symm
subroutine, public invmat_symm(a, potrf, uplo)
returns inverse of real symmetric, positive definite matrix
Definition mathlib.F:580
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:493
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:1086
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:176
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:1743
mathlib::diag_antisym
subroutine, public diag_antisym(matrix, evecs, evals)
Helper routine for diagonalizing anti symmetric matrices.
Definition mathlib.F:1795
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:286
mathlib::gcd
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1280
mathlib::symmetrize_matrix
subroutine, public symmetrize_matrix(a, option)
Symmetrize the matrix a.
Definition mathlib.F:1197
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:516
mathlib::invmat
subroutine, public invmat(a, info)
returns inverse of matrix using the lapack routines DGETRF and DGETRI
Definition mathlib.F:543
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:373
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:1496
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:1124
mathlib::erfc_cutoff
subroutine, public erfc_cutoff(eps, omg, r_cutoff)
compute a truncation radius for the shortrange operator
Definition mathlib.F:1686
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:151
mathlib::expint
elemental impure real(dp) function, public expint(n, x)
computes the exponential integral En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,...
Definition mathlib.F:1393
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:468
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:100
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:1263
mathlib::multinomial
pure real(kind=dp) function, public multinomial(n, k)
Calculates the multinomial coefficients.
Definition mathlib.F:255