(git:374b731)
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-2024 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 complex_diag
45
46 ! Public functions
47
48 PUBLIC :: angle, &
49 binomial, &
50 binomial_gen, &
51 multinomial, &
52 det_3x3, &
53 dihedral_angle, &
54 gcd, &
55 inv_3x3, &
56 lcm, &
57 vector_product, &
58 pswitch, &
59 rotate_vector, &
60 reflect_vector, &
61 expint, abnormal_value, &
62 get_diag, &
63 set_diag
64
65 INTERFACE det_3x3
66 MODULE PROCEDURE det_3x3_1, det_3x3_2
67 END INTERFACE
68
69 INTERFACE invert_matrix
70 MODULE PROCEDURE invert_matrix_d, invert_matrix_z
71 END INTERFACE
72
73 INTERFACE set_diag
74 MODULE PROCEDURE set_diag_scalar_d, set_diag_scalar_z
75 END INTERFACE
76
77 INTERFACE swap
78 MODULE PROCEDURE swap_scalar, swap_vector
79 END INTERFACE
80
81 INTERFACE unit_matrix
82 MODULE PROCEDURE unit_matrix_d, unit_matrix_z
83 END INTERFACE
84
85CONTAINS
86
87! **************************************************************************************************
88!> \brief Polynomial (5th degree) switching function
89!> f(a) = 1 .... f(b) = 0 with f'(a) = f"(a) = f'(b) = f"(b) = 0
90!> \param x ...
91!> \param a ...
92!> \param b ...
93!> \param order ...
94!> \return =0 : f(x)
95!> \return =1 : f'(x)
96!> \return =2 : f"(x)
97! **************************************************************************************************
98 FUNCTION pswitch(x, a, b, order) RESULT(fx)
99 REAL(kind=dp) :: x, a, b
100 INTEGER :: order
101 REAL(kind=dp) :: fx
102
103 REAL(kind=dp) :: u, u2, u3
104
105 cpassert(b > a)
106 IF (x < a .OR. x > b) THEN
107 ! outside switching intervall
108 IF (order > 0) THEN
109 ! derivatives are 0
110 fx = 0.0_dp
111 ELSE
112 IF (x < a) THEN
113 ! x < a => f(x) = 1
114 fx = 1.0_dp
115 ELSE
116 ! x > b => f(x) = 0
117 fx = 0.0_dp
118 END IF
119 END IF
120 ELSE
121 ! renormalized coordinate
122 u = (x - a)/(b - a)
123 SELECT CASE (order)
124 CASE (0)
125 u2 = u*u
126 u3 = u2*u
127 fx = 1._dp - 10._dp*u3 + 15._dp*u2*u2 - 6._dp*u2*u3
128 CASE (1)
129 u2 = u*u
130 fx = -30._dp*u2 + 60._dp*u*u2 - 30._dp*u2*u2
131 fx = fx/(b - a)
132 CASE (2)
133 u2 = u*u
134 fx = -60._dp*u + 180._dp*u2 - 120._dp*u*u2
135 fx = fx/(b - a)**2
136 CASE DEFAULT
137 cpabort('order not defined')
138 END SELECT
139 END IF
140
141 END FUNCTION pswitch
142
143! **************************************************************************************************
144!> \brief determines if a value is not normal (e.g. for Inf and Nan)
145!> based on IO to work also under optimization.
146!> \param a input value
147!> \return TRUE for NaN and Inf
148! **************************************************************************************************
149 LOGICAL FUNCTION abnormal_value(a)
150 REAL(kind=dp) :: a
151
152 CHARACTER(LEN=32) :: buffer
153
154 abnormal_value = .false.
155 ! the function should work when compiled with -ffast-math and similar
156 ! unfortunately, that option asserts that all numbers are normals,
157 ! which the compiler uses to optimize the function to .FALSE. if based on the IEEE module
158 ! therefore, pass this to the Fortran runtime/printf, if things are NaN or Inf, error out.
159 WRITE (buffer, *) a
160 IF (index(buffer, "N") .NE. 0 .OR. index(buffer, "n") .NE. 0) abnormal_value = .true.
161
162 END FUNCTION abnormal_value
163
164! **************************************************************************************************
165!> \brief Calculation of the angle between the vectors a and b.
166!> The angle is returned in radians.
167!> \param a ...
168!> \param b ...
169!> \return ...
170!> \date 14.10.1998
171!> \author MK
172!> \version 1.0
173! **************************************************************************************************
174 PURE FUNCTION angle(a, b) RESULT(angle_ab)
175 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: a, b
176 REAL(kind=dp) :: angle_ab
177
178 REAL(kind=dp) :: length_of_a, length_of_b
179 REAL(kind=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
180
181 length_of_a = sqrt(dot_product(a, a))
182 length_of_b = sqrt(dot_product(b, b))
183
184 IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
185 a_norm(:) = a(:)/length_of_a
186 b_norm(:) = b(:)/length_of_b
187 angle_ab = acos(min(max(dot_product(a_norm, b_norm), -1.0_dp), 1.0_dp))
188 ELSE
189 angle_ab = 0.0_dp
190 END IF
191
192 END FUNCTION angle
193
194! **************************************************************************************************
195!> \brief The binomial coefficient n over k for 0 <= k <= n is calculated,
196!> otherwise zero is returned.
197!> \param n ...
198!> \param k ...
199!> \return ...
200!> \date 08.03.1999
201!> \author MK
202!> \version 1.0
203! **************************************************************************************************
204 ELEMENTAL FUNCTION binomial(n, k) RESULT(n_over_k)
205 INTEGER, INTENT(IN) :: n, k
206 REAL(kind=dp) :: n_over_k
207
208 IF ((k >= 0) .AND. (k <= n)) THEN
209 n_over_k = fac(n)/(fac(n - k)*fac(k))
210 ELSE
211 n_over_k = 0.0_dp
212 END IF
213
214 END FUNCTION binomial
215
216! **************************************************************************************************
217!> \brief The generalized binomial coefficient z over k for 0 <= k <= n is calculated.
218!> (z) z*(z-1)*...*(z-k+2)*(z-k+1)
219!> ( ) = ---------------------------
220!> (k) k!
221!> \param z ...
222!> \param k ...
223!> \return ...
224!> \date 11.11.2019
225!> \author FS
226!> \version 1.0
227! **************************************************************************************************
228 ELEMENTAL FUNCTION binomial_gen(z, k) RESULT(z_over_k)
229 REAL(kind=dp), INTENT(IN) :: z
230 INTEGER, INTENT(IN) :: k
231 REAL(kind=dp) :: z_over_k
232
233 INTEGER :: i
234
235 IF (k >= 0) THEN
236 z_over_k = 1.0_dp
237 DO i = 1, k
238 z_over_k = z_over_k*(z - i + 1)/real(i, dp)
239 END DO
240 ELSE
241 z_over_k = 0.0_dp
242 END IF
243
244 END FUNCTION binomial_gen
245
246! **************************************************************************************************
247!> \brief Calculates the multinomial coefficients
248!> \param n ...
249!> \param k ...
250!> \return ...
251!> \author Ole Schuett
252! **************************************************************************************************
253 PURE FUNCTION multinomial(n, k) RESULT(res)
254 INTEGER, INTENT(IN) :: n
255 INTEGER, DIMENSION(:), INTENT(IN) :: k
256 REAL(kind=dp) :: res
257
258 INTEGER :: i
259 REAL(kind=dp) :: denom
260
261 IF (all(k >= 0) .AND. sum(k) == n) THEN
262 denom = 1.0_dp
263 DO i = 1, SIZE(k)
264 denom = denom*fac(k(i))
265 END DO
266 res = fac(n)/denom
267 ELSE
268 res = 0.0_dp
269 END IF
270
271 END FUNCTION multinomial
272
273! **************************************************************************************************
274!> \brief The rotation matrix rotmat which rotates a vector about a
275!> rotation axis defined by the vector a is build up.
276!> The rotation angle is phi (radians).
277!> \param phi ...
278!> \param a ...
279!> \param rotmat ...
280!> \date 16.10.1998
281!> \author MK
282!> \version 1.0
283! **************************************************************************************************
284 PURE SUBROUTINE build_rotmat(phi, a, rotmat)
285 REAL(kind=dp), INTENT(IN) :: phi
286 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
287 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: rotmat
288
289 REAL(kind=dp) :: cosp, cost, length_of_a, sinp
290 REAL(kind=dp), DIMENSION(3) :: d
291
292 length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
293 ! Check the length of the vector a
294 IF (length_of_a > eps_geo) THEN
295
296 d(:) = a(:)/length_of_a
297
298 cosp = cos(phi)
299 sinp = sin(phi)
300 cost = 1.0_dp - cosp
301
302 rotmat(1, 1) = d(1)*d(1)*cost + cosp
303 rotmat(1, 2) = d(1)*d(2)*cost - d(3)*sinp
304 rotmat(1, 3) = d(1)*d(3)*cost + d(2)*sinp
305 rotmat(2, 1) = d(2)*d(1)*cost + d(3)*sinp
306 rotmat(2, 2) = d(2)*d(2)*cost + cosp
307 rotmat(2, 3) = d(2)*d(3)*cost - d(1)*sinp
308 rotmat(3, 1) = d(3)*d(1)*cost - d(2)*sinp
309 rotmat(3, 2) = d(3)*d(2)*cost + d(1)*sinp
310 rotmat(3, 3) = d(3)*d(3)*cost + cosp
311 ELSE
312 CALL unit_matrix(rotmat)
313 END IF
314
315 END SUBROUTINE build_rotmat
316
317! **************************************************************************************************
318!> \brief Returns the determinante of the 3x3 matrix a.
319!> \param a ...
320!> \return ...
321!> \date 13.03.2004
322!> \author MK
323!> \version 1.0
324! **************************************************************************************************
325 PURE FUNCTION det_3x3_1(a) RESULT(det_a)
326 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: a
327 REAL(kind=dp) :: det_a
328
329 det_a = a(1, 1)*(a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2)) + &
330 a(1, 2)*(a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3)) + &
331 a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
332
333 END FUNCTION det_3x3_1
334
335! **************************************************************************************************
336!> \brief Returns the determinante of the 3x3 matrix a given by its columns.
337!> \param a1 ...
338!> \param a2 ...
339!> \param a3 ...
340!> \return ...
341!> \date 13.03.2004
342!> \author MK
343!> \version 1.0
344! **************************************************************************************************
345 PURE FUNCTION det_3x3_2(a1, a2, a3) RESULT(det_a)
346 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a1, a2, a3
347 REAL(kind=dp) :: det_a
348
349 det_a = a1(1)*(a2(2)*a3(3) - a3(2)*a2(3)) + &
350 a2(1)*(a3(2)*a1(3) - a1(2)*a3(3)) + &
351 a3(1)*(a1(2)*a2(3) - a2(2)*a1(3))
352
353 END FUNCTION det_3x3_2
354
355! **************************************************************************************************
356!> \brief Diagonalize the symmetric n by n matrix a using the LAPACK
357!> library. Only the upper triangle of matrix a is used.
358!> Externals (LAPACK 3.0)
359!> \param a ...
360!> \param eigval ...
361!> \param dac ...
362!> \date 29.03.1999
363!> \par Variables
364!> - a : Symmetric matrix to be diagonalized (input; upper triangle) ->
365!> - eigenvectors of the matrix a (output).
366!> - dac : If true, then the divide-and-conquer algorithm is applied.
367!> - eigval : Eigenvalues of the matrix a (output).
368!> \author MK
369!> \version 1.0
370! **************************************************************************************************
371 SUBROUTINE diamat_all(a, eigval, dac)
372 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
373 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigval
374 LOGICAL, INTENT(IN), OPTIONAL :: dac
375
376 CHARACTER(len=*), PARAMETER :: routinen = 'diamat_all'
377
378 INTEGER :: handle, info, liwork, lwork, n, nb
379 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
380 INTEGER, EXTERNAL :: ilaenv
381 LOGICAL :: divide_and_conquer
382 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: work
383
384 EXTERNAL dsyev, dsyevd
385
386 CALL timeset(routinen, handle)
387
388 ! Get the size of the matrix a
389 n = SIZE(a, 1)
390
391 ! Check the size of matrix a
392 IF (SIZE(a, 2) /= n) THEN
393 cpabort("Check the size of matrix a (parameter #1)")
394 END IF
395
396 ! Check the size of vector eigval
397 IF (SIZE(eigval) /= n) THEN
398 cpabort("The dimension of vector eigval is too small")
399 END IF
400
401 ! Check, if the divide-and-conquer algorithm is requested
402
403 IF (PRESENT(dac)) THEN
404 divide_and_conquer = dac
405 ELSE
406 divide_and_conquer = .false.
407 END IF
408
409 ! Get the optimal work storage size
410
411 IF (divide_and_conquer) THEN
412 lwork = 2*n**2 + 6*n + 1
413 liwork = 5*n + 3
414 ELSE
415 nb = ilaenv(1, "DSYTRD", "U", n, -1, -1, -1)
416 lwork = (nb + 2)*n
417 END IF
418
419 ! Allocate work storage
420
421 ALLOCATE (work(lwork))
422 IF (divide_and_conquer) THEN
423 ALLOCATE (iwork(liwork))
424 END IF
425
426 ! Diagonalize the matrix a
427
428 info = 0
429 IF (divide_and_conquer) THEN
430 CALL dsyevd("V", "U", n, a, n, eigval, work, lwork, iwork, liwork, info)
431 ELSE
432 CALL dsyev("V", "U", n, a, n, eigval, work, lwork, info)
433 END IF
434
435 IF (info /= 0) THEN
436 IF (divide_and_conquer) THEN
437 cpabort("The matrix diagonalization with dsyevd failed")
438 ELSE
439 cpabort("The matrix diagonalization with dsyev failed")
440 END IF
441 END IF
442
443 ! Release work storage
444
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 INTEGER :: lwork, n
547 INTEGER, ALLOCATABLE :: ipiv(:)
548 REAL(kind=dp), ALLOCATABLE :: work(:)
549
550 n = SIZE(a, 1)
551 lwork = 20*n
552 ALLOCATE (ipiv(n))
553 ALLOCATE (work(lwork))
554 ipiv = 0
555 work = 0._dp
556 info = 0
557 CALL dgetrf(n, n, a, n, ipiv, info)
558 IF (info == 0) THEN
559 CALL dgetri(n, a, n, ipiv, work, lwork, info)
560 END IF
561 DEALLOCATE (ipiv, work)
562 END SUBROUTINE invmat
563
564! **************************************************************************************************
565!> \brief returns inverse of real symmetric, positive definite matrix
566!> \param a matrix
567!> \param cholesky_triangle if cholesky decomposition of a was already done
568!> using dpotrf, indicating if the upper or lower triangle of a is
569!> stored. If not given, cholesky decomposition of a will be done
570!> before inversion.
571!> \author Dorothea Golze [02.2015]
572! **************************************************************************************************
573 SUBROUTINE invmat_symm(a, cholesky_triangle)
574 REAL(kind=dp), INTENT(INOUT) :: a(:, :)
575 CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: cholesky_triangle
576
577 CHARACTER(LEN=*), PARAMETER :: routinen = 'invmat_symm'
578
579 CHARACTER(LEN=1) :: my_triangle
580 INTEGER :: handle, i, info, n
581
582 CALL timeset(routinen, handle)
583
584 n = SIZE(a, 1)
585 info = 0
586
587 IF (PRESENT(cholesky_triangle)) THEN
588 my_triangle = cholesky_triangle
589 ELSE
590 my_triangle = "U"
591 END IF
592
593 ! do cholesky decomposition
594 IF (.NOT. PRESENT(cholesky_triangle)) THEN
595 CALL dpotrf(my_triangle, n, a, n, info)
596 IF (info /= 0) THEN
597 cpabort("DPOTRF failed")
598 END IF
599 END IF
600
601 ! do inversion using the cholesky decomposition
602 CALL dpotri(my_triangle, n, a, n, info)
603 IF (info /= 0) THEN
604 cpabort("Matrix inversion failed")
605 END IF
606
607 IF (my_triangle == "U") THEN
608 DO i = 1, n - 1
609 a(i + 1:n, i) = a(i, i + 1:n)
610 END DO
611 ELSE
612 DO i = 1, n - 1
613 a(i, i + 1:n) = a(i + 1:n, i)
614 END DO
615 END IF
616
617 CALL timestop(handle)
618
619 END SUBROUTINE invmat_symm
620
621! **************************************************************************************************
622!> \brief Compute the inverse of the n by n real matrix a using the LAPACK
623!> library
624!> \param a ...
625!> \param a_inverse ...
626!> \param eval_error ...
627!> \param option ...
628!> \param improve ...
629!> \date 23.03.1999
630!> \par Variables
631!> - a : Real matrix to be inverted (input).
632!> - a_inverse: Inverse of the matrix a (output).
633!> - a_lu : LU factorization of matrix a.
634!> - a_norm : Norm of matrix a.
635!> - error : Estimated error of the inversion.
636!> - r_cond : Reciprocal condition number of the matrix a.
637!> - trans : "N" => invert a
638!> - "T" => invert transpose(a)
639!> \author MK
640!> \version 1.0
641!> \note NB add improve argument, used to disable call to dgerfs
642! **************************************************************************************************
643 SUBROUTINE invert_matrix_d(a, a_inverse, eval_error, option, improve)
644 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
645 REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: a_inverse
646 REAL(KIND=dp), INTENT(OUT) :: eval_error
647 CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: option
648 LOGICAL, INTENT(IN), OPTIONAL :: improve
649
650 CHARACTER(LEN=1) :: norm, trans
651 CHARACTER(LEN=default_string_length) :: message
652 INTEGER :: info, iter, n
653 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv, iwork
654 LOGICAL :: do_improve
655 REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
656 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, work
657 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
658 REAL(KIND=dp), EXTERNAL :: dlange
659
660 EXTERNAL dgecon, dgerfs, dgetrf, dgetrs
661
662 ! Check for optional parameter
663 IF (PRESENT(option)) THEN
664 trans = option
665 ELSE
666 trans = "N"
667 END IF
668
669 IF (PRESENT(improve)) THEN
670 do_improve = improve
671 ELSE
672 do_improve = .true.
673 END IF
674
675 ! Get the dimension of matrix a
676 n = SIZE(a, 1)
677
678 ! Check array dimensions
679 IF (n == 0) THEN
680 cpabort("Matrix to be inverted of zero size")
681 END IF
682
683 IF (n /= SIZE(a, 2)) THEN
684 cpabort("Check the array bounds of parameter #1")
685 END IF
686
687 IF ((n /= SIZE(a_inverse, 1)) .OR. &
688 (n /= SIZE(a_inverse, 2))) THEN
689 cpabort("Check the array bounds of parameter #2")
690 END IF
691
692 ! Allocate work storage
693 ALLOCATE (a_lu(n, n))
694
695 ALLOCATE (b(n, n))
696
697 ALLOCATE (berr(n))
698
699 ALLOCATE (ferr(n))
700
701 ALLOCATE (ipiv(n))
702
703 ALLOCATE (iwork(n))
704
705 ALLOCATE (work(4*n))
706
707 a_lu(1:n, 1:n) = a(1:n, 1:n)
708
709 ! Compute the LU factorization of the matrix a
710 CALL dgetrf(n, n, a_lu, n, ipiv, info)
711
712 IF (info /= 0) THEN
713 cpabort("The LU factorization in dgetrf failed")
714 END IF
715
716 ! Compute the norm of the matrix a
717
718 IF (trans == "N") THEN
719 norm = '1'
720 ELSE
721 norm = 'I'
722 END IF
723
724 a_norm = dlange(norm, n, n, a, n, work)
725
726 ! Compute the reciprocal of the condition number of a
727
728 CALL dgecon(norm, n, a_lu, n, a_norm, r_cond, work, iwork, info)
729
730 IF (info /= 0) THEN
731 cpabort("The computation of the condition number in dgecon failed")
732 END IF
733
734 IF (r_cond < epsilon(0.0_dp)) THEN
735 WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
736 CALL cp_abort(__location__, &
737 "Bad condition number "//trim(message)//" (smaller than the machine "// &
738 "working precision)")
739 END IF
740
741 ! Solve a system of linear equations using the LU factorization computed by dgetrf
742
743 CALL unit_matrix(a_inverse)
744
745 CALL dgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
746
747 IF (info /= 0) THEN
748 cpabort("Solving the system of linear equations in dgetrs failed")
749 END IF
750
751 ! Improve the computed solution iteratively
752 CALL unit_matrix(b) ! Initialize right-hand sides
753
754 eval_error = 0.0_dp
755
756 IF (do_improve) THEN
757 DO iter = 1, 10
758
759 CALL dgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
760 work, iwork, info)
761
762 IF (info /= 0) THEN
763 cpabort("Improving the computed solution in dgerfs failed")
764 END IF
765
766 old_eval_error = eval_error
767 eval_error = maxval(ferr)
768
769 IF (abs(eval_error - old_eval_error) <= epsilon(1.0_dp)) EXIT
770
771 END DO
772 END IF
773
774 ! Release work storage
775 DEALLOCATE (work)
776 DEALLOCATE (iwork)
777 DEALLOCATE (ipiv)
778 DEALLOCATE (ferr)
779 DEALLOCATE (berr)
780 DEALLOCATE (b)
781 DEALLOCATE (a_lu)
782
783 END SUBROUTINE invert_matrix_d
784
785! **************************************************************************************************
786!> \brief Compute the inverse of the n by n complex matrix a using the LAPACK
787!> library
788!> \param a ...
789!> \param a_inverse ...
790!> \param eval_error ...
791!> \param option ...
792!> \date 08.06.2009
793!> \par Variables
794!> - a : Complex matrix to be inverted (input).
795!> - a_inverse: Inverse of the matrix a (output).
796!> - a_lu : LU factorization of matrix a.
797!> - a_norm : Norm of matrix a.
798!> - error : Estimated error of the inversion.
799!> - r_cond : Reciprocal condition number of the matrix a.
800!> - trans : "N" => invert a
801!> - "T" => invert transpose(a)
802!> \author MK
803!> \version 1.0
804! **************************************************************************************************
805 SUBROUTINE invert_matrix_z(a, a_inverse, eval_error, option)
806 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: a
807 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: a_inverse
808 REAL(KIND=dp), INTENT(OUT) :: eval_error
809 CHARACTER(LEN=1), INTENT(IN), OPTIONAL :: option
810
811 CHARACTER(LEN=1) :: norm, trans
812 CHARACTER(LEN=default_string_length) :: message
813 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
814 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
815 INTEGER :: info, iter, n
816 INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
817 REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
818 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, rwork
819 REAL(KIND=dp), EXTERNAL :: zlange
820
821 EXTERNAL zgecon, zgerfs, zgetrf, zgetrs
822
823 ! Check for optional parameter
824 IF (PRESENT(option)) THEN
825 trans = option
826 ELSE
827 trans = "N"
828 END IF
829
830 ! Get the dimension of matrix a
831 n = SIZE(a, 1)
832
833 ! Check array dimensions
834 IF (n == 0) THEN
835 cpabort("Matrix to be inverted of zero size")
836 END IF
837
838 IF (n /= SIZE(a, 2)) THEN
839 cpabort("Check the array bounds of parameter #1")
840 END IF
841
842 IF ((n /= SIZE(a_inverse, 1)) .OR. &
843 (n /= SIZE(a_inverse, 2))) THEN
844 cpabort("Check the array bounds of parameter #2")
845 END IF
846
847 ! Allocate work storage
848 ALLOCATE (a_lu(n, n))
849
850 ALLOCATE (b(n, n))
851
852 ALLOCATE (berr(n))
853
854 ALLOCATE (ferr(n))
855
856 ALLOCATE (ipiv(n))
857
858 ALLOCATE (rwork(2*n))
859
860 ALLOCATE (work(2*n))
861
862 a_lu(1:n, 1:n) = a(1:n, 1:n)
863
864 ! Compute the LU factorization of the matrix a
865 CALL zgetrf(n, n, a_lu, n, ipiv, info)
866
867 IF (info /= 0) THEN
868 cpabort("The LU factorization in dgetrf failed")
869 END IF
870
871 ! Compute the norm of the matrix a
872
873 IF (trans == "N") THEN
874 norm = '1'
875 ELSE
876 norm = 'I'
877 END IF
878
879 a_norm = zlange(norm, n, n, a, n, work)
880
881 ! Compute the reciprocal of the condition number of a
882
883 CALL zgecon(norm, n, a_lu, n, a_norm, r_cond, work, rwork, info)
884
885 IF (info /= 0) THEN
886 cpabort("The computation of the condition number in dgecon failed")
887 END IF
888
889 IF (r_cond < epsilon(0.0_dp)) THEN
890 WRITE (message, "(A,ES10.3)") "R_COND =", r_cond
891 CALL cp_abort(__location__, &
892 "Bad condition number "//trim(message)//" (smaller than the machine "// &
893 "working precision)")
894 END IF
895
896 ! Solve a system of linear equations using the LU factorization computed by dgetrf
897
898 CALL unit_matrix(a_inverse)
899
900 CALL zgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
901
902 IF (info /= 0) THEN
903 cpabort("Solving the system of linear equations in dgetrs failed")
904 END IF
905
906 ! Improve the computed solution iteratively
907 CALL unit_matrix(b) ! Initialize right-hand sides
908
909 eval_error = 0.0_dp
910
911 DO iter = 1, 10
912
913 CALL zgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
914 work, rwork, info)
915
916 IF (info /= 0) THEN
917 cpabort("Improving the computed solution in dgerfs failed")
918 END IF
919
920 old_eval_error = eval_error
921 eval_error = maxval(ferr)
922
923 IF (abs(eval_error - old_eval_error) <= epsilon(1.0_dp)) EXIT
924
925 END DO
926
927 ! Release work storage
928 DEALLOCATE (work)
929 DEALLOCATE (rwork)
930 DEALLOCATE (ipiv)
931 DEALLOCATE (ferr)
932 DEALLOCATE (berr)
933 DEALLOCATE (b)
934 DEALLOCATE (a_lu)
935
936 END SUBROUTINE invert_matrix_z
937
938! **************************************************************************************************
939!> \brief returns the pseudoinverse of a real, square matrix using singular
940!> value decomposition
941!> \param a matrix a
942!> \param a_pinverse pseudoinverse of matrix a
943!> \param rskip parameter for setting small singular values to zero
944!> \param determinant determinant of matrix a (optional output)
945!> \param sval array holding singular values of matrix a (optional output)
946!> \author Dorothea Golze [02.2015]
947! **************************************************************************************************
948 SUBROUTINE get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
949 REAL(kind=dp), DIMENSION(:, :) :: a, a_pinverse
950 REAL(kind=dp), INTENT(IN) :: rskip
951 REAL(kind=dp), INTENT(OUT), OPTIONAL :: determinant
952 REAL(kind=dp), DIMENSION(:), INTENT(INOUT), &
953 OPTIONAL, POINTER :: sval
954
955 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_pseudo_inverse_svd'
956
957 INTEGER :: handle, i, info, lwork, n
958 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
959 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: sig, work
960 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: sig_plus, temp_mat, u, vt
961
962 CALL timeset(routinen, handle)
963
964 n = SIZE(a, 1)
965 ALLOCATE (u(n, n), vt(n, n), sig(n), sig_plus(n, n), iwork(8*n), work(1), temp_mat(n, n))
966 u(:, :) = 0.0_dp
967 vt(:, :) = 0.0_dp
968 sig(:) = 0.0_dp
969 sig_plus = 0.0_dp
970 work = 0.0_dp
971 iwork = 0
972 IF (PRESENT(determinant)) determinant = 1.0_dp
973
974 ! work size query
975 lwork = -1
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("ERROR in DGESDD: Could not retrieve work array sizes")
981 END IF
982 lwork = int(work(1))
983 DEALLOCATE (work)
984 ALLOCATE (work(lwork))
985
986 ! do SVD
987 CALL dgesdd('A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
988 lwork, iwork(1), info)
989
990 IF (info /= 0) THEN
991 cpabort("SVD failed")
992 END IF
993
994 IF (PRESENT(sval)) THEN
995 cpassert(.NOT. ASSOCIATED(sval))
996 ALLOCATE (sval(n))
997 sval(:) = sig
998 END IF
999
1000 ! set singular values that are too small to zero
1001 DO i = 1, n
1002 IF (sig(i) > rskip*maxval(sig)) THEN
1003 IF (PRESENT(determinant)) &
1004 determinant = determinant*sig(i)
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 .NE. 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 exponential integral
1394!> En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,1,..
1395!> Note: Ei(-x) = -E1(x)
1396!> \param n ...
1397!> \param x ...
1398!> \return ...
1399!> \par History
1400!> 05.2007 Created
1401!> \author Manuel Guidon (adapted from Numerical recipies)
1402! **************************************************************************************************
1403 ELEMENTAL IMPURE FUNCTION expint(n, x)
1404 INTEGER, INTENT(IN) :: n
1405 REAL(dp), INTENT(IN) :: x
1406 REAL(dp) :: expint
1407
1408 INTEGER, PARAMETER :: maxit = 100
1409 REAL(dp), PARAMETER :: eps = 6.e-14_dp, euler = 0.5772156649015328606065120_dp, &
1410 fpmin = tiny(0.0_dp)
1411
1412 INTEGER :: i, ii, nm1
1413 REAL(dp) :: a, b, c, d, del, fact, h, psi
1414
1415 nm1 = n - 1
1416
1417 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
1418 cpabort("Invalid argument")
1419 ELSE IF (n .EQ. 0) THEN !Special case.
1420 expint = exp(-x)/x
1421 ELSE IF (x .EQ. 0.0_dp) THEN !Another special case.
1422 expint = 1.0_dp/nm1
1423 ELSE IF (x .GT. 1.0_dp) THEN !Lentz's algorithm (5.2).
1424 b = x + n
1425 c = 1.0_dp/fpmin
1426 d = 1.0_dp/b
1427 h = d
1428 DO i = 1, maxit
1429 a = -i*(nm1 + i)
1430 b = b + 2.0_dp
1431 d = 1.0_dp/(a*d + b)
1432 c = b + a/c
1433 del = c*d
1434 h = h*del
1435 IF (abs(del - 1.0_dp) .LT. eps) THEN
1436 expint = h*exp(-x)
1437 RETURN
1438 END IF
1439 END DO
1440 cpabort("continued fraction failed in expint")
1441 ELSE !Evaluate series.
1442 IF (nm1 .NE. 0) THEN !Set first term.
1443 expint = 1.0_dp/nm1
1444 ELSE
1445 expint = -log(x) - euler
1446 END IF
1447 fact = 1.0_dp
1448 DO i = 1, maxit
1449 fact = -fact*x/i
1450 IF (i .NE. nm1) THEN
1451 del = -fact/(i - nm1)
1452 ELSE
1453 psi = -euler !Compute I(n).
1454 DO ii = 1, nm1
1455 psi = psi + 1.0_dp/ii
1456 END DO
1457 del = fact*(-log(x) + psi)
1458 END IF
1459 expint = expint + del
1460 IF (abs(del) .LT. abs(expint)*eps) RETURN
1461 END DO
1462 cpabort("series failed in expint")
1463 END IF
1464
1465 END FUNCTION expint
1466
1467! **************************************************************************************************
1468!> \brief Jacobi matrix diagonalization. The eigenvalues are returned in
1469!> vector d and the eigenvectors are returned in matrix v in ascending
1470!> order.
1471!>
1472!> \param a ...
1473!> \param d ...
1474!> \param v ...
1475!> \par History
1476!> - Creation (20.11.98, Matthias Krack)
1477! **************************************************************************************************
1478 SUBROUTINE jacobi(a, d, v)
1479 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1480 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: d
1481 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: v
1482
1483 INTEGER :: n
1484
1485 n = SIZE(d(:))
1486
1487 ! Diagonalize matrix a
1488 CALL diag(n, a, d, v)
1489
1490 ! Sort eigenvalues and eigenvector in ascending order
1491 CALL eigsrt(n, d, v)
1492
1493 END SUBROUTINE jacobi
1494
1495! **************************************************************************************************
1496!> \brief Diagonalize matrix a. The eigenvalues are returned in vector d
1497!> and the eigenvectors are returned in matrix v.
1498!>
1499!> \param n matrix/vector extent (problem size)
1500!> \param a matrix to be diagonalised
1501!> \param d vector of eigenvalues
1502!> \param v matrix of eigenvectors
1503!> \par History
1504!> - Creation (20.11.98, Matthias Krack)
1505! **************************************************************************************************
1506 SUBROUTINE diag(n, a, d, v)
1507 INTEGER, INTENT(IN) :: n
1508 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1509 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: d
1510 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: v
1511
1512 REAL(kind=dp), PARAMETER :: a_eps = 1.0e-10_dp, d_eps = 1.0e-3_dp
1513
1514 INTEGER :: i, ip, iq
1515 REAL(kind=dp) :: a_max, apq, c, d_min, dip, diq, g, h, s, &
1516 t, tau, theta, tresh
1517 REAL(kind=dp), DIMENSION(n) :: b, z
1518
1519 a_max = 0.0_dp
1520 DO ip = 1, n - 1
1521 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1522 b(ip) = a(ip, ip) ! get_diag(a)
1523 END DO
1524 b(n) = a(n, n)
1525
1526 CALL unit_matrix(v)
1527
1528 ! Go for 50 iterations
1529 DO i = 1, 50
1530 d = b
1531 d_min = max(d_eps, minval(abs(b)))
1532 IF (a_max < a_eps*d_min) RETURN
1533 tresh = merge(a_max, 0.0_dp, (i < 4))
1534 z = 0.0_dp
1535 DO ip = 1, n - 1
1536 DO iq = ip + 1, n
1537 dip = d(ip)
1538 diq = d(iq)
1539 apq = a(ip, iq)
1540 g = 100.0_dp*abs(apq)
1541 IF (tresh < abs(apq)) THEN
1542 h = diq - dip
1543 IF ((abs(h) + g) .NE. abs(h)) THEN
1544 theta = 0.5_dp*h/apq
1545 t = 1.0_dp/(abs(theta) + sqrt(1.0_dp + theta**2))
1546 IF (theta < 0.0_dp) t = -t
1547 ELSE
1548 t = apq/h
1549 END IF
1550 c = 1.0_dp/sqrt(1.0_dp + t**2)
1551 s = t*c
1552 tau = s/(1.0_dp + c)
1553 h = t*apq
1554 z(ip) = z(ip) - h
1555 z(iq) = z(iq) + h
1556 d(ip) = dip - h
1557 d(iq) = diq + h
1558 a(ip, iq) = 0.0_dp
1559 CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1560 CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1561 CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1562 CALL jrotate(v(:, ip), v(:, iq), s, tau)
1563 ELSE IF ((4 < i) .AND. &
1564 ((abs(dip) + g) == abs(dip)) .AND. &
1565 ((abs(diq) + g) == abs(diq))) THEN
1566 a(ip, iq) = 0.0_dp
1567 END IF
1568 END DO
1569 END DO
1570 b = b + z
1571 a_max = 0.0_dp
1572 DO ip = 1, n - 1
1573 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1574 END DO
1575 END DO
1576 WRITE (*, '(/,T2,A,/)') 'Too many iterations in jacobi'
1577
1578 END SUBROUTINE diag
1579
1580! **************************************************************************************************
1581!> \brief Perform a Jacobi rotation of the vectors a and b.
1582!>
1583!> \param a ...
1584!> \param b ...
1585!> \param ss ...
1586!> \param tt ...
1587!> \par History
1588!> - Creation (20.11.98, Matthias Krack)
1589! **************************************************************************************************
1590 PURE SUBROUTINE jrotate(a, b, ss, tt)
1591 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1592 REAL(kind=dp), INTENT(IN) :: ss, tt
1593
1594 REAL(kind=dp) :: u, v
1595
1596 u = 1.0_dp - ss*tt
1597 v = ss/u
1598
1599 a = a*u - b*ss
1600 b = b*(u + ss*v) + a*v
1601
1602 END SUBROUTINE jrotate
1603
1604! **************************************************************************************************
1605!> \brief Sort the values in vector d in ascending order and swap the
1606!> corresponding columns of matrix v.
1607!>
1608!> \param n ...
1609!> \param d ...
1610!> \param v ...
1611!> \par History
1612!> - Creation (20.11.98, Matthias Krack)
1613! **************************************************************************************************
1614 SUBROUTINE eigsrt(n, d, v)
1615 INTEGER, INTENT(IN) :: n
1616 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: d
1617 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: v
1618
1619 INTEGER :: i, j
1620
1621 DO i = 1, n - 1
1622 j = sum(minloc(d(i:n))) + i - 1
1623 IF (j /= i) THEN
1624 CALL swap(d(i), d(j))
1625 CALL swap(v(:, i), v(:, j))
1626 END IF
1627 END DO
1628
1629 END SUBROUTINE eigsrt
1630
1631! **************************************************************************
1632!> \brief Swap two scalars
1633!>
1634!> \param a ...
1635!> \param b ...
1636!> \par History
1637!> - Creation (20.11.98, Matthias Krack)
1638! **************************************************************************************************
1639 ELEMENTAL SUBROUTINE swap_scalar(a, b)
1640 REAL(kind=dp), INTENT(INOUT) :: a, b
1641
1642 REAL(kind=dp) :: c
1643
1644 c = a
1645 a = b
1646 b = c
1647
1648 END SUBROUTINE swap_scalar
1649
1650! **************************************************************************
1651!> \brief Swap two vectors
1652!>
1653!> \param a ...
1654!> \param b ...
1655!> \par History
1656!> - Creation (20.11.98, Matthias Krack)
1657! **************************************************************************************************
1658 SUBROUTINE swap_vector(a, b)
1659 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1660
1661 INTEGER :: i, n
1662 REAL(kind=dp) :: c
1663
1664 n = SIZE(a)
1665
1666 IF (n /= SIZE(b)) THEN
1667 cpabort("Check the array bounds of the parameters")
1668 END IF
1669
1670 DO i = 1, n
1671 c = a(i)
1672 a(i) = b(i)
1673 b(i) = c
1674 END DO
1675
1676 END SUBROUTINE swap_vector
1677
1678! **************************************************************************************************
1679!> \brief - compute a truncation radius for the shortrange operator
1680!> \param eps target accuracy!> \param omg screening parameter
1681!> \param omg ...
1682!> \param r_cutoff cutoff radius
1683!> \par History
1684!> 10.2012 created [Hossein Banihashemian]
1685!> 05.2019 moved here from hfx_types (A. Bussy)
1686!> \author Hossein Banihashemian
1687! **************************************************************************************************
1688 SUBROUTINE erfc_cutoff(eps, omg, r_cutoff)
1689 IMPLICIT NONE
1690
1691 REAL(dp), INTENT(in) :: eps, omg
1692 REAL(dp), INTENT(out) :: r_cutoff
1693
1694 CHARACTER(LEN=*), PARAMETER :: routinen = 'erfc_cutoff'
1695
1696 REAL(dp), PARAMETER :: abstol = 1e-10_dp, soltol = 1e-16_dp
1697 REAL(dp) :: r0, f0, fprime0, delta_r
1698 INTEGER :: iter, handle
1699 INTEGER, PARAMETER :: itermax = 1000
1700
1701 CALL timeset(routinen, handle)
1702
1703 ! initial guess assuming that we are in the asymptotic regime of the erf, and the solution is about 10.
1704 r0 = sqrt(-log(eps*omg*10**2))/omg
1705 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1706
1707 DO iter = 1, itermax
1708 delta_r = f0/fprime0
1709 r0 = r0 - delta_r
1710 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1711 IF (abs(delta_r) .LT. abstol .OR. abs(f0) .LT. soltol) EXIT
1712 END DO
1713 cpassert(iter <= itermax)
1714 r_cutoff = r0
1715
1716 CALL timestop(handle)
1717 CONTAINS
1718! **************************************************************************************************
1719!> \brief ...
1720!> \param r ...
1721!> \param eps ...
1722!> \param omega ...
1723!> \param fn ...
1724!> \param df ...
1725! **************************************************************************************************
1726 ELEMENTAL SUBROUTINE eval_transc_func(r, eps, omega, fn, df)
1727 REAL(dp), INTENT(in) :: r, eps, omega
1728 REAL(dp), INTENT(out) :: fn, df
1729
1730 REAL(dp) :: qr
1731
1732 qr = omega*r
1733 fn = erfc(qr) - r*eps
1734 df = -2.0_dp*oorootpi*omega*exp(-qr**2) - eps
1735 END SUBROUTINE eval_transc_func
1736 END SUBROUTINE erfc_cutoff
1737
1738! **************************************************************************************************
1739!> \brief Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd
1740!> \param matrix ...
1741!> \param eigenvectors ...
1742!> \param eigenvalues ...
1743!> \author A. Bussy
1744! **************************************************************************************************
1745 SUBROUTINE complex_diag(matrix, eigenvectors, eigenvalues)
1746 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: matrix, eigenvectors
1747 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1748
1749 COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1750 INTEGER :: info, liwork, lrwork, lwork, n
1751 INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1752 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: rwork
1753
1754 eigenvectors(:, :) = matrix(:, :)
1755
1756 n = SIZE(matrix, 1)
1757 ALLOCATE (iwork(1), rwork(1), work(1))
1758
1759 ! work space query
1760 lwork = -1
1761 lrwork = -1
1762 liwork = -1
1763
1764 CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1765
1766 lwork = ceiling(real(work(1), kind=dp))
1767 lrwork = ceiling(real(rwork(1), kind=dp))
1768 liwork = iwork(1)
1769
1770 DEALLOCATE (iwork, rwork, work)
1771 ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
1772
1773 ! diagonalization proper
1774 CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1775
1776 DEALLOCATE (iwork, rwork, work)
1777 IF (info /= 0) &
1778 cpabort("Diagonalisation of a complex matrix failed")
1779
1780 END SUBROUTINE complex_diag
1781
1782END 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:195
mathlib::det_3x3
Definition mathlib.F:65
mathlib::invert_matrix
Definition mathlib.F:69
mathlib::set_diag
Definition mathlib.F:73
mathlib::unit_matrix
Definition mathlib.F:81
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:949
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:1479
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:205
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::invmat_symm
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
Definition mathlib.F:574
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:229
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: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:175
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:285
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::complex_diag
subroutine, public complex_diag(matrix, eigenvectors, eigenvalues)
Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd.
Definition mathlib.F:1746
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:372
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:1507
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:1689
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:150
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:1404
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:99
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:254