Line data Source code
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 : ! **************************************************************************************************
15 : MODULE 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 :
93 : CONTAINS
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 29870 : 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 29870 : CPASSERT(b > a)
114 29870 : IF (x < a .OR. x > b) THEN
115 : ! outside switching intervall
116 27934 : IF (order > 0) THEN
117 : ! derivatives are 0
118 : fx = 0.0_dp
119 : ELSE
120 13967 : IF (x < a) THEN
121 : ! x < a => f(x) = 1
122 : fx = 1.0_dp
123 : ELSE
124 : ! x > b => f(x) = 0
125 13640 : fx = 0.0_dp
126 : END IF
127 : END IF
128 : ELSE
129 : ! renormalized coordinate
130 1936 : u = (x - a)/(b - a)
131 2904 : SELECT CASE (order)
132 : CASE (0)
133 968 : u2 = u*u
134 968 : u3 = u2*u
135 968 : fx = 1._dp - 10._dp*u3 + 15._dp*u2*u2 - 6._dp*u2*u3
136 : CASE (1)
137 968 : u2 = u*u
138 968 : fx = -30._dp*u2 + 60._dp*u*u2 - 30._dp*u2*u2
139 968 : fx = fx/(b - a)
140 : CASE (2)
141 0 : u2 = u*u
142 0 : fx = -60._dp*u + 180._dp*u2 - 120._dp*u*u2
143 0 : fx = fx/(b - a)**2
144 : CASE DEFAULT
145 1936 : CPABORT('order not defined')
146 : END SELECT
147 : END IF
148 :
149 29870 : 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 367551 : LOGICAL FUNCTION abnormal_value(a)
158 : REAL(KIND=dp) :: a
159 :
160 : CHARACTER(LEN=32) :: buffer
161 :
162 367551 : 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 367551 : WRITE (buffer, *) a
168 367551 : IF (INDEX(buffer, "N") /= 0 .OR. INDEX(buffer, "n") /= 0) abnormal_value = .TRUE.
169 :
170 367551 : 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 754822 : 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 754822 : REAL(KIND=dp), DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
188 :
189 3019288 : length_of_a = SQRT(DOT_PRODUCT(a, a))
190 3019288 : length_of_b = SQRT(DOT_PRODUCT(b, b))
191 :
192 754822 : IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo)) THEN
193 3019288 : a_norm(:) = a(:)/length_of_a
194 3019288 : b_norm(:) = b(:)/length_of_b
195 3019288 : 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 754822 : 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 529061863 : ELEMENTAL FUNCTION binomial(n, k) RESULT(n_over_k)
213 : INTEGER, INTENT(IN) :: n, k
214 : REAL(KIND=dp) :: n_over_k
215 :
216 529061863 : IF ((k >= 0) .AND. (k <= n)) THEN
217 525977677 : n_over_k = fac(n)/(fac(n - k)*fac(k))
218 : ELSE
219 : n_over_k = 0.0_dp
220 : END IF
221 :
222 529061863 : 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 171640 : 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 171640 : IF (k >= 0) THEN
244 : z_over_k = 1.0_dp
245 257460 : DO i = 1, k
246 257460 : 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 171640 : 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 8982 : 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 71856 : IF (ALL(k >= 0) .AND. SUM(k) == n) THEN
270 5280 : denom = 1.0_dp
271 21120 : DO i = 1, SIZE(k)
272 21120 : denom = denom*fac(k(i))
273 : END DO
274 5280 : res = fac(n)/denom
275 : ELSE
276 : res = 0.0_dp
277 : END IF
278 :
279 8982 : 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 15759 : 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 15759 : 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 15759 : IF (length_of_a > eps_geo) THEN
303 :
304 63036 : d(:) = a(:)/length_of_a
305 :
306 15759 : cosp = COS(phi)
307 15759 : sinp = SIN(phi)
308 15759 : cost = 1.0_dp - cosp
309 :
310 15759 : rotmat(1, 1) = d(1)*d(1)*cost + cosp
311 15759 : rotmat(1, 2) = d(1)*d(2)*cost - d(3)*sinp
312 15759 : rotmat(1, 3) = d(1)*d(3)*cost + d(2)*sinp
313 15759 : rotmat(2, 1) = d(2)*d(1)*cost + d(3)*sinp
314 15759 : rotmat(2, 2) = d(2)*d(2)*cost + cosp
315 15759 : rotmat(2, 3) = d(2)*d(3)*cost - d(1)*sinp
316 15759 : rotmat(3, 1) = d(3)*d(1)*cost - d(2)*sinp
317 15759 : rotmat(3, 2) = d(3)*d(2)*cost + d(1)*sinp
318 15759 : rotmat(3, 3) = d(3)*d(3)*cost + cosp
319 : ELSE
320 0 : CALL unit_matrix(rotmat)
321 : END IF
322 :
323 15759 : 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 19082358 : 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 19082358 : a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
340 :
341 19082358 : 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 5460 : 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 5460 : a3(1)*(a1(2)*a2(3) - a2(2)*a1(3))
360 :
361 5460 : 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 85030 : 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 85030 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
388 : INTEGER, EXTERNAL :: ilaenv
389 : LOGICAL :: divide_and_conquer
390 85030 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
391 :
392 : EXTERNAL dsyev, dsyevd
393 :
394 85030 : CALL timeset(routineN, handle)
395 :
396 : ! Get the size of the matrix a
397 85030 : n = SIZE(a, 1)
398 :
399 : ! Check the size of matrix a
400 85030 : IF (SIZE(a, 2) /= n) THEN
401 0 : CPABORT("Check the size of matrix a (parameter #1)")
402 : END IF
403 :
404 : ! Check the size of vector eigval
405 85030 : IF (SIZE(eigval) /= n) THEN
406 0 : 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 85030 : IF (PRESENT(dac)) THEN
412 217 : divide_and_conquer = dac
413 : ELSE
414 : divide_and_conquer = .FALSE.
415 : END IF
416 :
417 : ! Get the optimal work storage size
418 :
419 217 : IF (divide_and_conquer) THEN
420 217 : lwork = 2*n**2 + 6*n + 1
421 217 : liwork = 5*n + 3
422 : ELSE
423 84813 : nb = ilaenv(1, "DSYTRD", "U", n, -1, -1, -1)
424 84813 : lwork = (nb + 2)*n
425 : END IF
426 :
427 : ! Allocate work storage
428 :
429 255090 : ALLOCATE (work(lwork))
430 85030 : IF (divide_and_conquer) THEN
431 651 : ALLOCATE (iwork(liwork))
432 : END IF
433 :
434 : ! Diagonalize the matrix a
435 :
436 85030 : info = 0
437 : IF (divide_and_conquer) THEN
438 217 : CALL dsyevd("V", "U", n, a, n, eigval, work, lwork, iwork, liwork, info)
439 : ELSE
440 85009 : CALL dsyev("V", "U", n, a, n, eigval, work, lwork, info)
441 : END IF
442 :
443 85030 : IF (info /= 0) THEN
444 0 : IF (divide_and_conquer) THEN
445 0 : CPABORT("The matrix diagonalization with dsyevd failed")
446 : ELSE
447 0 : CPABORT("The matrix diagonalization with dsyev failed")
448 : END IF
449 : END IF
450 :
451 : ! Release work storage
452 85030 : DEALLOCATE (work)
453 :
454 85030 : IF (divide_and_conquer) THEN
455 217 : DEALLOCATE (iwork)
456 : END IF
457 :
458 85030 : CALL timestop(handle)
459 :
460 170060 : 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 116 : 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 464 : abc = vector_product(bc, -ab)
482 464 : 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 464 : det_abcd = det_3x3(abc, bcd, -bc)
487 116 : dihedral_angle_abcd = SIGN(1.0_dp, det_abcd)*angle(abc, bcd)
488 :
489 116 : 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 63538 : 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 63538 : n = MIN(SIZE(a, 1), SIZE(a, 2))
507 :
508 6369658 : DO i = 1, n
509 6369658 : a_diag(i) = a(i, i)
510 : END DO
511 :
512 63538 : 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 18318003 : 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 18318003 : det_a = 1.0_dp/det_3x3(a)
529 :
530 18318003 : a_inv(1, 1) = (a(2, 2)*a(3, 3) - a(3, 2)*a(2, 3))*det_a
531 18318003 : a_inv(2, 1) = (a(2, 3)*a(3, 1) - a(3, 3)*a(2, 1))*det_a
532 18318003 : a_inv(3, 1) = (a(2, 1)*a(3, 2) - a(3, 1)*a(2, 2))*det_a
533 :
534 18318003 : a_inv(1, 2) = (a(1, 3)*a(3, 2) - a(3, 3)*a(1, 2))*det_a
535 18318003 : a_inv(2, 2) = (a(1, 1)*a(3, 3) - a(3, 1)*a(1, 3))*det_a
536 18318003 : a_inv(3, 2) = (a(1, 2)*a(3, 1) - a(3, 2)*a(1, 1))*det_a
537 :
538 18318003 : a_inv(1, 3) = (a(1, 2)*a(2, 3) - a(2, 2)*a(1, 3))*det_a
539 18318003 : a_inv(2, 3) = (a(1, 3)*a(2, 1) - a(2, 3)*a(1, 1))*det_a
540 18318003 : a_inv(3, 3) = (a(1, 1)*a(2, 2) - a(2, 1)*a(1, 2))*det_a
541 :
542 18318003 : 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 8443 : 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 8443 : INTEGER, ALLOCATABLE :: ipiv(:)
557 8443 : REAL(KIND=dp), ALLOCATABLE :: work(:)
558 :
559 8443 : CALL timeset(routineN, handle)
560 :
561 8443 : n = SIZE(a, 1)
562 8443 : lwork = 20*n
563 25329 : ALLOCATE (ipiv(n))
564 25329 : ALLOCATE (work(lwork))
565 8443 : ipiv = 0
566 8443 : work = 0._dp
567 8443 : info = 0
568 551377 : CALL dgetrf(n, n, a, n, ipiv, info)
569 8443 : IF (info == 0) THEN
570 550793 : CALL dgetri(n, a, n, ipiv, work, lwork, info)
571 : END IF
572 8443 : DEALLOCATE (ipiv, work)
573 :
574 8443 : CALL timestop(handle)
575 :
576 8443 : 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 3026 : 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 3026 : CALL timeset(routineN, handle)
598 :
599 3026 : do_potrf = .TRUE.
600 3026 : IF (PRESENT(potrf)) do_potrf = potrf
601 :
602 3026 : myuplo = 'U'
603 3026 : IF (PRESENT(uplo)) myuplo = uplo
604 :
605 3026 : n = SIZE(a, 1)
606 3026 : info = 0
607 :
608 : ! do cholesky decomposition
609 3026 : IF (do_potrf) THEN
610 4009 : CALL dpotrf(myuplo, n, a, n, info)
611 2809 : IF (info /= 0) CPABORT("DPOTRF failed")
612 : END IF
613 :
614 : ! do inversion using the cholesky decomposition
615 4226 : CALL dpotri(myuplo, n, a, n, info)
616 3026 : IF (info /= 0) CPABORT("Matrix inversion failed")
617 :
618 : ! complete the matrix
619 3026 : IF ((myuplo == "U") .OR. (myuplo == "u")) THEN
620 3026 : CALL symmetrize_matrix(a, "upper_to_lower")
621 : ELSE
622 0 : CALL symmetrize_matrix(a, "lower_to_upper")
623 : END IF
624 :
625 3026 : CALL timestop(handle)
626 :
627 3026 : 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 4676 : 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 4676 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv, iwork
662 : LOGICAL :: do_improve
663 : REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
664 4676 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: berr, ferr, work
665 4676 : 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 4676 : IF (PRESENT(option)) THEN
672 524 : trans = option
673 : ELSE
674 4152 : trans = "N"
675 : END IF
676 :
677 4676 : IF (PRESENT(improve)) THEN
678 648 : do_improve = improve
679 : ELSE
680 : do_improve = .TRUE.
681 : END IF
682 :
683 : ! Get the dimension of matrix a
684 4676 : n = SIZE(a, 1)
685 :
686 : ! Check array dimensions
687 4676 : IF (n == 0) THEN
688 0 : CPABORT("Matrix to be inverted of zero size")
689 : END IF
690 :
691 4676 : IF (n /= SIZE(a, 2)) THEN
692 0 : CPABORT("Check the array bounds of parameter #1")
693 : END IF
694 :
695 4676 : IF ((n /= SIZE(a_inverse, 1)) .OR. &
696 : (n /= SIZE(a_inverse, 2))) THEN
697 0 : CPABORT("Check the array bounds of parameter #2")
698 : END IF
699 :
700 : ! Allocate work storage
701 18704 : ALLOCATE (a_lu(n, n))
702 14028 : ALLOCATE (b(n, n))
703 14028 : ALLOCATE (berr(n))
704 9352 : ALLOCATE (ferr(n))
705 14028 : ALLOCATE (ipiv(n))
706 9352 : ALLOCATE (iwork(n))
707 14028 : ALLOCATE (work(4*n))
708 :
709 2932192 : a_lu(1:n, 1:n) = a(1:n, 1:n)
710 :
711 : ! Compute the LU factorization of the matrix a
712 4676 : CALL dgetrf(n, n, a_lu, n, ipiv, info)
713 :
714 4676 : IF (info /= 0) THEN
715 0 : CPABORT("The LU factorization in dgetrf failed")
716 : END IF
717 :
718 : ! Compute the norm of the matrix a
719 :
720 4676 : IF (trans == "N") THEN
721 4426 : norm = '1'
722 : ELSE
723 250 : norm = 'I'
724 : END IF
725 :
726 4676 : a_norm = dlange(norm, n, n, a, n, work)
727 :
728 : ! Compute the reciprocal of the condition number of a
729 :
730 4676 : CALL dgecon(norm, n, a_lu, n, a_norm, r_cond, work, iwork, info)
731 :
732 4676 : IF (info /= 0) THEN
733 0 : CPABORT("The computation of the condition number in dgecon failed")
734 : END IF
735 :
736 4676 : IF (r_cond < EPSILON(0.0_dp)) THEN
737 0 : 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 0 : "working precision)")
741 : END IF
742 :
743 : ! Solve a system of linear equations using the LU factorization computed by dgetrf
744 :
745 4676 : CALL unit_matrix(a_inverse)
746 :
747 4676 : CALL dgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
748 :
749 4676 : IF (info /= 0) THEN
750 0 : CPABORT("Solving the system of linear equations in dgetrs failed")
751 : END IF
752 :
753 : ! Improve the computed solution iteratively
754 4676 : CALL unit_matrix(b) ! Initialize right-hand sides
755 :
756 4676 : eval_error = 0.0_dp
757 :
758 4676 : IF (do_improve) THEN
759 11948 : 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 11676 : work, iwork, info)
763 :
764 11676 : IF (info /= 0) THEN
765 0 : CPABORT("Improving the computed solution in dgerfs failed")
766 : END IF
767 :
768 11676 : old_eval_error = eval_error
769 230818 : eval_error = MAXVAL(ferr)
770 :
771 11948 : 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 4676 : DEALLOCATE (work)
778 4676 : DEALLOCATE (iwork)
779 4676 : DEALLOCATE (ipiv)
780 4676 : DEALLOCATE (ferr)
781 4676 : DEALLOCATE (berr)
782 4676 : DEALLOCATE (b)
783 4676 : DEALLOCATE (a_lu)
784 :
785 4676 : 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 0 : 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 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work
816 0 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a_lu, b
817 : INTEGER :: info, iter, n
818 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
819 : REAL(KIND=dp) :: a_norm, old_eval_error, r_cond
820 0 : 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 0 : IF (PRESENT(option)) THEN
827 0 : trans = option
828 : ELSE
829 0 : trans = "N"
830 : END IF
831 :
832 : ! Get the dimension of matrix a
833 0 : n = SIZE(a, 1)
834 :
835 : ! Check array dimensions
836 0 : IF (n == 0) THEN
837 0 : CPABORT("Matrix to be inverted of zero size")
838 : END IF
839 :
840 0 : IF (n /= SIZE(a, 2)) THEN
841 0 : CPABORT("Check the array bounds of parameter #1")
842 : END IF
843 :
844 0 : IF ((n /= SIZE(a_inverse, 1)) .OR. &
845 : (n /= SIZE(a_inverse, 2))) THEN
846 0 : CPABORT("Check the array bounds of parameter #2")
847 : END IF
848 :
849 : ! Allocate work storage
850 0 : ALLOCATE (a_lu(n, n))
851 0 : ALLOCATE (b(n, n))
852 0 : ALLOCATE (berr(n))
853 0 : ALLOCATE (ferr(n))
854 0 : ALLOCATE (ipiv(n))
855 0 : ALLOCATE (rwork(2*n))
856 0 : ALLOCATE (work(2*n))
857 :
858 0 : a_lu(1:n, 1:n) = a(1:n, 1:n)
859 :
860 : ! Compute the LU factorization of the matrix a
861 0 : CALL zgetrf(n, n, a_lu, n, ipiv, info)
862 :
863 0 : IF (info /= 0) THEN
864 0 : CPABORT("The LU factorization in dgetrf failed")
865 : END IF
866 :
867 : ! Compute the norm of the matrix a
868 :
869 0 : IF (trans == "N") THEN
870 0 : norm = '1'
871 : ELSE
872 0 : norm = 'I'
873 : END IF
874 :
875 0 : a_norm = zlange(norm, n, n, a, n, work)
876 :
877 : ! Compute the reciprocal of the condition number of a
878 :
879 0 : CALL zgecon(norm, n, a_lu, n, a_norm, r_cond, work, rwork, info)
880 :
881 0 : IF (info /= 0) THEN
882 0 : CPABORT("The computation of the condition number in dgecon failed")
883 : END IF
884 :
885 0 : IF (r_cond < EPSILON(0.0_dp)) THEN
886 0 : 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 0 : "working precision)")
890 : END IF
891 :
892 : ! Solve a system of linear equations using the LU factorization computed by dgetrf
893 :
894 0 : CALL unit_matrix(a_inverse)
895 :
896 0 : CALL zgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
897 :
898 0 : IF (info /= 0) THEN
899 0 : CPABORT("Solving the system of linear equations in dgetrs failed")
900 : END IF
901 :
902 : ! Improve the computed solution iteratively
903 0 : CALL unit_matrix(b) ! Initialize right-hand sides
904 :
905 0 : eval_error = 0.0_dp
906 :
907 0 : 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 0 : work, rwork, info)
911 :
912 0 : IF (info /= 0) THEN
913 0 : CPABORT("Improving the computed solution in dgerfs failed")
914 : END IF
915 :
916 0 : old_eval_error = eval_error
917 0 : eval_error = MAXVAL(ferr)
918 :
919 0 : IF (ABS(eval_error - old_eval_error) <= EPSILON(1.0_dp)) EXIT
920 :
921 : END DO
922 :
923 : ! Release work storage
924 0 : DEALLOCATE (work)
925 0 : DEALLOCATE (rwork)
926 0 : DEALLOCATE (ipiv)
927 0 : DEALLOCATE (ferr)
928 0 : DEALLOCATE (berr)
929 0 : DEALLOCATE (b)
930 0 : DEALLOCATE (a_lu)
931 :
932 0 : 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 861 : 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 861 : INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
955 861 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sig, work
956 861 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sig_plus, temp_mat, u, vt
957 :
958 861 : CALL timeset(routineN, handle)
959 :
960 861 : n = SIZE(a, 1)
961 12054 : ALLOCATE (u(n, n), vt(n, n), sig(n), sig_plus(n, n), iwork(8*n), work(1), temp_mat(n, n))
962 861 : u(:, :) = 0.0_dp
963 861 : vt(:, :) = 0.0_dp
964 861 : sig(:) = 0.0_dp
965 861 : sig_plus = 0.0_dp
966 861 : work = 0.0_dp
967 861 : iwork = 0
968 861 : IF (PRESENT(determinant)) determinant = 1.0_dp
969 :
970 : ! work size query
971 861 : 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 861 : lwork, iwork(1), info)
974 :
975 861 : IF (info /= 0) THEN
976 0 : CPABORT("ERROR in DGESDD: Could not retrieve work array sizes")
977 : END IF
978 861 : lwork = INT(work(1))
979 861 : DEALLOCATE (work)
980 2583 : 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 861 : lwork, iwork(1), info)
985 :
986 861 : IF (info /= 0) THEN
987 0 : CPABORT("SVD failed")
988 : END IF
989 :
990 861 : IF (PRESENT(sval)) THEN
991 0 : CPASSERT(.NOT. ASSOCIATED(sval))
992 0 : ALLOCATE (sval(n))
993 0 : sval(:) = sig
994 : END IF
995 :
996 : ! set singular values that are too small to zero
997 5141 : DO i = 1, n
998 63417 : IF (sig(i) > rskip*MAXVAL(sig)) THEN
999 4280 : IF (PRESENT(determinant)) &
1000 0 : determinant = determinant*sig(i)
1001 4280 : sig_plus(i, i) = 1._dp/sig(i)
1002 : ELSE
1003 0 : sig_plus(i, i) = 0.0_dp
1004 : END IF
1005 : END DO
1006 :
1007 : ! build pseudoinverse: V*sig_plus*UT
1008 861 : CALL dgemm("N", "T", n, n, n, 1._dp, sig_plus, n, u, n, 0._dp, temp_mat, n)
1009 861 : CALL dgemm("T", "N", n, n, n, 1._dp, vt, n, temp_mat, n, 0._dp, a_pinverse, n)
1010 :
1011 861 : DEALLOCATE (u, vt, sig, iwork, work, sig_plus, temp_mat)
1012 :
1013 861 : CALL timestop(handle)
1014 :
1015 861 : 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 1177 : 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 1177 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eig, work
1033 1177 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dinv, temp_mat
1034 :
1035 1177 : CALL timeset(routineN, handle)
1036 :
1037 1177 : info = 0
1038 1177 : n = SIZE(a, 1)
1039 9416 : ALLOCATE (dinv(n, n), eig(n), work(1), temp_mat(n, n))
1040 1177 : dinv(:, :) = 0.0_dp
1041 1177 : eig(:) = 0.0_dp
1042 1177 : work(:) = 0.0_dp
1043 1177 : temp_mat = 0.0_dp
1044 :
1045 : ! work size query
1046 1177 : lwork = -1
1047 1177 : CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1048 1177 : IF (info /= 0) THEN
1049 0 : CPABORT("ERROR in DSYEV: Could not retrieve work array sizes")
1050 : END IF
1051 1177 : lwork = INT(work(1))
1052 1177 : DEALLOCATE (work)
1053 3531 : ALLOCATE (work(lwork))
1054 1177 : work = 0.0_dp
1055 :
1056 : ! get eigenvalues and eigenvectors
1057 1177 : CALL dsyev('V', 'U', n, a, n, eig(1), work(1), lwork, info)
1058 :
1059 1177 : IF (info /= 0) THEN
1060 0 : CPABORT("Matrix diagonalization failed")
1061 : END IF
1062 :
1063 : ! set eigenvalues that are too small to zero
1064 15283 : DO i = 1, n
1065 189031 : IF (eig(i) > rskip*MAXVAL(eig)) THEN
1066 12616 : dinv(i, i) = 1.0_dp/eig(i)
1067 : ELSE
1068 1490 : dinv(i, i) = 0._dp
1069 : END IF
1070 : END DO
1071 :
1072 : ! build pseudoinverse: U*dinv*UT
1073 1177 : CALL dgemm("N", "T", n, n, n, 1._dp, dinv, n, a, n, 0._dp, temp_mat, n)
1074 1177 : CALL dgemm("N", "N", n, n, n, 1._dp, a, n, temp_mat, n, 0._dp, a_pinverse, n)
1075 :
1076 1177 : DEALLOCATE (eig, work, dinv, temp_mat)
1077 :
1078 1177 : CALL timestop(handle)
1079 :
1080 1177 : 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 19487 : 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 19487 : length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1100 :
1101 19487 : IF (length_of_b > eps_geo) THEN
1102 :
1103 77948 : d(:) = b(:)/length_of_b
1104 :
1105 : ! Calculate the mirror image a_mirror of the vector a
1106 19487 : scapro = a(1)*d(1) + a(2)*d(2) + a(3)*d(3)
1107 :
1108 77948 : a_mirror(:) = a(:) - 2.0_dp*scapro*d(:)
1109 :
1110 : ELSE
1111 :
1112 0 : a_mirror(:) = 0.0_dp
1113 :
1114 : END IF
1115 :
1116 19487 : 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 4098 : 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 4098 : length_of_b = SQRT(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1140 4098 : IF (length_of_b > eps_geo) THEN
1141 :
1142 : ! Build up the rotation matrix rotmat
1143 3964 : CALL build_rotmat(phi, b, rotmat)
1144 :
1145 : ! Rotate the vector a by phi about the axis defined by vector b
1146 63424 : a_rot(:) = MATMUL(rotmat, a)
1147 :
1148 : ELSE
1149 :
1150 536 : a_rot(:) = 0.0_dp
1151 :
1152 : END IF
1153 :
1154 4098 : 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 43888 : 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 43888 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1171 244259 : DO i = 1, n
1172 244259 : a(i, i) = b
1173 : END DO
1174 :
1175 43888 : END SUBROUTINE set_diag_scalar_d
1176 :
1177 : ! **************************************************************************************************
1178 : !> \brief ...
1179 : !> \param a ...
1180 : !> \param b ...
1181 : ! **************************************************************************************************
1182 0 : 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 0 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1189 0 : DO i = 1, n
1190 0 : a(i, i) = b
1191 : END DO
1192 :
1193 0 : 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 19696 : 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 19696 : n = MIN(SIZE(a, 1), SIZE(a, 2))
1210 :
1211 19696 : IF (option == "lower_to_upper") THEN
1212 0 : DO i = 1, n - 1
1213 0 : a(i, i + 1:n) = a(i + 1:n, i)
1214 : END DO
1215 19696 : ELSE IF (option == "upper_to_lower") THEN
1216 522697 : DO i = 1, n - 1
1217 29908250 : a(i + 1:n, i) = a(i, i + 1:n)
1218 : END DO
1219 36 : ELSE IF (option == "anti_lower_to_upper") THEN
1220 0 : DO i = 1, n - 1
1221 0 : a(i, i + 1:n) = -a(i + 1:n, i)
1222 : END DO
1223 36 : ELSE IF (option == "anti_upper_to_lower") THEN
1224 564 : DO i = 1, n - 1
1225 4716 : a(i + 1:n, i) = -a(i, i + 1:n)
1226 : END DO
1227 : ELSE
1228 0 : CPABORT("Invalid option <"//TRIM(option)//"> was specified for parameter #2")
1229 : END IF
1230 :
1231 19696 : 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 43888 : PURE SUBROUTINE unit_matrix_d(a)
1241 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1242 :
1243 7727648 : a(:, :) = 0.0_dp
1244 43888 : CALL set_diag(a, 1.0_dp)
1245 :
1246 43888 : END SUBROUTINE unit_matrix_d
1247 :
1248 : ! **************************************************************************************************
1249 : !> \brief ...
1250 : !> \param a ...
1251 : ! **************************************************************************************************
1252 0 : PURE SUBROUTINE unit_matrix_z(a)
1253 : COMPLEX(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: a
1254 :
1255 0 : a(:, :) = (0.0_dp, 0.0_dp)
1256 0 : CALL set_diag(a, (1.0_dp, 0.0_dp))
1257 :
1258 0 : 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 8011 : 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 8011 : c(1) = a(2)*b(3) - a(3)*b(2)
1274 8011 : c(2) = a(3)*b(1) - a(1)*b(3)
1275 8011 : c(3) = a(1)*b(2) - a(2)*b(1)
1276 :
1277 8011 : 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 1070975 : ELEMENTAL FUNCTION gcd(a, b)
1287 : INTEGER, INTENT(IN) :: a, b
1288 : INTEGER :: gcd
1289 :
1290 : INTEGER :: aa, ab, l, rem, s
1291 :
1292 1070975 : aa = ABS(a)
1293 1070975 : ab = ABS(b)
1294 1070975 : IF (aa < ab) THEN
1295 : s = aa
1296 : l = ab
1297 : ELSE
1298 1047227 : s = ab
1299 1047227 : l = aa
1300 : END IF
1301 1070975 : IF (s /= 0) THEN
1302 : DO
1303 1071113 : rem = MOD(l, s)
1304 1071113 : IF (rem == 0) EXIT
1305 : l = s
1306 1070975 : s = rem
1307 : END DO
1308 : GCD = s
1309 : ELSE
1310 : GCD = l
1311 : END IF
1312 1070975 : 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 511454 : ELEMENTAL FUNCTION lcm(a, b)
1322 : INTEGER, INTENT(IN) :: a, b
1323 : INTEGER :: lcm
1324 :
1325 : INTEGER :: tmp
1326 :
1327 511454 : tmp = gcd(a, b)
1328 511454 : IF (tmp == 0) THEN
1329 : lcm = 0
1330 : ELSE
1331 : ! could still overflow if the true lcm is larger than maxint
1332 511454 : lcm = ABS((a/tmp)*b)
1333 : END IF
1334 511454 : 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 0 : 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 0 : IF (x <= 0._dp) THEN
1354 0 : CPABORT("Invalid argument")
1355 : END IF
1356 :
1357 0 : IF (x < fpmin) THEN
1358 0 : ei = LOG(x) + euler
1359 0 : ELSE IF (x <= -LOG(EPS)) THEN
1360 : sum1 = 0._dp
1361 : fact = 1._dp
1362 0 : DO k = 1, maxit
1363 0 : fact = fact*x/REAL(k, dp)
1364 0 : term = fact/REAL(k, dp)
1365 0 : sum1 = sum1 + term
1366 0 : IF (term < eps*sum1) EXIT
1367 : END DO
1368 0 : ei = sum1 + LOG(x) + euler
1369 : ELSE
1370 : sum1 = 0._dp
1371 : term = 1._dp
1372 0 : DO k = 1, maxit
1373 0 : prev = term
1374 0 : term = term*REAL(k, dp)/x
1375 0 : IF (term < eps) EXIT
1376 0 : IF (term < prev) THEN
1377 0 : sum1 = sum1 + term
1378 : ELSE
1379 0 : sum1 = sum1 - prev
1380 0 : EXIT
1381 : END IF
1382 : END DO
1383 0 : ei = EXP(x)*(1._dp + sum1)/x
1384 : END IF
1385 :
1386 0 : 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 298362910 : 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 298362910 : nm1 = n - 1
1412 :
1413 298362910 : IF (n < 0 .OR. x < 0.0_dp .OR. (x == 0.0_dp .AND. (n == 0 .OR. n == 1))) THEN
1414 0 : CPABORT("Invalid argument")
1415 298362910 : ELSE IF (n == 0) THEN !Special case.
1416 0 : expint = EXP(-x)/x
1417 298362910 : ELSE IF (x == 0.0_dp) THEN !Another special case.
1418 0 : expint = 1.0_dp/nm1
1419 298362910 : ELSE IF (x > 1.0_dp) THEN !Lentz's algorithm (5.2).
1420 240863916 : b = x + n
1421 240863916 : c = 1.0_dp/FPMIN
1422 240863916 : d = 1.0_dp/b
1423 240863916 : h = d
1424 6974042655 : DO i = 1, MAXIT
1425 6974042655 : a = -i*(nm1 + i)
1426 6974042655 : b = b + 2.0_dp
1427 6974042655 : d = 1.0_dp/(a*d + b)
1428 6974042655 : c = b + a/c
1429 6974042655 : del = c*d
1430 6974042655 : h = h*del
1431 6974042655 : IF (ABS(del - 1.0_dp) < EPS) THEN
1432 240863916 : expint = h*EXP(-x)
1433 240863916 : RETURN
1434 : END IF
1435 : END DO
1436 0 : CPABORT("continued fraction failed in expint")
1437 : ELSE !Evaluate series.
1438 57498994 : IF (nm1 /= 0) THEN !Set first term.
1439 0 : expint = 1.0_dp/nm1
1440 : ELSE
1441 57498994 : expint = -LOG(x) - euler
1442 : END IF
1443 : fact = 1.0_dp
1444 491523421 : DO i = 1, MAXIT
1445 491523421 : fact = -fact*x/i
1446 491523421 : IF (i /= nm1) THEN
1447 491523421 : del = -fact/(i - nm1)
1448 : ELSE
1449 : psi = -euler !Compute I(n).
1450 0 : DO ii = 1, nm1
1451 0 : psi = psi + 1.0_dp/ii
1452 : END DO
1453 0 : del = fact*(-LOG(x) + psi)
1454 : END IF
1455 491523421 : expint = expint + del
1456 491523421 : IF (ABS(del) < ABS(expint)*EPS) RETURN
1457 : END DO
1458 0 : 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 34478 : 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 34478 : n = SIZE(d(:))
1482 :
1483 : ! Diagonalize matrix a
1484 34478 : CALL diag(n, a, d, v)
1485 :
1486 : ! Sort eigenvalues and eigenvector in ascending order
1487 34478 : CALL eigsrt(n, d, v)
1488 :
1489 34478 : 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 34478 : 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 34478 : REAL(KIND=dp), DIMENSION(n) :: b, z
1515 :
1516 34478 : CALL timeset(routineN, handle)
1517 :
1518 34478 : a_max = 0.0_dp
1519 121765 : DO ip = 1, n - 1
1520 879538 : a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip + 1:n))))
1521 121765 : b(ip) = a(ip, ip) ! get_diag(a)
1522 : END DO
1523 34478 : b(n) = a(n, n)
1524 :
1525 34478 : CALL unit_matrix(v)
1526 :
1527 : ! Go for 50 iterations
1528 160078 : DO i = 1, 50
1529 975669 : d = b
1530 975669 : d_min = MAX(d_eps, MINVAL(ABS(b)))
1531 160078 : IF (a_max < a_eps*d_min) THEN
1532 34478 : CALL timestop(handle)
1533 : RETURN
1534 : END IF
1535 125600 : tresh = MERGE(a_max, 0.0_dp, (i < 4))
1536 819426 : z = 0.0_dp
1537 693826 : DO ip = 1, n - 1
1538 6401612 : DO iq = ip + 1, n
1539 5707786 : dip = d(ip)
1540 5707786 : diq = d(iq)
1541 5707786 : apq = a(ip, iq)
1542 5707786 : g = 100.0_dp*ABS(apq)
1543 6276012 : IF (tresh < ABS(apq)) THEN
1544 3334999 : h = diq - dip
1545 3334999 : IF ((ABS(h) + g) /= ABS(h)) THEN
1546 2523614 : theta = 0.5_dp*h/apq
1547 2523614 : t = 1.0_dp/(ABS(theta) + SQRT(1.0_dp + theta**2))
1548 2523614 : IF (theta < 0.0_dp) t = -t
1549 : ELSE
1550 811385 : t = apq/h
1551 : END IF
1552 3334999 : c = 1.0_dp/SQRT(1.0_dp + t**2)
1553 3334999 : s = t*c
1554 3334999 : tau = s/(1.0_dp + c)
1555 3334999 : h = t*apq
1556 3334999 : z(ip) = z(ip) - h
1557 3334999 : z(iq) = z(iq) + h
1558 3334999 : d(ip) = dip - h
1559 3334999 : d(iq) = diq + h
1560 3334999 : a(ip, iq) = 0.0_dp
1561 3334999 : CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1562 3334999 : CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1563 3334999 : CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1564 3334999 : CALL jrotate(v(:, ip), v(:, iq), s, tau)
1565 : ELSE IF ((4 < i) .AND. &
1566 0 : ((ABS(dip) + g) == ABS(dip)) .AND. &
1567 2372787 : ((ABS(diq) + g) == ABS(diq))) THEN
1568 0 : a(ip, iq) = 0.0_dp
1569 : END IF
1570 : END DO
1571 : END DO
1572 819426 : b = b + z
1573 : a_max = 0.0_dp
1574 693826 : DO ip = 1, n - 1
1575 6401612 : a_max = MAX(a_max, MAXVAL(ABS(a(ip, ip + 1:n))))
1576 : END DO
1577 : END DO
1578 0 : WRITE (*, '(/,T2,A,/)') 'Too many iterations in jacobi'
1579 :
1580 0 : 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 13339996 : 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 13339996 : u = 1.0_dp - ss*tt
1601 13339996 : v = ss/u
1602 :
1603 228402348 : a = a*u - b*ss
1604 228402348 : b = b*(u + ss*v) + a*v
1605 :
1606 13339996 : 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 34478 : 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 121765 : DO i = 1, n - 1
1626 1054112 : j = SUM(MINLOC(d(i:n))) + i - 1
1627 121765 : IF (j /= i) THEN
1628 42828 : CALL swap(d(i), d(j))
1629 42828 : CALL swap(v(:, i), v(:, j))
1630 : END IF
1631 : END DO
1632 :
1633 34478 : 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 42828 : ELEMENTAL SUBROUTINE swap_scalar(a, b)
1644 : REAL(KIND=dp), INTENT(INOUT) :: a, b
1645 :
1646 : REAL(KIND=dp) :: c
1647 :
1648 42828 : c = a
1649 42828 : a = b
1650 42828 : b = c
1651 :
1652 42828 : 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 42828 : 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 42828 : n = SIZE(a)
1669 :
1670 42828 : IF (n /= SIZE(b)) THEN
1671 0 : CPABORT("Check the array bounds of the parameters")
1672 : END IF
1673 :
1674 956900 : DO i = 1, n
1675 914072 : c = a(i)
1676 914072 : a(i) = b(i)
1677 956900 : b(i) = c
1678 : END DO
1679 :
1680 42828 : 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 66 : 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 66 : 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 66 : r0 = SQRT(-LOG(eps*omg*10**2))/omg
1707 66 : CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1708 :
1709 638 : DO iter = 1, iterMAX
1710 638 : delta_r = f0/fprime0
1711 638 : r0 = r0 - delta_r
1712 638 : CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1713 638 : IF (ABS(delta_r) < abstol .OR. ABS(f0) < soltol) EXIT
1714 : END DO
1715 66 : CPASSERT(iter <= itermax)
1716 66 : r_cutoff = r0
1717 :
1718 66 : 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 704 : 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 704 : qr = omega*r
1735 704 : fn = ERFC(qr) - r*eps
1736 704 : df = -2.0_dp*oorootpi*omega*EXP(-qr**2) - eps
1737 704 : 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 68919 : 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 68919 : COMPLEX(KIND=dp), DIMENSION(:), ALLOCATABLE :: work
1755 : INTEGER :: handle, info, liwork, lrwork, lwork, n
1756 68919 : INTEGER, DIMENSION(:), ALLOCATABLE :: iwork
1757 68919 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: rwork
1758 :
1759 68919 : CALL timeset(routineN, handle)
1760 :
1761 68919 : 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 68919 : n = SIZE(matrix, 1)
1765 68919 : ALLOCATE (iwork(1), rwork(1), work(1))
1766 :
1767 : ! work space query
1768 68919 : lwork = -1
1769 68919 : lrwork = -1
1770 68919 : liwork = -1
1771 :
1772 68919 : CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1773 :
1774 68919 : lwork = CEILING(REAL(work(1), KIND=dp))
1775 68919 : lrwork = CEILING(rwork(1))
1776 68919 : liwork = iwork(1)
1777 :
1778 68919 : DEALLOCATE (iwork, rwork, work)
1779 482433 : ALLOCATE (iwork(liwork), rwork(lrwork), work(lwork))
1780 2507197 : eigenvectors(:, :) = matrix(:, :)
1781 :
1782 : ! final diagonalization
1783 68919 : CALL zheevd('V', 'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1784 :
1785 68919 : DEALLOCATE (iwork, rwork, work)
1786 :
1787 68919 : IF (info /= 0) CPABORT("Diagonalisation of a complex matrix failed")
1788 :
1789 68919 : CALL timestop(handle)
1790 :
1791 68919 : 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 3539 : 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 3539 : 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 3539 : n = SIZE(matrix, 1)
1812 21234 : ALLOCATE (matrix_c(n, n), eigenvalues(n))
1813 :
1814 235869 : matrix_c(:, :) = CMPLX(0.0_dp, -matrix, kind=dp)
1815 3539 : CALL diag_complex(matrix_c, evecs, eigenvalues)
1816 27874 : evals = CMPLX(0.0_dp, eigenvalues, kind=dp)
1817 :
1818 3539 : DEALLOCATE (matrix_c, eigenvalues)
1819 3539 : 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 2570 : 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 2570 : n = SIZE(A_in, 1)
1842 2570 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1843 2570 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1844 2570 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1845 2570 : IF (n /= SIZE(C_out, 1)) CPABORT("Incompatible (rows) result array 3 (C).")
1846 2570 : IF (n /= SIZE(C_out, 2)) CPABORT("Incompatible (cols) result array 3 (C).")
1847 2570 : 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 0 : CPABORT("Unknown transpose character for array 1 (A).")
1851 2570 : 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 0 : CPABORT("Unknown transpose character for array 2 (B).")
1855 :
1856 2570 : CALL timeset(routineN, handle)
1857 :
1858 2570 : CALL ZGEMM(A_trans, B_trans, n, n, n, z_one, A_in, n, B_in, n, z_zero, C_out, n)
1859 :
1860 2570 : CALL timestop(handle)
1861 :
1862 2570 : 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 2 : 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 2 : n = SIZE(A_in, 1)
1885 2 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1886 2 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1887 2 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1888 2 : IF (n /= SIZE(C_out, 1)) CPABORT("Incompatible (rows) result array 3 (C).")
1889 2 : IF (n /= SIZE(C_out, 2)) CPABORT("Incompatible (cols) result array 3 (C).")
1890 2 : 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 0 : CPABORT("Unknown transpose character for array 1 (A).")
1894 2 : 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 0 : CPABORT("Unknown transpose character for array 2 (B).")
1898 :
1899 2 : CALL timeset(routineN, handle)
1900 :
1901 2 : 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 2 : CALL timestop(handle)
1904 :
1905 2 : 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 36386 : 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 36386 : COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1930 : INTEGER :: handle, n
1931 :
1932 36386 : n = SIZE(A_in, 1)
1933 36386 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1934 36386 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1935 36386 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1936 36386 : IF (n /= SIZE(C_in, 1)) CPABORT("Incompatible (rows) array 3 (C).")
1937 36386 : IF (n /= SIZE(C_in, 2)) CPABORT("Non-square array 3 (C).")
1938 36386 : IF (n /= SIZE(D_out, 1)) CPABORT("Incompatible (rows) result array 4 (D).")
1939 36386 : IF (n /= SIZE(D_out, 2)) CPABORT("Incompatible (cols) result array 4 (D).")
1940 36386 : 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 0 : CPABORT("Unknown transpose character for array 1 (A).")
1944 36386 : 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 0 : CPABORT("Unknown transpose character for array 2 (B).")
1948 36386 : 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 0 : CPABORT("Unknown transpose character for array 3 (C).")
1952 :
1953 36386 : CALL timeset(routineN, handle)
1954 :
1955 145544 : ALLOCATE (work(n, n), source=z_zero)
1956 :
1957 36386 : CALL ZGEMM(A_trans, B_trans, n, n, n, z_one, A_in, n, B_in, n, z_zero, work, n)
1958 36386 : CALL ZGEMM('N', C_trans, n, n, n, z_one, work, n, C_in, n, z_zero, D_out, n)
1959 :
1960 36386 : DEALLOCATE (work)
1961 :
1962 36386 : CALL timestop(handle)
1963 :
1964 36386 : 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 4 : 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 4 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1990 :
1991 4 : n = SIZE(A_in, 1)
1992 4 : IF (n /= SIZE(A_in, 2)) CPABORT("Non-square array 1 (A).")
1993 4 : IF (n /= SIZE(B_in, 1)) CPABORT("Incompatible (rows) array 2 (B).")
1994 4 : IF (n /= SIZE(B_in, 2)) CPABORT("Non-square array 2 (B).")
1995 4 : IF (n /= SIZE(C_in, 1)) CPABORT("Incompatible (rows) array 3 (C).")
1996 4 : IF (n /= SIZE(C_in, 2)) CPABORT("Non-square array 3 (C).")
1997 4 : IF (n /= SIZE(D_out, 1)) CPABORT("Incompatible (rows) result array 4 (D).")
1998 4 : IF (n /= SIZE(D_out, 2)) CPABORT("Incompatible (cols) result array 4 (D).")
1999 4 : 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 0 : CPABORT("Unknown transpose character for array 1 (A).")
2003 4 : 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 0 : CPABORT("Unknown transpose character for array 2 (B).")
2007 4 : 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 0 : CPABORT("Unknown transpose character for array 3 (C).")
2011 :
2012 4 : CALL timeset(routineN, handle)
2013 :
2014 16 : ALLOCATE (work(n, n), source=0.0_dp)
2015 :
2016 4 : CALL DGEMM(A_trans, B_trans, n, n, n, 1.0_dp, A_in, n, B_in, n, 0.0_dp, work, n)
2017 4 : CALL DGEMM('N', C_trans, n, n, n, 1.0_dp, work, n, C_in, n, 0.0_dp, D_out, n)
2018 :
2019 4 : DEALLOCATE (work)
2020 :
2021 4 : CALL timestop(handle)
2022 :
2023 4 : 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 4 : 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 4 : 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 4 : CALL timeset(routineN, handle)
2049 :
2050 : ! Test sizes
2051 4 : nao = SIZE(A_in, 1)
2052 4 : nmo = SIZE(eigenvalues)
2053 20 : ALLOCATE (evals(nao), cevals(nao))
2054 32 : ALLOCATE (work(nao, nao), B(nao, nao), A(nao, nao))
2055 21988 : A(:, :) = A_in(:, :)
2056 21988 : B(:, :) = B_in(:, :)
2057 :
2058 : ! Diagonalize -S matrix, this way the NULL space is at the end of the spectrum
2059 21988 : B = -B
2060 4 : CALL diag_complex(B, work, evals)
2061 228 : evals(:) = -evals(:)
2062 4 : nc = nao
2063 228 : DO i = 1, nao
2064 228 : IF (evals(i) < -1.0_dp) THEN
2065 0 : nc = i - 1
2066 0 : EXIT
2067 : END IF
2068 : END DO
2069 4 : CPASSERT(nc /= 0)
2070 :
2071 4 : IF (nc /= nao) THEN
2072 0 : IF (nc < nmo) THEN
2073 : ! Copy NULL space definition to last vectors of eigenvectors (if needed)
2074 0 : ncol = nmo - nc
2075 0 : 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 0 : DO icol = nc + 1, nao
2079 0 : DO irow = 1, nao
2080 0 : work(irow, icol) = czero
2081 : END DO
2082 : END DO
2083 : ! Set small eigenvalues to a dummy save value
2084 0 : evals(nc + 1:nao) = 1.0_dp
2085 : END IF
2086 : ! Calculate U*s**(-1/2)
2087 228 : cevals(:) = CMPLX(1.0_dp/SQRT(evals(:)), 0.0_dp, KIND=dp)
2088 228 : DO i = 1, MIN(SIZE(work, 2), SIZE(cevals))
2089 228 : 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 4 : CALL gemm_square(work, 'C', A, 'N', B)
2093 4 : CALL gemm_square(B, 'N', work, 'N', A)
2094 4 : IF (nc /= nao) THEN
2095 : ! set diagonal values to save large value
2096 0 : DO icol = nc + 1, nao
2097 0 : A(icol, icol) = 10000*cone
2098 : END DO
2099 : END IF
2100 : ! Diagonalize
2101 4 : CALL diag_complex(A, B, evals)
2102 228 : eigenvalues(1:nmo) = evals(1:nmo)
2103 4 : nx = MIN(nc, nmo)
2104 : ! Restore vectors C = U^(-1) * C*
2105 4 : lda = SIZE(work, 1)
2106 4 : ldb = SIZE(B, 1)
2107 4 : ldc = SIZE(eigenvectors, 1)
2108 : CALL ZGEMM("N", "N", nao, nx, nc, cone, work, &
2109 4 : lda, B, ldb, czero, eigenvectors, ldc)
2110 :
2111 4 : DEALLOCATE (evals, cevals, work, B, A)
2112 :
2113 4 : CALL timestop(handle)
2114 4 : END SUBROUTINE geeig_right
2115 :
2116 : END MODULE mathlib
|
