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