(git:08825b9)
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!> \param exp_mx optional precomputed EXP(-x)
1525!> \return ...
1526!> \par History
1527!> 05.2007 Created
1528!> \author Manuel Guidon (adapted from Numerical recipies)
1529! **************************************************************************************************
1530 ELEMENTAL IMPURE FUNCTION expint(n, x, exp_mx)
1531 INTEGER, INTENT(IN) :: n
1532 REAL(dp), INTENT(IN) :: x
1533 REAL(dp), INTENT(IN), OPTIONAL :: exp_mx
1534 REAL(dp) :: expint
1535
1536 INTEGER, PARAMETER :: maxit = 100
1537 REAL(dp), PARAMETER :: cf_rescale_inv = 1.0e-100_dp, cf_rescale_max = 1.0e100_dp, &
1538 eps = 6.e-14_dp, euler = 0.5772156649015328606065120_dp
1539
1540 INTEGER :: i, ii, nm1
1541 REAL(dp) :: a, a0, a1, a2, b, b0, b1, b2, del, fact, &
1542 my_exp_mx, psi
1543
1544 nm1 = n - 1
1545
1546 IF (n < 0 .OR. x < 0.0_dp .OR. (x == 0.0_dp .AND. (n == 0 .OR. n == 1))) THEN
1547 cpabort("Invalid argument")
1548 ELSE IF (n == 0) THEN !Special case.
1549 IF (PRESENT(exp_mx)) THEN
1550 expint = exp_mx/x
1551 ELSE
1552 expint = exp(-x)/x
1553 END IF
1554 ELSE IF (x == 0.0_dp) THEN !Another special case.
1555 expint = 1.0_dp/nm1
1556 ELSE IF (x > 750.0_dp) THEN !E_n(x) <= exp(-x)/x underflows to zero in IEEE-754 dp.
1557 expint = 0.0_dp
1558 ELSE IF (x > 1.0_dp) THEN !3-term continued-fraction recurrence.
1559 IF (PRESENT(exp_mx)) THEN
1560 my_exp_mx = exp_mx
1561 ELSE
1562 my_exp_mx = exp(-x)
1563 END IF
1564 b = x + n
1565 a0 = 0.0_dp
1566 b0 = 1.0_dp
1567 a1 = 1.0_dp
1568 b1 = b
1569 DO i = 1, maxit
1570 a = -i*(nm1 + i)
1571 b = b + 2.0_dp
1572 a2 = b*a1 + a*a0
1573 b2 = b*b1 + a*b0
1574 ! Rescale 3-term continued-fraction coefficients to prevent overflow
1575 IF (abs(b2) > cf_rescale_max) THEN
1576 a2 = a2*cf_rescale_inv
1577 b2 = b2*cf_rescale_inv
1578 a1 = a1*cf_rescale_inv
1579 b1 = b1*cf_rescale_inv
1580 END IF
1581 IF (abs(a2*b1 - a1*b2) < eps*abs(a2*b1) .OR. i == maxit) THEN
1582 IF (abs(a2*b1 - a1*b2) >= eps*abs(a2*b1)) THEN
1583 cpabort("continued fraction failed in expint")
1584 END IF
1585 expint = (a2/b2)*my_exp_mx
1586 RETURN
1587 END IF
1588 a0 = a1
1589 b0 = b1
1590 a1 = a2
1591 b1 = b2
1592 END DO
1593 ELSE !Evaluate series.
1594 IF (nm1 /= 0) THEN !Set first term.
1595 expint = 1.0_dp/nm1
1596 ELSE
1597 expint = -log(x) - euler
1598 END IF
1599 fact = 1.0_dp
1600 DO i = 1, maxit
1601 fact = -fact*x/i
1602 IF (i /= nm1) THEN
1603 del = -fact/(i - nm1)
1604 ELSE
1605 psi = -euler !Compute I(n).
1606 DO ii = 1, nm1
1607 psi = psi + 1.0_dp/ii
1608 END DO
1609 del = fact*(-log(x) + psi)
1610 END IF
1611 expint = expint + del
1612 IF (abs(del) < abs(expint)*eps) RETURN
1613 END DO
1614 cpabort("series failed in expint")
1615 END IF
1616 END FUNCTION expint
1617
1618! **************************************************************************************************
1619!> \brief Jacobi matrix diagonalization. The eigenvalues are returned in
1620!> vector d and the eigenvectors are returned in matrix v in ascending
1621!> order.
1622!>
1623!> \param a ...
1624!> \param d ...
1625!> \param v ...
1626!> \par History
1627!> - Creation (20.11.98, Matthias Krack)
1628! **************************************************************************************************
1629 SUBROUTINE jacobi(a, d, v)
1630 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1631 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: d
1632 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: v
1633
1634 INTEGER :: n
1635
1636 n = SIZE(d(:))
1637
1638 ! Diagonalize matrix a
1639 CALL diag(n, a, d, v)
1640
1641 ! Sort eigenvalues and eigenvector in ascending order
1642 CALL eigsrt(n, d, v)
1643
1644 END SUBROUTINE jacobi
1645
1646! **************************************************************************************************
1647!> \brief Diagonalize matrix a. The eigenvalues are returned in vector d
1648!> and the eigenvectors are returned in matrix v.
1649!>
1650!> \param n matrix/vector extent (problem size)
1651!> \param a matrix to be diagonalised
1652!> \param d vector of eigenvalues
1653!> \param v matrix of eigenvectors
1654!> \par History
1655!> - Creation (20.11.98, Matthias Krack)
1656! **************************************************************************************************
1657 SUBROUTINE diag(n, a, d, v)
1658 INTEGER, INTENT(IN) :: n
1659 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1660 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: d
1661 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: v
1662
1663 CHARACTER(len=*), PARAMETER :: routinen = 'diag'
1664 INTEGER, PARAMETER :: maxit = 100
1665 REAL(kind=dp), PARAMETER :: a_eps = 1.0e-10_dp, d_eps = 1.0e-3_dp
1666
1667 CHARACTER(LEN=default_string_length) :: max_iter
1668 INTEGER :: handle, i, ip, iq
1669 REAL(kind=dp) :: a_max, apq, c, d_min, dip, diq, g, h, s, &
1670 t, tau, theta, tresh
1671 REAL(kind=dp), DIMENSION(n) :: b, z
1672
1673 CALL timeset(routinen, handle)
1674
1675 a_max = 0.0_dp
1676 DO ip = 1, n - 1
1677 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1678 b(ip) = a(ip, ip) ! get_diag(a)
1679 END DO
1680 b(n) = a(n, n)
1681
1682 CALL unit_matrix(v)
1683
1684 DO i = 1, maxit
1685 d = b
1686 d_min = max(d_eps, minval(abs(b)))
1687 IF (a_max < a_eps*d_min) THEN
1688 CALL timestop(handle)
1689 RETURN
1690 END IF
1691 tresh = merge(a_max, 0.0_dp, (i < 4))
1692 z = 0.0_dp
1693 DO ip = 1, n - 1
1694 DO iq = ip + 1, n
1695 dip = d(ip)
1696 diq = d(iq)
1697 apq = a(ip, iq)
1698 g = 100.0_dp*abs(apq)
1699 IF (tresh < abs(apq)) THEN
1700 h = diq - dip
1701 IF ((abs(h) + g) /= abs(h)) THEN
1702 theta = 0.5_dp*h/apq
1703 t = 1.0_dp/(abs(theta) + sqrt(1.0_dp + theta**2))
1704 IF (theta < 0.0_dp) t = -t
1705 ELSE
1706 t = apq/h
1707 END IF
1708 c = 1.0_dp/sqrt(1.0_dp + t**2)
1709 s = t*c
1710 tau = s/(1.0_dp + c)
1711 h = t*apq
1712 z(ip) = z(ip) - h
1713 z(iq) = z(iq) + h
1714 d(ip) = dip - h
1715 d(iq) = diq + h
1716 a(ip, iq) = 0.0_dp
1717 CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1718 CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1719 CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1720 CALL jrotate(v(:, ip), v(:, iq), s, tau)
1721 ELSE IF ((4 < i) .AND. &
1722 ((abs(dip) + g) == abs(dip)) .AND. &
1723 ((abs(diq) + g) == abs(diq))) THEN
1724 a(ip, iq) = 0.0_dp
1725 END IF
1726 END DO
1727 END DO
1728 b = b + z
1729 a_max = 0.0_dp
1730 DO ip = 1, n - 1
1731 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1732 END DO
1733 END DO
1734 WRITE (max_iter, *) maxit
1735 CALL cp_warn(__location__, &
1736 "Jacobi matrix diagonalization does not converge after "// &
1737 trim(adjustl(max_iter))//" iterations; stop here anyway")
1738
1739 CALL timestop(handle)
1740
1741 END SUBROUTINE diag
1742
1743! **************************************************************************************************
1744!> \brief Perform a Jacobi rotation of the vectors a and b.
1745!>
1746!> \param a ...
1747!> \param b ...
1748!> \param ss ...
1749!> \param tt ...
1750!> \par History
1751!> - Creation (20.11.98, Matthias Krack)
1752! **************************************************************************************************
1753 PURE SUBROUTINE jrotate(a, b, ss, tt)
1754 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1755 REAL(kind=dp), INTENT(IN) :: ss, tt
1756
1757 REAL(kind=dp) :: u, v
1758
1759 u = 1.0_dp - ss*tt
1760 v = ss/u
1761
1762 a = a*u - b*ss
1763 b = b*(u + ss*v) + a*v
1764
1765 END SUBROUTINE jrotate
1766
1767! **************************************************************************************************
1768!> \brief Sort the values in vector d in ascending order and swap the
1769!> corresponding columns of matrix v.
1770!>
1771!> \param n ...
1772!> \param d ...
1773!> \param v ...
1774!> \par History
1775!> - Creation (20.11.98, Matthias Krack)
1776! **************************************************************************************************
1777 SUBROUTINE eigsrt(n, d, v)
1778 INTEGER, INTENT(IN) :: n
1779 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: d
1780 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: v
1781
1782 INTEGER :: i, j
1783
1784 DO i = 1, n - 1
1785 j = sum(minloc(d(i:n))) + i - 1
1786 IF (j /= i) THEN
1787 CALL swap(d(i), d(j))
1788 CALL swap(v(:, i), v(:, j))
1789 END IF
1790 END DO
1791
1792 END SUBROUTINE eigsrt
1793
1794! **************************************************************************
1795!> \brief Swap two scalars
1796!>
1797!> \param a ...
1798!> \param b ...
1799!> \par History
1800!> - Creation (20.11.98, Matthias Krack)
1801! **************************************************************************************************
1802 ELEMENTAL SUBROUTINE swap_scalar(a, b)
1803 REAL(kind=dp), INTENT(INOUT) :: a, b
1804
1805 REAL(kind=dp) :: c
1806
1807 c = a
1808 a = b
1809 b = c
1810
1811 END SUBROUTINE swap_scalar
1812
1813! **************************************************************************
1814!> \brief Swap two vectors
1815!>
1816!> \param a ...
1817!> \param b ...
1818!> \par History
1819!> - Creation (20.11.98, Matthias Krack)
1820! **************************************************************************************************
1821 SUBROUTINE swap_vector(a, b)
1822 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1823
1824 INTEGER :: i, n
1825 REAL(kind=dp) :: c
1826
1827 n = SIZE(a)
1828
1829 IF (n /= SIZE(b)) THEN
1830 cpabort("Check the array bounds of the parameters")
1831 END IF
1832
1833 DO i = 1, n
1834 c = a(i)
1835 a(i) = b(i)
1836 b(i) = c
1837 END DO
1838
1839 END SUBROUTINE swap_vector
1840
1841! **************************************************************************************************
1842!> \brief - compute a truncation radius for the shortrange operator
1843!> \param eps target accuracy!> \param omg screening parameter
1844!> \param omg ...
1845!> \param r_cutoff cutoff radius
1846!> \par History
1847!> 10.2012 created [Hossein Banihashemian]
1848!> 05.2019 moved here from hfx_types (A. Bussy)
1849!> \author Hossein Banihashemian
1850! **************************************************************************************************
1851 SUBROUTINE erfc_cutoff(eps, omg, r_cutoff)
1852 REAL(dp), INTENT(in) :: eps, omg
1853 REAL(dp), INTENT(out) :: r_cutoff
1854
1855 CHARACTER(LEN=*), PARAMETER :: routinen = 'erfc_cutoff'
1856
1857 REAL(dp), PARAMETER :: abstol = 1e-10_dp, soltol = 1e-16_dp
1858 REAL(dp) :: r0, f0, fprime0, delta_r
1859 INTEGER :: iter, handle
1860 INTEGER, PARAMETER :: itermax = 1000
1861
1862 CALL timeset(routinen, handle)
1863
1864 ! initial guess assuming that we are in the asymptotic regime of the erf, and the solution is about 10.
1865 r0 = sqrt(-log(eps*omg*10**2))/omg
1866 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1867
1868 DO iter = 1, itermax
1869 delta_r = f0/fprime0
1870 r0 = r0 - delta_r
1871 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1872 IF (abs(delta_r) < abstol .OR. abs(f0) < soltol) EXIT
1873 END DO
1874 cpassert(iter <= itermax)
1875 r_cutoff = r0
1876
1877 CALL timestop(handle)
1878 CONTAINS
1879! **************************************************************************************************
1880!> \brief ...
1881!> \param r ...
1882!> \param eps ...
1883!> \param omega ...
1884!> \param fn ...
1885!> \param df ...
1886! **************************************************************************************************
1887 ELEMENTAL SUBROUTINE eval_transc_func(r, eps, omega, fn, df)
1888 REAL(dp), INTENT(in) :: r, eps, omega
1889 REAL(dp), INTENT(out) :: fn, df
1890
1891 REAL(dp) :: qr
1892
1893 qr = omega*r
1894 fn = erfc(qr) - r*eps
1895 df = -2.0_dp*oorootpi*omega*exp(-qr**2) - eps
1896 END SUBROUTINE eval_transc_func
1897 END SUBROUTINE erfc_cutoff
1898
1899! **************************************************************************************************
1900!> \brief Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd
1901!> \param matrix Hermitian matrix is preserved
1902!> \param eigenvectors ...
1903!> \param eigenvalues ...
1904!> \author A. Bussy
1905! **************************************************************************************************
1906 SUBROUTINE diag_complex(matrix, eigenvectors, eigenvalues)
1907 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
1908 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
1909 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1910
1911 CHARACTER(len=*), PARAMETER :: routinen = 'diag_complex'
1912
1913 COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1914 INTEGER :: handle, info, liwork, lrwork, lwork, n
1915 INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1916 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: rwork
1917
1918 CALL timeset(routinen, handle)
1919
1920 IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) cpabort("Expected square matrix")
1921 ! IF (MAXVAL(ABS(matrix - CONJG(TRANSPOSE(matrix)))) > 1e-14_dp) CPABORT("Expected hermitian matrix")
1922
1923 n = SIZE(matrix, 1)
1924 ALLOCATE (iwork(1), rwork(1), work(1))
1925
1926 ! work space query
1927 lwork = -1
1928 lrwork = -1
1929 liwork = -1
1930
1931 CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1932
1933 lwork = ceiling(real(work(1), kind=dp))
1934 lrwork = ceiling(rwork(1))
1935 liwork = iwork(1)
1936
1937 DEALLOCATE (iwork, rwork, work)
1938 ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
1939 eigenvectors(:, :) = matrix(:, :)
1940
1941 ! final diagonalization
1942 CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1943
1944 DEALLOCATE (iwork, rwork, work)
1945
1946 IF (info /= 0) cpabort("Diagonalisation of a complex matrix failed")
1947
1948 CALL timestop(handle)
1949
1950 END SUBROUTINE diag_complex
1951
1952! **************************************************************************************************
1953!> \brief Helper routine for diagonalizing anti symmetric matrices
1954!> \param matrix ...
1955!> \param evecs ...
1956!> \param evals ...
1957! **************************************************************************************************
1958 SUBROUTINE diag_antisym(matrix, evecs, evals)
1959 REAL(dp), DIMENSION(:, :) :: matrix
1960 COMPLEX(dp), DIMENSION(:, :) :: evecs
1961 COMPLEX(dp), DIMENSION(:) :: evals
1962
1963 COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: matrix_c
1964 INTEGER :: n
1965 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1966
1967 IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) cpabort("Expected square matrix")
1968 ! IF (MAXVAL(ABS(matrix + TRANSPOSE(matrix))) > 1e-14_dp) CPABORT("Expected anti-symmetric matrix")
1969
1970 n = SIZE(matrix, 1)
1971 ALLOCATE (matrix_c(n, n), eigenvalues(n))
1972
1973 matrix_c(:, :) = cmplx(0.0_dp, -matrix, kind=dp)
1974 CALL diag_complex(matrix_c, evecs, eigenvalues)
1975 evals = cmplx(0.0_dp, eigenvalues, kind=dp)
1976
1977 DEALLOCATE (matrix_c, eigenvalues)
1978 END SUBROUTINE diag_antisym
1979! **************************************************************************************************
1980!> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged
1981!> \param A_in Input matrix 1
1982!> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1983!> \param B_in Input matrix 2
1984!> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1985!> \param C_out Output matrix
1986!> \par History
1987!> 11.2025 created [Stepan Marek]
1988! **************************************************************************************************
1989 SUBROUTINE zgemm_square_2(A_in, A_trans, B_in, B_trans, C_out)
1990 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1991 CHARACTER, INTENT(IN) :: A_trans
1992 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1993 CHARACTER, INTENT(IN) :: B_trans
1994 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: C_out
1995
1996 CHARACTER(len=*), PARAMETER :: routineN = 'zgemm_square_2'
1997
1998 INTEGER :: handle, n
1999
2000 n = SIZE(a_in, 1)
2001 IF (n /= SIZE(a_in, 2)) cpabort("Non-square array 1 (A).")
2002 IF (n /= SIZE(b_in, 1)) cpabort("Incompatible (rows) array 2 (B).")
2003 IF (n /= SIZE(b_in, 2)) cpabort("Non-square array 2 (B).")
2004 IF (n /= SIZE(c_out, 1)) cpabort("Incompatible (rows) result array 3 (C).")
2005 IF (n /= SIZE(c_out, 2)) cpabort("Incompatible (cols) result array 3 (C).")
2006 IF (.NOT. (a_trans == 'N' .OR. a_trans == 'n' .OR. &
2007 a_trans == 'T' .OR. a_trans == 't' .OR. &
2008 a_trans == 'C' .OR. a_trans == 'c')) THEN
2009 cpabort("Unknown transpose character for array 1 (A).")
2010 END IF
2011 IF (.NOT. (b_trans == 'N' .OR. b_trans == 'n' .OR. &
2012 b_trans == 'T' .OR. b_trans == 't' .OR. &
2013 b_trans == 'C' .OR. b_trans == 'c')) THEN
2014 cpabort("Unknown transpose character for array 2 (B).")
2015 END IF
2016
2017 CALL timeset(routinen, handle)
2018
2019 CALL zgemm(a_trans, b_trans, n, n, n, z_one, a_in, n, b_in, n, z_zero, c_out, n)
2020
2021 CALL timestop(handle)
2022
2023 END SUBROUTINE zgemm_square_2
2024! **************************************************************************************************
2025!> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, real matrices
2026!> \param A_in Input matrix 1
2027!> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
2028!> \param B_in Input matrix 2
2029!> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
2030!> \param C_out Output matrix
2031!> \par History
2032!> 11.2025 created [Stepan Marek]
2033! **************************************************************************************************
2034 SUBROUTINE dgemm_square_2(A_in, A_trans, B_in, B_trans, C_out)
2035 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a_in
2036 CHARACTER, INTENT(IN) :: A_trans
2037 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: b_in
2038 CHARACTER, INTENT(IN) :: B_trans
2039 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: c_out
2040
2041 CHARACTER(len=*), PARAMETER :: routineN = 'dgemm_square_2'
2042
2043 INTEGER :: handle, n
2044
2045 n = SIZE(a_in, 1)
2046 IF (n /= SIZE(a_in, 2)) cpabort("Non-square array 1 (A).")
2047 IF (n /= SIZE(b_in, 1)) cpabort("Incompatible (rows) array 2 (B).")
2048 IF (n /= SIZE(b_in, 2)) cpabort("Non-square array 2 (B).")
2049 IF (n /= SIZE(c_out, 1)) cpabort("Incompatible (rows) result array 3 (C).")
2050 IF (n /= SIZE(c_out, 2)) cpabort("Incompatible (cols) result array 3 (C).")
2051 IF (.NOT. (a_trans == 'N' .OR. a_trans == 'n' .OR. &
2052 a_trans == 'T' .OR. a_trans == 't' .OR. &
2053 a_trans == 'C' .OR. a_trans == 'c')) THEN
2054 cpabort("Unknown transpose character for array 1 (A).")
2055 END IF
2056 IF (.NOT. (b_trans == 'N' .OR. b_trans == 'n' .OR. &
2057 b_trans == 'T' .OR. b_trans == 't' .OR. &
2058 b_trans == 'C' .OR. b_trans == 'c')) THEN
2059 cpabort("Unknown transpose character for array 2 (B).")
2060 END IF
2061
2062 CALL timeset(routinen, handle)
2063
2064 CALL dgemm(a_trans, b_trans, n, n, n, 1.0_dp, a_in, n, b_in, n, 0.0_dp, c_out, n)
2065
2066 CALL timestop(handle)
2067
2068 END SUBROUTINE dgemm_square_2
2069! **************************************************************************************************
2070!> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, for 3 matrices
2071!> \param A_in Input matrix 1
2072!> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
2073!> \param B_in Input matrix 2
2074!> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
2075!> \param C_in Input matrix 3
2076!> \param C_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 3
2077!> \param D_out Output matrix
2078!> \par History
2079!> 11.2025 created [Stepan Marek]
2080! **************************************************************************************************
2081 SUBROUTINE zgemm_square_3(A_in, A_trans, B_in, B_trans, C_in, C_trans, D_out)
2082 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
2083 CHARACTER, INTENT(IN) :: A_trans
2084 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
2085 CHARACTER, INTENT(IN) :: B_trans
2086 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: C_in
2087 CHARACTER, INTENT(IN) :: C_trans
2088 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: D_out
2089
2090 CHARACTER(len=*), PARAMETER :: routineN = 'zgemm_square_3'
2091
2092 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
2093 INTEGER :: handle, n
2094
2095 n = SIZE(a_in, 1)
2096 IF (n /= SIZE(a_in, 2)) cpabort("Non-square array 1 (A).")
2097 IF (n /= SIZE(b_in, 1)) cpabort("Incompatible (rows) array 2 (B).")
2098 IF (n /= SIZE(b_in, 2)) cpabort("Non-square array 2 (B).")
2099 IF (n /= SIZE(c_in, 1)) cpabort("Incompatible (rows) array 3 (C).")
2100 IF (n /= SIZE(c_in, 2)) cpabort("Non-square array 3 (C).")
2101 IF (n /= SIZE(d_out, 1)) cpabort("Incompatible (rows) result array 4 (D).")
2102 IF (n /= SIZE(d_out, 2)) cpabort("Incompatible (cols) result array 4 (D).")
2103 IF (.NOT. (a_trans == 'N' .OR. a_trans == 'n' .OR. &
2104 a_trans == 'T' .OR. a_trans == 't' .OR. &
2105 a_trans == 'C' .OR. a_trans == 'c')) THEN
2106 cpabort("Unknown transpose character for array 1 (A).")
2107 END IF
2108 IF (.NOT. (b_trans == 'N' .OR. b_trans == 'n' .OR. &
2109 b_trans == 'T' .OR. b_trans == 't' .OR. &
2110 b_trans == 'C' .OR. b_trans == 'c')) THEN
2111 cpabort("Unknown transpose character for array 2 (B).")
2112 END IF
2113 IF (.NOT. (c_trans == 'N' .OR. c_trans == 'n' .OR. &
2114 c_trans == 'T' .OR. c_trans == 't' .OR. &
2115 c_trans == 'C' .OR. c_trans == 'c')) THEN
2116 cpabort("Unknown transpose character for array 3 (C).")
2117 END IF
2118
2119 CALL timeset(routinen, handle)
2120
2121 ALLOCATE (work(n, n), source=z_zero)
2122
2123 CALL zgemm(a_trans, b_trans, n, n, n, z_one, a_in, n, b_in, n, z_zero, work, n)
2124 CALL zgemm('N', c_trans, n, n, n, z_one, work, n, c_in, n, z_zero, d_out, n)
2125
2126 DEALLOCATE (work)
2127
2128 CALL timestop(handle)
2129
2130 END SUBROUTINE zgemm_square_3
2131! **************************************************************************************************
2132!> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, for 3 matrices
2133!> \param A_in Input matrix 1
2134!> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
2135!> \param B_in Input matrix 2
2136!> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
2137!> \param C_in Input matrix 3
2138!> \param C_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 3
2139!> \param D_out Output matrix
2140!> \par History
2141!> 11.2025 created [Stepan Marek]
2142! **************************************************************************************************
2143 SUBROUTINE dgemm_square_3(A_in, A_trans, B_in, B_trans, C_in, C_trans, D_out)
2144 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a_in
2145 CHARACTER, INTENT(IN) :: A_trans
2146 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: b_in
2147 CHARACTER, INTENT(IN) :: B_trans
2148 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: c_in
2149 CHARACTER, INTENT(IN) :: C_trans
2150 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: d_out
2151
2152 CHARACTER(len=*), PARAMETER :: routineN = 'dgemm_square_3'
2153
2154 INTEGER :: handle, n
2155 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
2156
2157 n = SIZE(a_in, 1)
2158 IF (n /= SIZE(a_in, 2)) cpabort("Non-square array 1 (A).")
2159 IF (n /= SIZE(b_in, 1)) cpabort("Incompatible (rows) array 2 (B).")
2160 IF (n /= SIZE(b_in, 2)) cpabort("Non-square array 2 (B).")
2161 IF (n /= SIZE(c_in, 1)) cpabort("Incompatible (rows) array 3 (C).")
2162 IF (n /= SIZE(c_in, 2)) cpabort("Non-square array 3 (C).")
2163 IF (n /= SIZE(d_out, 1)) cpabort("Incompatible (rows) result array 4 (D).")
2164 IF (n /= SIZE(d_out, 2)) cpabort("Incompatible (cols) result array 4 (D).")
2165 IF (.NOT. (a_trans == 'N' .OR. a_trans == 'n' .OR. &
2166 a_trans == 'T' .OR. a_trans == 't' .OR. &
2167 a_trans == 'C' .OR. a_trans == 'c')) THEN
2168 cpabort("Unknown transpose character for array 1 (A).")
2169 END IF
2170 IF (.NOT. (b_trans == 'N' .OR. b_trans == 'n' .OR. &
2171 b_trans == 'T' .OR. b_trans == 't' .OR. &
2172 b_trans == 'C' .OR. b_trans == 'c')) THEN
2173 cpabort("Unknown transpose character for array 2 (B).")
2174 END IF
2175 IF (.NOT. (c_trans == 'N' .OR. c_trans == 'n' .OR. &
2176 c_trans == 'T' .OR. c_trans == 't' .OR. &
2177 c_trans == 'C' .OR. c_trans == 'c')) THEN
2178 cpabort("Unknown transpose character for array 3 (C).")
2179 END IF
2180
2181 CALL timeset(routinen, handle)
2182
2183 ALLOCATE (work(n, n), source=0.0_dp)
2184
2185 CALL dgemm(a_trans, b_trans, n, n, n, 1.0_dp, a_in, n, b_in, n, 0.0_dp, work, n)
2186 CALL dgemm('N', c_trans, n, n, n, 1.0_dp, work, n, c_in, n, 0.0_dp, d_out, n)
2187
2188 DEALLOCATE (work)
2189
2190 CALL timestop(handle)
2191
2192 END SUBROUTINE dgemm_square_3
2193
2194! **************************************************************************************************
2195!> \brief Solve the generalized eigenvalue equation for complex matrices A*v = B*v*λ
2196!> \param A_in ...
2197!> \param B_in ...
2198!> \param eigenvalues ...
2199!> \param eigenvectors ...
2200!> \author Shridhar Shanbhag
2201! **************************************************************************************************
2202 SUBROUTINE geeig_right(A_in, B_in, eigenvalues, eigenvectors)
2203 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: a_in, b_in
2204 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
2205 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
2206
2207 CHARACTER(len=*), PARAMETER :: routinen = 'geeig_right'
2208 COMPLEX(KIND=dp), PARAMETER :: cone = cmplx(1.0_dp, 0.0_dp, kind=dp), &
2209 czero = cmplx(0.0_dp, 0.0_dp, kind=dp)
2210
2211 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
2212 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, work
2213 INTEGER :: handle, i, icol, irow, lda, ldb, ldc, &
2214 nao, nc, ncol, nmo, nx
2215 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
2216
2217 CALL timeset(routinen, handle)
2218
2219 ! Test sizes
2220 nao = SIZE(a_in, 1)
2221 nmo = SIZE(eigenvalues)
2222 ALLOCATE (evals(nao), cevals(nao))
2223 ALLOCATE (work(nao, nao), b(nao, nao), a(nao, nao))
2224 a(:, :) = a_in(:, :)
2225 b(:, :) = b_in(:, :)
2226
2227 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
2228 b = -b
2229 CALL diag_complex(b, work, evals)
2230 evals(:) = -evals(:)
2231 nc = nao
2232 DO i = 1, nao
2233 IF (evals(i) < -1.0_dp) THEN
2234 nc = i - 1
2235 EXIT
2236 END IF
2237 END DO
2238 cpassert(nc /= 0)
2239
2240 IF (nc /= nao) THEN
2241 IF (nc < nmo) THEN
2242 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
2243 ncol = nmo - nc
2244 CALL zcopy(ncol*nao, work(1, nc + 1), 1, eigenvectors(1, nc + 1), 1)
2245 END IF
2246 ! Set NULL space in eigenvector matrix of S to zero
2247 DO icol = nc + 1, nao
2248 DO irow = 1, nao
2249 work(irow, icol) = czero
2250 END DO
2251 END DO
2252 ! Set small eigenvalues to a dummy save value
2253 evals(nc + 1:nao) = 1.0_dp
2254 END IF
2255 ! Calculate U*s**(-1/2)
2256 cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
2257 DO i = 1, min(SIZE(work, 2), SIZE(cevals))
2258 work(1:min(SIZE(work, 2), SIZE(cevals)), i) = cevals(i)*work(1:min(SIZE(work, 2), SIZE(cevals)), i)
2259 END DO
2260 ! Reduce to get U^(-C) * H * U^(-1)
2261 CALL gemm_square(work, 'C', a, 'N', b)
2262 CALL gemm_square(b, 'N', work, 'N', a)
2263 IF (nc /= nao) THEN
2264 ! set diagonal values to save large value
2265 DO icol = nc + 1, nao
2266 a(icol, icol) = 10000*cone
2267 END DO
2268 END IF
2269 ! Diagonalize
2270 CALL diag_complex(a, b, evals)
2271 eigenvalues(1:nmo) = evals(1:nmo)
2272 nx = min(nc, nmo)
2273 ! Restore vectors C = U^(-1) * C*
2274 lda = SIZE(work, 1)
2275 ldb = SIZE(b, 1)
2276 ldc = SIZE(eigenvectors, 1)
2277 CALL zgemm("N", "N", nao, nx, nc, cone, work, &
2278 lda, b, ldb, czero, eigenvectors, ldc)
2279
2280 DEALLOCATE (evals, cevals, work, b, a)
2281
2282 CALL timestop(handle)
2283 END SUBROUTINE geeig_right
2284
2285END 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:1630
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::expint
elemental impure real(dp) function, public expint(n, x, exp_mx)
computes the exponential integral En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,...
Definition mathlib.F:1531
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:1907
mathlib::diag_antisym
subroutine, public diag_antisym(matrix, evecs, evals)
Helper routine for diagonalizing anti symmetric matrices.
Definition mathlib.F:1959
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:1658
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:2203
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:1852
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::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