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