(git:339af6c)
Loading...
Searching...
No Matches
mathlib.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Collection of simple mathematical functions and subroutines
10!> \par History
11!> FUNCTION angle updated and FUNCTION dihedral angle added; cleaned
12!> (13.03.2004,MK)
13!> \author MK (15.11.1998)
14! **************************************************************************************************
15MODULE mathlib
16
17 USE kinds, ONLY: default_string_length,&
18 dp
19 USE mathconstants, ONLY: euler,&
20 fac,&
21 oorootpi,&
22 z_one,&
23 z_zero
24#include "../base/base_uses.f90"
25
26 IMPLICIT NONE
27
28 PRIVATE
29
30 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mathlib'
31 REAL(KIND=dp), PARAMETER :: eps_geo = 1.0e-6_dp
32
33 ! Public subroutines
34
35 PUBLIC :: build_rotmat, &
36 jacobi, &
37 diamat_all, &
38 invmat, &
39 invmat_symm, &
40 invert_matrix, &
41 get_pseudo_inverse_svd, &
42 get_pseudo_inverse_diag, &
43 symmetrize_matrix, &
44 unit_matrix, diag, &
45 erfc_cutoff, &
46 diag_antisym, &
47 diag_complex, &
48 gemm_square, &
49 geeig_right
50
51 ! Public functions
52
53 PUBLIC :: angle, &
54 binomial, &
55 binomial_gen, &
56 multinomial, &
57 det_3x3, &
58 dihedral_angle, &
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 = sqrt(dot_product(a, a))
191 length_of_b = sqrt(dot_product(b, 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)) &
1001 determinant = determinant*sig(i)
1002 sig_plus(i, i) = 1._dp/sig(i)
1003 ELSE
1004 sig_plus(i, i) = 0.0_dp
1005 END IF
1006 END DO
1007
1008 ! build pseudoinverse: V*sig_plus*UT
1009 CALL dgemm("N", "T", n, n, n, 1._dp, sig_plus, n, u, n, 0._dp, temp_mat, n)
1010 CALL dgemm("T", "N", n, n, n, 1._dp, vt, n, temp_mat, n, 0._dp, a_pinverse, n)
1011
1012 DEALLOCATE (u, vt, sig, iwork, work, sig_plus, temp_mat)
1013
1014 CALL timestop(handle)
1015
1016 END SUBROUTINE get_pseudo_inverse_svd
1017
1018! **************************************************************************************************
1019!> \brief returns the pseudoinverse of a real, symmetric and positive definite
1020!> matrix using diagonalization.
1021!> \param a matrix a
1022!> \param a_pinverse pseudoinverse of matrix a
1023!> \param rskip parameter for setting small eigenvalues to zero
1024!> \author Dorothea Golze [02.2015]
1025! **************************************************************************************************
1026 SUBROUTINE get_pseudo_inverse_diag(a, a_pinverse, rskip)
1027 REAL(kind=dp), DIMENSION(:, :) :: a, a_pinverse
1028 REAL(kind=dp), INTENT(IN) :: rskip
1029
1030 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_pseudo_inverse_diag'
1031
1032 INTEGER :: handle, i, info, lwork, n
1033 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
1034 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: dinv, temp_mat
1035
1036 CALL timeset(routinen, handle)
1037
1038 info = 0
1039 n = SIZE(a, 1)
1040 ALLOCATE (dinv(n, n), eig(n), work(1), temp_mat(n, n))
1041 dinv(:, :) = 0.0_dp
1042 eig(:) = 0.0_dp
1043 work(:) = 0.0_dp
1044 temp_mat = 0.0_dp
1045
1046 ! work size query
1047 lwork = -1
1048 CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1049 IF (info /= 0) THEN
1050 cpabort("ERROR in DSYEV: Could not retrieve work array sizes")
1051 END IF
1052 lwork = int(work(1))
1053 DEALLOCATE (work)
1054 ALLOCATE (work(lwork))
1055 work = 0.0_dp
1056
1057 ! get eigenvalues and eigenvectors
1058 CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1059
1060 IF (info /= 0) THEN
1061 cpabort("Matrix diagonalization failed")
1062 END IF
1063
1064 ! set eigenvalues that are too small to zero
1065 DO i = 1, n
1066 IF (eig(i) > rskip*maxval(eig)) THEN
1067 dinv(i, i) = 1.0_dp/eig(i)
1068 ELSE
1069 dinv(i, i) = 0._dp
1070 END IF
1071 END DO
1072
1073 ! build pseudoinverse: U*dinv*UT
1074 CALL dgemm("N", "T", n, n, n, 1._dp, dinv, n, a, n, 0._dp, temp_mat, n)
1075 CALL dgemm("N", "N", n, n, n, 1._dp, a, n, temp_mat, n, 0._dp, a_pinverse, n)
1076
1077 DEALLOCATE (eig, work, dinv, temp_mat)
1078
1079 CALL timestop(handle)
1080
1081 END SUBROUTINE get_pseudo_inverse_diag
1082
1083! **************************************************************************************************
1084!> \brief Reflection of the vector a through a mirror plane defined by the
1085!> normal vector b. The reflected vector a is stored in a_mirror.
1086!> \param a ...
1087!> \param b ...
1088!> \return ...
1089!> \date 16.10.1998
1090!> \author MK
1091!> \version 1.0
1092! **************************************************************************************************
1093 PURE FUNCTION reflect_vector(a, b) RESULT(a_mirror)
1094 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a, b
1095 REAL(kind=dp), DIMENSION(3) :: a_mirror
1096
1097 REAL(kind=dp) :: length_of_b, scapro
1098 REAL(kind=dp), DIMENSION(3) :: d
1099
1100 length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1101
1102 IF (length_of_b > eps_geo) THEN
1103
1104 d(:) = b(:)/length_of_b
1105
1106 ! Calculate the mirror image a_mirror of the vector a
1107 scapro = a(1)*d(1) + a(2)*d(2) + a(3)*d(3)
1108
1109 a_mirror(:) = a(:) - 2.0_dp*scapro*d(:)
1110
1111 ELSE
1112
1113 a_mirror(:) = 0.0_dp
1114
1115 END IF
1116
1117 END FUNCTION reflect_vector
1118
1119! **************************************************************************************************
1120!> \brief Rotation of the vector a about an rotation axis defined by the
1121!> vector b. The rotation angle is phi (radians). The rotated vector
1122!> a is stored in a_rot.
1123!> \param a ...
1124!> \param phi ...
1125!> \param b ...
1126!> \return ...
1127!> \date 16.10.1998
1128!> \author MK
1129!> \version 1.0
1130! **************************************************************************************************
1131 PURE FUNCTION rotate_vector(a, phi, b) RESULT(a_rot)
1132 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a
1133 REAL(kind=dp), INTENT(IN) :: phi
1134 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: b
1135 REAL(kind=dp), DIMENSION(3) :: a_rot
1136
1137 REAL(kind=dp) :: length_of_b
1138 REAL(kind=dp), DIMENSION(3, 3) :: rotmat
1139
1140 length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1141 IF (length_of_b > eps_geo) THEN
1142
1143 ! Build up the rotation matrix rotmat
1144 CALL build_rotmat(phi, b, rotmat)
1145
1146 ! Rotate the vector a by phi about the axis defined by vector b
1147 a_rot(:) = matmul(rotmat, a)
1148
1149 ELSE
1150
1151 a_rot(:) = 0.0_dp
1152
1153 END IF
1154
1155 END FUNCTION rotate_vector
1156
1157! **************************************************************************************************
1158!> \brief Set the diagonal elements of matrix a to b.
1159!> \param a ...
1160!> \param b ...
1161!> \date 20.11.1998
1162!> \author MK
1163!> \version 1.0
1164! **************************************************************************************************
1165 PURE SUBROUTINE set_diag_scalar_d(a, b)
1166 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1167 REAL(kind=dp), INTENT(IN) :: b
1168
1169 INTEGER :: i, n
1170
1171 n = min(SIZE(a, 1), SIZE(a, 2))
1172 DO i = 1, n
1173 a(i, i) = b
1174 END DO
1175
1176 END SUBROUTINE set_diag_scalar_d
1177
1178! **************************************************************************************************
1179!> \brief ...
1180!> \param a ...
1181!> \param b ...
1182! **************************************************************************************************
1183 PURE SUBROUTINE set_diag_scalar_z(a, b)
1184 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1185 COMPLEX(KIND=dp), INTENT(IN) :: b
1186
1187 INTEGER :: i, n
1188
1189 n = min(SIZE(a, 1), SIZE(a, 2))
1190 DO i = 1, n
1191 a(i, i) = b
1192 END DO
1193
1194 END SUBROUTINE set_diag_scalar_z
1195
1196! **************************************************************************************************
1197!> \brief Symmetrize the matrix a.
1198!> \param a ...
1199!> \param option ...
1200!> \date 16.10.1998
1201!> \author MK
1202!> \version 1.0
1203! **************************************************************************************************
1204 SUBROUTINE symmetrize_matrix(a, option)
1205 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1206 CHARACTER(LEN=*), INTENT(IN) :: option
1207
1208 INTEGER :: i, n
1209
1210 n = min(SIZE(a, 1), SIZE(a, 2))
1211
1212 IF (option == "lower_to_upper") THEN
1213 DO i = 1, n - 1
1214 a(i, i + 1:n) = a(i + 1:n, i)
1215 END DO
1216 ELSE IF (option == "upper_to_lower") THEN
1217 DO i = 1, n - 1
1218 a(i + 1:n, i) = a(i, i + 1:n)
1219 END DO
1220 ELSE IF (option == "anti_lower_to_upper") THEN
1221 DO i = 1, n - 1
1222 a(i, i + 1:n) = -a(i + 1:n, i)
1223 END DO
1224 ELSE IF (option == "anti_upper_to_lower") THEN
1225 DO i = 1, n - 1
1226 a(i + 1:n, i) = -a(i, i + 1:n)
1227 END DO
1228 ELSE
1229 cpabort("Invalid option <"//trim(option)//"> was specified for parameter #2")
1230 END IF
1231
1232 END SUBROUTINE symmetrize_matrix
1233
1234! **************************************************************************************************
1235!> \brief Set the matrix a to be a unit matrix.
1236!> \param a ...
1237!> \date 16.10.1998
1238!> \author MK
1239!> \version 1.0
1240! **************************************************************************************************
1241 PURE SUBROUTINE unit_matrix_d(a)
1242 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1243
1244 a(:, :) = 0.0_dp
1245 CALL set_diag(a, 1.0_dp)
1246
1247 END SUBROUTINE unit_matrix_d
1248
1249! **************************************************************************************************
1250!> \brief ...
1251!> \param a ...
1252! **************************************************************************************************
1253 PURE SUBROUTINE unit_matrix_z(a)
1254 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1255
1256 a(:, :) = (0.0_dp, 0.0_dp)
1257 CALL set_diag(a, (1.0_dp, 0.0_dp))
1258
1259 END SUBROUTINE unit_matrix_z
1260
1261! **************************************************************************************************
1262!> \brief Calculation of the vector product c = a x b.
1263!> \param a ...
1264!> \param b ...
1265!> \return ...
1266!> \date 16.10.1998
1267!> \author MK
1268!> \version 1.0
1269! **************************************************************************************************
1270 PURE FUNCTION vector_product(a, b) RESULT(c)
1271 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: a, b
1272 REAL(kind=dp), DIMENSION(3) :: c
1273
1274 c(1) = a(2)*b(3) - a(3)*b(2)
1275 c(2) = a(3)*b(1) - a(1)*b(3)
1276 c(3) = a(1)*b(2) - a(2)*b(1)
1277
1278 END FUNCTION vector_product
1279
1280! **************************************************************************************************
1281!> \brief computes the greatest common divisor of two number
1282!> \param a ...
1283!> \param b ...
1284!> \return ...
1285!> \author Joost VandeVondele
1286! **************************************************************************************************
1287 ELEMENTAL FUNCTION gcd(a, b)
1288 INTEGER, INTENT(IN) :: a, b
1289 INTEGER :: gcd
1290
1291 INTEGER :: aa, ab, l, rem, s
1292
1293 aa = abs(a)
1294 ab = abs(b)
1295 IF (aa < ab) THEN
1296 s = aa
1297 l = ab
1298 ELSE
1299 s = ab
1300 l = aa
1301 END IF
1302 IF (s /= 0) THEN
1303 DO
1304 rem = mod(l, s)
1305 IF (rem == 0) EXIT
1306 l = s
1307 s = rem
1308 END DO
1309 gcd = s
1310 ELSE
1311 gcd = l
1312 END IF
1313 END FUNCTION gcd
1314
1315! **************************************************************************************************
1316!> \brief computes the least common multiplier of two numbers
1317!> \param a ...
1318!> \param b ...
1319!> \return ...
1320!> \author Joost VandeVondele
1321! **************************************************************************************************
1322 ELEMENTAL FUNCTION lcm(a, b)
1323 INTEGER, INTENT(IN) :: a, b
1324 INTEGER :: lcm
1325
1326 INTEGER :: tmp
1327
1328 tmp = gcd(a, b)
1329 IF (tmp == 0) THEN
1330 lcm = 0
1331 ELSE
1332 ! could still overflow if the true lcm is larger than maxint
1333 lcm = abs((a/tmp)*b)
1334 END IF
1335 END FUNCTION lcm
1336
1337! **************************************************************************************************
1338!> \brief computes the exponential integral
1339!> Ei(x) = Int(exp(-x*t)/t,t=1..infinity) x>0
1340!> \param x ...
1341!> \return ...
1342!> \author JGH (adapted from Numerical recipies)
1343! **************************************************************************************************
1344 FUNCTION ei(x)
1345 REAL(dp) :: x, ei
1346
1347 INTEGER, PARAMETER :: maxit = 100
1348 REAL(dp), PARAMETER :: eps = epsilon(0.0_dp), &
1349 fpmin = tiny(0.0_dp)
1350
1351 INTEGER :: k
1352 REAL(dp) :: fact, prev, sum1, term
1353
1354 IF (x <= 0._dp) THEN
1355 cpabort("Invalid argument")
1356 END IF
1357
1358 IF (x < fpmin) THEN
1359 ei = log(x) + euler
1360 ELSE IF (x <= -log(eps)) THEN
1361 sum1 = 0._dp
1362 fact = 1._dp
1363 DO k = 1, maxit
1364 fact = fact*x/real(k, dp)
1365 term = fact/real(k, dp)
1366 sum1 = sum1 + term
1367 IF (term < eps*sum1) EXIT
1368 END DO
1369 ei = sum1 + log(x) + euler
1370 ELSE
1371 sum1 = 0._dp
1372 term = 1._dp
1373 DO k = 1, maxit
1374 prev = term
1375 term = term*real(k, dp)/x
1376 IF (term < eps) EXIT
1377 IF (term < prev) THEN
1378 sum1 = sum1 + term
1379 ELSE
1380 sum1 = sum1 - prev
1381 EXIT
1382 END IF
1383 END DO
1384 ei = exp(x)*(1._dp + sum1)/x
1385 END IF
1386
1387 END FUNCTION ei
1388
1389! **************************************************************************************************
1390!> \brief computes the exponential integral
1391!> En(x) = Int(exp(-x*t)/t^n,t=1..infinity) x>0, n=0,1,..
1392!> Note: Ei(-x) = -E1(x)
1393!> \param n ...
1394!> \param x ...
1395!> \return ...
1396!> \par History
1397!> 05.2007 Created
1398!> \author Manuel Guidon (adapted from Numerical recipies)
1399! **************************************************************************************************
1400 ELEMENTAL IMPURE FUNCTION expint(n, x)
1401 INTEGER, INTENT(IN) :: n
1402 REAL(dp), INTENT(IN) :: x
1403 REAL(dp) :: expint
1404
1405 INTEGER, PARAMETER :: maxit = 100
1406 REAL(dp), PARAMETER :: eps = 6.e-14_dp, euler = 0.5772156649015328606065120_dp, &
1407 fpmin = tiny(0.0_dp)
1408
1409 INTEGER :: i, ii, nm1
1410 REAL(dp) :: a, b, c, d, del, fact, h, psi
1411
1412 nm1 = n - 1
1413
1414 IF (n < 0 .OR. x < 0.0_dp .OR. (x == 0.0_dp .AND. (n == 0 .OR. n == 1))) THEN
1415 cpabort("Invalid argument")
1416 ELSE IF (n == 0) THEN !Special case.
1417 expint = exp(-x)/x
1418 ELSE IF (x == 0.0_dp) THEN !Another special case.
1419 expint = 1.0_dp/nm1
1420 ELSE IF (x > 1.0_dp) THEN !Lentz's algorithm (5.2).
1421 b = x + n
1422 c = 1.0_dp/fpmin
1423 d = 1.0_dp/b
1424 h = d
1425 DO i = 1, maxit
1426 a = -i*(nm1 + i)
1427 b = b + 2.0_dp
1428 d = 1.0_dp/(a*d + b)
1429 c = b + a/c
1430 del = c*d
1431 h = h*del
1432 IF (abs(del - 1.0_dp) < eps) THEN
1433 expint = h*exp(-x)
1434 RETURN
1435 END IF
1436 END DO
1437 cpabort("continued fraction failed in expint")
1438 ELSE !Evaluate series.
1439 IF (nm1 /= 0) THEN !Set first term.
1440 expint = 1.0_dp/nm1
1441 ELSE
1442 expint = -log(x) - euler
1443 END IF
1444 fact = 1.0_dp
1445 DO i = 1, maxit
1446 fact = -fact*x/i
1447 IF (i /= nm1) THEN
1448 del = -fact/(i - nm1)
1449 ELSE
1450 psi = -euler !Compute I(n).
1451 DO ii = 1, nm1
1452 psi = psi + 1.0_dp/ii
1453 END DO
1454 del = fact*(-log(x) + psi)
1455 END IF
1456 expint = expint + del
1457 IF (abs(del) < abs(expint)*eps) RETURN
1458 END DO
1459 cpabort("series failed in expint")
1460 END IF
1461
1462 END FUNCTION expint
1463
1464! **************************************************************************************************
1465!> \brief Jacobi matrix diagonalization. The eigenvalues are returned in
1466!> vector d and the eigenvectors are returned in matrix v in ascending
1467!> order.
1468!>
1469!> \param a ...
1470!> \param d ...
1471!> \param v ...
1472!> \par History
1473!> - Creation (20.11.98, Matthias Krack)
1474! **************************************************************************************************
1475 SUBROUTINE jacobi(a, d, v)
1476 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1477 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: d
1478 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: v
1479
1480 INTEGER :: n
1481
1482 n = SIZE(d(:))
1483
1484 ! Diagonalize matrix a
1485 CALL diag(n, a, d, v)
1486
1487 ! Sort eigenvalues and eigenvector in ascending order
1488 CALL eigsrt(n, d, v)
1489
1490 END SUBROUTINE jacobi
1491
1492! **************************************************************************************************
1493!> \brief Diagonalize matrix a. The eigenvalues are returned in vector d
1494!> and the eigenvectors are returned in matrix v.
1495!>
1496!> \param n matrix/vector extent (problem size)
1497!> \param a matrix to be diagonalised
1498!> \param d vector of eigenvalues
1499!> \param v matrix of eigenvectors
1500!> \par History
1501!> - Creation (20.11.98, Matthias Krack)
1502! **************************************************************************************************
1503 SUBROUTINE diag(n, a, d, v)
1504 INTEGER, INTENT(IN) :: n
1505 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1506 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: d
1507 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: v
1508
1509 CHARACTER(len=*), PARAMETER :: routinen = 'diag'
1510 REAL(kind=dp), PARAMETER :: a_eps = 1.0e-10_dp, d_eps = 1.0e-3_dp
1511
1512 INTEGER :: handle, i, ip, iq
1513 REAL(kind=dp) :: a_max, apq, c, d_min, dip, diq, g, h, s, &
1514 t, tau, theta, tresh
1515 REAL(kind=dp), DIMENSION(n) :: b, z
1516
1517 CALL timeset(routinen, handle)
1518
1519 a_max = 0.0_dp
1520 DO ip = 1, n - 1
1521 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1522 b(ip) = a(ip, ip) ! get_diag(a)
1523 END DO
1524 b(n) = a(n, n)
1525
1526 CALL unit_matrix(v)
1527
1528 ! Go for 50 iterations
1529 DO i = 1, 50
1530 d = b
1531 d_min = max(d_eps, minval(abs(b)))
1532 IF (a_max < a_eps*d_min) THEN
1533 CALL timestop(handle)
1534 RETURN
1535 END IF
1536 tresh = merge(a_max, 0.0_dp, (i < 4))
1537 z = 0.0_dp
1538 DO ip = 1, n - 1
1539 DO iq = ip + 1, n
1540 dip = d(ip)
1541 diq = d(iq)
1542 apq = a(ip, iq)
1543 g = 100.0_dp*abs(apq)
1544 IF (tresh < abs(apq)) THEN
1545 h = diq - dip
1546 IF ((abs(h) + g) /= abs(h)) THEN
1547 theta = 0.5_dp*h/apq
1548 t = 1.0_dp/(abs(theta) + sqrt(1.0_dp + theta**2))
1549 IF (theta < 0.0_dp) t = -t
1550 ELSE
1551 t = apq/h
1552 END IF
1553 c = 1.0_dp/sqrt(1.0_dp + t**2)
1554 s = t*c
1555 tau = s/(1.0_dp + c)
1556 h = t*apq
1557 z(ip) = z(ip) - h
1558 z(iq) = z(iq) + h
1559 d(ip) = dip - h
1560 d(iq) = diq + h
1561 a(ip, iq) = 0.0_dp
1562 CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1563 CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1564 CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1565 CALL jrotate(v(:, ip), v(:, iq), s, tau)
1566 ELSE IF ((4 < i) .AND. &
1567 ((abs(dip) + g) == abs(dip)) .AND. &
1568 ((abs(diq) + g) == abs(diq))) THEN
1569 a(ip, iq) = 0.0_dp
1570 END IF
1571 END DO
1572 END DO
1573 b = b + z
1574 a_max = 0.0_dp
1575 DO ip = 1, n - 1
1576 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1577 END DO
1578 END DO
1579 WRITE (*, '(/,T2,A,/)') 'Too many iterations in jacobi'
1580
1581 CALL timestop(handle)
1582
1583 END SUBROUTINE diag
1584
1585! **************************************************************************************************
1586!> \brief Perform a Jacobi rotation of the vectors a and b.
1587!>
1588!> \param a ...
1589!> \param b ...
1590!> \param ss ...
1591!> \param tt ...
1592!> \par History
1593!> - Creation (20.11.98, Matthias Krack)
1594! **************************************************************************************************
1595 PURE SUBROUTINE jrotate(a, b, ss, tt)
1596 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1597 REAL(kind=dp), INTENT(IN) :: ss, tt
1598
1599 REAL(kind=dp) :: u, v
1600
1601 u = 1.0_dp - ss*tt
1602 v = ss/u
1603
1604 a = a*u - b*ss
1605 b = b*(u + ss*v) + a*v
1606
1607 END SUBROUTINE jrotate
1608
1609! **************************************************************************************************
1610!> \brief Sort the values in vector d in ascending order and swap the
1611!> corresponding columns of matrix v.
1612!>
1613!> \param n ...
1614!> \param d ...
1615!> \param v ...
1616!> \par History
1617!> - Creation (20.11.98, Matthias Krack)
1618! **************************************************************************************************
1619 SUBROUTINE eigsrt(n, d, v)
1620 INTEGER, INTENT(IN) :: n
1621 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: d
1622 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: v
1623
1624 INTEGER :: i, j
1625
1626 DO i = 1, n - 1
1627 j = sum(minloc(d(i:n))) + i - 1
1628 IF (j /= i) THEN
1629 CALL swap(d(i), d(j))
1630 CALL swap(v(:, i), v(:, j))
1631 END IF
1632 END DO
1633
1634 END SUBROUTINE eigsrt
1635
1636! **************************************************************************
1637!> \brief Swap two scalars
1638!>
1639!> \param a ...
1640!> \param b ...
1641!> \par History
1642!> - Creation (20.11.98, Matthias Krack)
1643! **************************************************************************************************
1644 ELEMENTAL SUBROUTINE swap_scalar(a, b)
1645 REAL(kind=dp), INTENT(INOUT) :: a, b
1646
1647 REAL(kind=dp) :: c
1648
1649 c = a
1650 a = b
1651 b = c
1652
1653 END SUBROUTINE swap_scalar
1654
1655! **************************************************************************
1656!> \brief Swap two vectors
1657!>
1658!> \param a ...
1659!> \param b ...
1660!> \par History
1661!> - Creation (20.11.98, Matthias Krack)
1662! **************************************************************************************************
1663 SUBROUTINE swap_vector(a, b)
1664 REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: a, b
1665
1666 INTEGER :: i, n
1667 REAL(kind=dp) :: c
1668
1669 n = SIZE(a)
1670
1671 IF (n /= SIZE(b)) THEN
1672 cpabort("Check the array bounds of the parameters")
1673 END IF
1674
1675 DO i = 1, n
1676 c = a(i)
1677 a(i) = b(i)
1678 b(i) = c
1679 END DO
1680
1681 END SUBROUTINE swap_vector
1682
1683! **************************************************************************************************
1684!> \brief - compute a truncation radius for the shortrange operator
1685!> \param eps target accuracy!> \param omg screening parameter
1686!> \param omg ...
1687!> \param r_cutoff cutoff radius
1688!> \par History
1689!> 10.2012 created [Hossein Banihashemian]
1690!> 05.2019 moved here from hfx_types (A. Bussy)
1691!> \author Hossein Banihashemian
1692! **************************************************************************************************
1693 SUBROUTINE erfc_cutoff(eps, omg, r_cutoff)
1694 IMPLICIT NONE
1695
1696 REAL(dp), INTENT(in) :: eps, omg
1697 REAL(dp), INTENT(out) :: r_cutoff
1698
1699 CHARACTER(LEN=*), PARAMETER :: routinen = 'erfc_cutoff'
1700
1701 REAL(dp), PARAMETER :: abstol = 1e-10_dp, soltol = 1e-16_dp
1702 REAL(dp) :: r0, f0, fprime0, delta_r
1703 INTEGER :: iter, handle
1704 INTEGER, PARAMETER :: itermax = 1000
1705
1706 CALL timeset(routinen, handle)
1707
1708 ! initial guess assuming that we are in the asymptotic regime of the erf, and the solution is about 10.
1709 r0 = sqrt(-log(eps*omg*10**2))/omg
1710 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1711
1712 DO iter = 1, itermax
1713 delta_r = f0/fprime0
1714 r0 = r0 - delta_r
1715 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1716 IF (abs(delta_r) < abstol .OR. abs(f0) < soltol) EXIT
1717 END DO
1718 cpassert(iter <= itermax)
1719 r_cutoff = r0
1720
1721 CALL timestop(handle)
1722 CONTAINS
1723! **************************************************************************************************
1724!> \brief ...
1725!> \param r ...
1726!> \param eps ...
1727!> \param omega ...
1728!> \param fn ...
1729!> \param df ...
1730! **************************************************************************************************
1731 ELEMENTAL SUBROUTINE eval_transc_func(r, eps, omega, fn, df)
1732 REAL(dp), INTENT(in) :: r, eps, omega
1733 REAL(dp), INTENT(out) :: fn, df
1734
1735 REAL(dp) :: qr
1736
1737 qr = omega*r
1738 fn = erfc(qr) - r*eps
1739 df = -2.0_dp*oorootpi*omega*exp(-qr**2) - eps
1740 END SUBROUTINE eval_transc_func
1741 END SUBROUTINE erfc_cutoff
1742
1743! **************************************************************************************************
1744!> \brief Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd
1745!> \param matrix Hermitian matrix is preserved
1746!> \param eigenvectors ...
1747!> \param eigenvalues ...
1748!> \author A. Bussy
1749! **************************************************************************************************
1750 SUBROUTINE diag_complex(matrix, eigenvectors, eigenvalues)
1751 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(IN) :: matrix
1752 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
1753 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
1754
1755 CHARACTER(len=*), PARAMETER :: routinen = 'diag_complex'
1756
1757 COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1758 INTEGER :: handle, info, liwork, lrwork, lwork, n
1759 INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1760 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: rwork
1761
1762 CALL timeset(routinen, handle)
1763
1764 IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) cpabort("Expected square matrix")
1765 ! IF (MAXVAL(ABS(matrix - CONJG(TRANSPOSE(matrix)))) > 1e-14_dp) CPABORT("Expected hermitian matrix")
1766
1767 n = SIZE(matrix, 1)
1768 ALLOCATE (iwork(1), rwork(1), work(1))
1769
1770 ! work space query
1771 lwork = -1
1772 lrwork = -1
1773 liwork = -1
1774
1775 CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1776
1777 lwork = ceiling(real(work(1), kind=dp))
1778 lrwork = ceiling(rwork(1))
1779 liwork = iwork(1)
1780
1781 DEALLOCATE (iwork, rwork, work)
1782 ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
1783 eigenvectors(:, :) = matrix(:, :)
1784
1785 ! final diagonalization
1786 CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1787
1788 DEALLOCATE (iwork, rwork, work)
1789
1790 IF (info /= 0) cpabort("Diagonalisation of a complex matrix failed")
1791
1792 CALL timestop(handle)
1793
1794 END SUBROUTINE diag_complex
1795
1796! **************************************************************************************************
1797!> \brief Helper routine for diagonalizing anti symmetric matrices
1798!> \param matrix ...
1799!> \param evecs ...
1800!> \param evals ...
1801! **************************************************************************************************
1802 SUBROUTINE diag_antisym(matrix, evecs, evals)
1803 REAL(dp), DIMENSION(:, :) :: matrix
1804 COMPLEX(dp), DIMENSION(:, :) :: evecs
1805 COMPLEX(dp), DIMENSION(:) :: evals
1806
1807 COMPLEX(dp), ALLOCATABLE, DIMENSION(:, :) :: matrix_c
1808 INTEGER :: n
1809 REAL(dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues
1810
1811 IF (SIZE(matrix, 1) /= SIZE(matrix, 2)) cpabort("Expected square matrix")
1812 ! IF (MAXVAL(ABS(matrix + TRANSPOSE(matrix))) > 1e-14_dp) CPABORT("Expected anti-symmetric matrix")
1813
1814 n = SIZE(matrix, 1)
1815 ALLOCATE (matrix_c(n, n), eigenvalues(n))
1816
1817 matrix_c(:, :) = cmplx(0.0_dp, -matrix, kind=dp)
1818 CALL diag_complex(matrix_c, evecs, eigenvalues)
1819 evals = cmplx(0.0_dp, eigenvalues, kind=dp)
1820
1821 DEALLOCATE (matrix_c, eigenvalues)
1822 END SUBROUTINE diag_antisym
1823! **************************************************************************************************
1824!> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged
1825!> \param A_in Input matrix 1
1826!> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1827!> \param B_in Input matrix 2
1828!> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1829!> \param C_out Output matrix
1830!> \par History
1831!> 11.2025 created [Stepan Marek]
1832! **************************************************************************************************
1833 SUBROUTINE zgemm_square_2(A_in, A_trans, B_in, B_trans, C_out)
1834 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1835 CHARACTER, INTENT(IN) :: A_trans
1836 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1837 CHARACTER, INTENT(IN) :: B_trans
1838 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: C_out
1839
1840 CHARACTER(len=*), PARAMETER :: routineN = 'zgemm_square_2'
1841
1842 INTEGER :: handle, n
1843
1844 n = SIZE(a_in, 1)
1845 IF (n /= SIZE(a_in, 2)) cpabort("Non-square array 1 (A).")
1846 IF (n /= SIZE(b_in, 1)) cpabort("Incompatible (rows) array 2 (B).")
1847 IF (n /= SIZE(b_in, 2)) cpabort("Non-square array 2 (B).")
1848 IF (n /= SIZE(c_out, 1)) cpabort("Incompatible (rows) result array 3 (C).")
1849 IF (n /= SIZE(c_out, 2)) cpabort("Incompatible (cols) result array 3 (C).")
1850 IF (.NOT. (a_trans == 'N' .OR. a_trans == 'n' .OR. &
1851 a_trans == 'T' .OR. a_trans == 't' .OR. &
1852 a_trans == 'C' .OR. a_trans == 'c')) &
1853 cpabort("Unknown transpose character for array 1 (A).")
1854 IF (.NOT. (b_trans == 'N' .OR. b_trans == 'n' .OR. &
1855 b_trans == 'T' .OR. b_trans == 't' .OR. &
1856 b_trans == 'C' .OR. b_trans == 'c')) &
1857 cpabort("Unknown transpose character for array 2 (B).")
1858
1859 CALL timeset(routinen, handle)
1860
1861 CALL zgemm(a_trans, b_trans, n, n, n, z_one, a_in, n, b_in, n, z_zero, c_out, n)
1862
1863 CALL timestop(handle)
1864
1865 END SUBROUTINE zgemm_square_2
1866! **************************************************************************************************
1867!> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, real matrices
1868!> \param A_in Input matrix 1
1869!> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1870!> \param B_in Input matrix 2
1871!> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1872!> \param C_out Output matrix
1873!> \par History
1874!> 11.2025 created [Stepan Marek]
1875! **************************************************************************************************
1876 SUBROUTINE dgemm_square_2(A_in, A_trans, B_in, B_trans, C_out)
1877 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a_in
1878 CHARACTER, INTENT(IN) :: A_trans
1879 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: b_in
1880 CHARACTER, INTENT(IN) :: B_trans
1881 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: c_out
1882
1883 CHARACTER(len=*), PARAMETER :: routineN = 'dgemm_square_2'
1884
1885 INTEGER :: handle, n
1886
1887 n = SIZE(a_in, 1)
1888 IF (n /= SIZE(a_in, 2)) cpabort("Non-square array 1 (A).")
1889 IF (n /= SIZE(b_in, 1)) cpabort("Incompatible (rows) array 2 (B).")
1890 IF (n /= SIZE(b_in, 2)) cpabort("Non-square array 2 (B).")
1891 IF (n /= SIZE(c_out, 1)) cpabort("Incompatible (rows) result array 3 (C).")
1892 IF (n /= SIZE(c_out, 2)) cpabort("Incompatible (cols) result array 3 (C).")
1893 IF (.NOT. (a_trans == 'N' .OR. a_trans == 'n' .OR. &
1894 a_trans == 'T' .OR. a_trans == 't' .OR. &
1895 a_trans == 'C' .OR. a_trans == 'c')) &
1896 cpabort("Unknown transpose character for array 1 (A).")
1897 IF (.NOT. (b_trans == 'N' .OR. b_trans == 'n' .OR. &
1898 b_trans == 'T' .OR. b_trans == 't' .OR. &
1899 b_trans == 'C' .OR. b_trans == 'c')) &
1900 cpabort("Unknown transpose character for array 2 (B).")
1901
1902 CALL timeset(routinen, handle)
1903
1904 CALL dgemm(a_trans, b_trans, n, n, n, 1.0_dp, a_in, n, b_in, n, 0.0_dp, c_out, n)
1905
1906 CALL timestop(handle)
1907
1908 END SUBROUTINE dgemm_square_2
1909! **************************************************************************************************
1910!> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, for 3 matrices
1911!> \param A_in Input matrix 1
1912!> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1913!> \param B_in Input matrix 2
1914!> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1915!> \param C_in Input matrix 3
1916!> \param C_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 3
1917!> \param D_out Output matrix
1918!> \par History
1919!> 11.2025 created [Stepan Marek]
1920! **************************************************************************************************
1921 SUBROUTINE zgemm_square_3(A_in, A_trans, B_in, B_trans, C_in, C_trans, D_out)
1922 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: A_in
1923 CHARACTER, INTENT(IN) :: A_trans
1924 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: B_in
1925 CHARACTER, INTENT(IN) :: B_trans
1926 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: C_in
1927 CHARACTER, INTENT(IN) :: C_trans
1928 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: D_out
1929
1930 CHARACTER(len=*), PARAMETER :: routineN = 'zgemm_square_3'
1931
1932 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1933 INTEGER :: handle, n
1934
1935 n = SIZE(a_in, 1)
1936 IF (n /= SIZE(a_in, 2)) cpabort("Non-square array 1 (A).")
1937 IF (n /= SIZE(b_in, 1)) cpabort("Incompatible (rows) array 2 (B).")
1938 IF (n /= SIZE(b_in, 2)) cpabort("Non-square array 2 (B).")
1939 IF (n /= SIZE(c_in, 1)) cpabort("Incompatible (rows) array 3 (C).")
1940 IF (n /= SIZE(c_in, 2)) cpabort("Non-square array 3 (C).")
1941 IF (n /= SIZE(d_out, 1)) cpabort("Incompatible (rows) result array 4 (D).")
1942 IF (n /= SIZE(d_out, 2)) cpabort("Incompatible (cols) result array 4 (D).")
1943 IF (.NOT. (a_trans == 'N' .OR. a_trans == 'n' .OR. &
1944 a_trans == 'T' .OR. a_trans == 't' .OR. &
1945 a_trans == 'C' .OR. a_trans == 'c')) &
1946 cpabort("Unknown transpose character for array 1 (A).")
1947 IF (.NOT. (b_trans == 'N' .OR. b_trans == 'n' .OR. &
1948 b_trans == 'T' .OR. b_trans == 't' .OR. &
1949 b_trans == 'C' .OR. b_trans == 'c')) &
1950 cpabort("Unknown transpose character for array 2 (B).")
1951 IF (.NOT. (c_trans == 'N' .OR. c_trans == 'n' .OR. &
1952 c_trans == 'T' .OR. c_trans == 't' .OR. &
1953 c_trans == 'C' .OR. c_trans == 'c')) &
1954 cpabort("Unknown transpose character for array 3 (C).")
1955
1956 CALL timeset(routinen, handle)
1957
1958 ALLOCATE (work(n, n), source=z_zero)
1959
1960 CALL zgemm(a_trans, b_trans, n, n, n, z_one, a_in, n, b_in, n, z_zero, work, n)
1961 CALL zgemm('N', c_trans, n, n, n, z_one, work, n, c_in, n, z_zero, d_out, n)
1962
1963 DEALLOCATE (work)
1964
1965 CALL timestop(handle)
1966
1967 END SUBROUTINE zgemm_square_3
1968! **************************************************************************************************
1969!> \brief Square array multiplication via LAPACK routines, leaves inputs unchanged, for 3 matrices
1970!> \param A_in Input matrix 1
1971!> \param A_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 1
1972!> \param B_in Input matrix 2
1973!> \param B_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 2
1974!> \param C_in Input matrix 3
1975!> \param C_trans 'N' - no transpose, 'T' - transpose, 'C' - hermitian conj. of matrix 3
1976!> \param D_out Output matrix
1977!> \par History
1978!> 11.2025 created [Stepan Marek]
1979! **************************************************************************************************
1980 SUBROUTINE dgemm_square_3(A_in, A_trans, B_in, B_trans, C_in, C_trans, D_out)
1981 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: a_in
1982 CHARACTER, INTENT(IN) :: A_trans
1983 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: b_in
1984 CHARACTER, INTENT(IN) :: B_trans
1985 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: c_in
1986 CHARACTER, INTENT(IN) :: C_trans
1987 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: d_out
1988
1989 CHARACTER(len=*), PARAMETER :: routineN = 'dgemm_square_3'
1990
1991 INTEGER :: handle, n
1992 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1993
1994 n = SIZE(a_in, 1)
1995 IF (n /= SIZE(a_in, 2)) cpabort("Non-square array 1 (A).")
1996 IF (n /= SIZE(b_in, 1)) cpabort("Incompatible (rows) array 2 (B).")
1997 IF (n /= SIZE(b_in, 2)) cpabort("Non-square array 2 (B).")
1998 IF (n /= SIZE(c_in, 1)) cpabort("Incompatible (rows) array 3 (C).")
1999 IF (n /= SIZE(c_in, 2)) cpabort("Non-square array 3 (C).")
2000 IF (n /= SIZE(d_out, 1)) cpabort("Incompatible (rows) result array 4 (D).")
2001 IF (n /= SIZE(d_out, 2)) cpabort("Incompatible (cols) result array 4 (D).")
2002 IF (.NOT. (a_trans == 'N' .OR. a_trans == 'n' .OR. &
2003 a_trans == 'T' .OR. a_trans == 't' .OR. &
2004 a_trans == 'C' .OR. a_trans == 'c')) &
2005 cpabort("Unknown transpose character for array 1 (A).")
2006 IF (.NOT. (b_trans == 'N' .OR. b_trans == 'n' .OR. &
2007 b_trans == 'T' .OR. b_trans == 't' .OR. &
2008 b_trans == 'C' .OR. b_trans == 'c')) &
2009 cpabort("Unknown transpose character for array 2 (B).")
2010 IF (.NOT. (c_trans == 'N' .OR. c_trans == 'n' .OR. &
2011 c_trans == 'T' .OR. c_trans == 't' .OR. &
2012 c_trans == 'C' .OR. c_trans == 'c')) &
2013 cpabort("Unknown transpose character for array 3 (C).")
2014
2015 CALL timeset(routinen, handle)
2016
2017 ALLOCATE (work(n, n), source=0.0_dp)
2018
2019 CALL dgemm(a_trans, b_trans, n, n, n, 1.0_dp, a_in, n, b_in, n, 0.0_dp, work, n)
2020 CALL dgemm('N', c_trans, n, n, n, 1.0_dp, work, n, c_in, n, 0.0_dp, d_out, n)
2021
2022 DEALLOCATE (work)
2023
2024 CALL timestop(handle)
2025
2026 END SUBROUTINE dgemm_square_3
2027
2028! **************************************************************************************************
2029!> \brief Solve the generalized eigenvalue equation for complex matrices A*v = B*v*λ
2030!> \param A_in ...
2031!> \param B_in ...
2032!> \param eigenvalues ...
2033!> \param eigenvectors ...
2034!> \author Shridhar Shanbhag
2035! **************************************************************************************************
2036 SUBROUTINE geeig_right(A_in, B_in, eigenvalues, eigenvectors)
2037 COMPLEX(kind=dp), DIMENSION(:, :), INTENT(IN) :: a_in, b_in
2038 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: eigenvalues
2039 COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: eigenvectors
2040
2041 CHARACTER(len=*), PARAMETER :: routinen = 'geeig_right'
2042 COMPLEX(KIND=dp), PARAMETER :: cone = cmplx(1.0_dp, 0.0_dp, kind=dp), &
2043 czero = cmplx(0.0_dp, 0.0_dp, kind=dp)
2044
2045 COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cevals
2046 COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, work
2047 INTEGER :: handle, i, icol, irow, lda, ldb, ldc, &
2048 nao, nc, ncol, nmo, nx
2049 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: evals
2050
2051 CALL timeset(routinen, handle)
2052
2053 ! Test sizes
2054 nao = SIZE(a_in, 1)
2055 nmo = SIZE(eigenvalues)
2056 ALLOCATE (evals(nao), cevals(nao))
2057 ALLOCATE (work(nao, nao), b(nao, nao), a(nao, nao))
2058 a(:, :) = a_in(:, :)
2059 b(:, :) = b_in(:, :)
2060
2061 ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
2062 b = -b
2063 CALL diag_complex(b, work, evals)
2064 evals(:) = -evals(:)
2065 nc = nao
2066 DO i = 1, nao
2067 IF (evals(i) < -1.0_dp) THEN
2068 nc = i - 1
2069 EXIT
2070 END IF
2071 END DO
2072 cpassert(nc /= 0)
2073
2074 IF (nc /= nao) THEN
2075 IF (nc < nmo) THEN
2076 ! Copy NULL space definition to last vectors of eigenvectors (if needed)
2077 ncol = nmo - nc
2078 CALL zcopy(ncol*nao, work(1, nc + 1), 1, eigenvectors(1, nc + 1), 1)
2079 END IF
2080 ! Set NULL space in eigenvector matrix of S to zero
2081 DO icol = nc + 1, nao
2082 DO irow = 1, nao
2083 work(irow, icol) = czero
2084 END DO
2085 END DO
2086 ! Set small eigenvalues to a dummy save value
2087 evals(nc + 1:nao) = 1.0_dp
2088 END IF
2089 ! Calculate U*s**(-1/2)
2090 cevals(:) = cmplx(1.0_dp/sqrt(evals(:)), 0.0_dp, kind=dp)
2091 DO i = 1, min(SIZE(work, 2), SIZE(cevals))
2092 CALL zscal(min(SIZE(work, 2), SIZE(cevals)), cevals(i), work(1, i), 1)
2093 END DO
2094 ! Reduce to get U^(-C) * H * U^(-1)
2095 CALL gemm_square(work, 'C', a, 'N', b)
2096 CALL gemm_square(b, 'N', work, 'N', a)
2097 IF (nc /= nao) THEN
2098 ! set diagonal values to save large value
2099 DO icol = nc + 1, nao
2100 a(icol, icol) = 10000*cone
2101 END DO
2102 END IF
2103 ! Diagonalize
2104 CALL diag_complex(a, b, evals)
2105 eigenvalues(1:nmo) = evals(1:nmo)
2106 nx = min(nc, nmo)
2107 ! Restore vectors C = U^(-1) * C*
2108 lda = SIZE(work, 1)
2109 ldb = SIZE(b, 1)
2110 ldc = SIZE(eigenvectors, 1)
2111 CALL zgemm("N", "N", nao, nx, nc, cone, work, &
2112 lda, b, ldb, czero, eigenvectors, ldc)
2113
2114 DEALLOCATE (evals, cevals, work, b, a)
2115
2116 CALL timestop(handle)
2117 END SUBROUTINE geeig_right
2118
2119END 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:214
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:1476
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:1323
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:1027
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:1094
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:1751
mathlib::diag_antisym
subroutine, public diag_antisym(matrix, evecs, evals)
Helper routine for diagonalizing anti symmetric matrices.
Definition mathlib.F:1803
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::gcd
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1288
mathlib::symmetrize_matrix
subroutine, public symmetrize_matrix(a, option)
Symmetrize the matrix a.
Definition mathlib.F:1205
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:1504
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:2037
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:1132
mathlib::erfc_cutoff
subroutine, public erfc_cutoff(eps, omg, r_cutoff)
compute a truncation radius for the shortrange operator
Definition mathlib.F:1694
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:1401
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:1271
mathlib::multinomial
pure real(kind=dp) function, public multinomial(n, k)
Calculates the multinomial coefficients.
Definition mathlib.F:263