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