22 #include "../base/base_uses.f90"
28 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'mathlib'
29 REAL(KIND=
dp),
PARAMETER :: eps_geo = 1.0e-6_dp
66 MODULE PROCEDURE det_3x3_1, det_3x3_2
69 INTERFACE invert_matrix
70 MODULE PROCEDURE invert_matrix_d, invert_matrix_z
74 MODULE PROCEDURE set_diag_scalar_d, set_diag_scalar_z
78 MODULE PROCEDURE swap_scalar, swap_vector
82 MODULE PROCEDURE unit_matrix_d, unit_matrix_z
98 FUNCTION pswitch(x, a, b, order)
RESULT(fx)
99 REAL(kind=
dp) :: x, a, b
103 REAL(kind=
dp) :: u, u2, u3
106 IF (x < a .OR. x > b)
THEN
127 fx = 1._dp - 10._dp*u3 + 15._dp*u2*u2 - 6._dp*u2*u3
130 fx = -30._dp*u2 + 60._dp*u*u2 - 30._dp*u2*u2
134 fx = -60._dp*u + 180._dp*u2 - 120._dp*u*u2
137 cpabort(
'order not defined')
152 CHARACTER(LEN=32) :: buffer
160 IF (index(buffer,
"N") .NE. 0 .OR. index(buffer,
"n") .NE. 0)
abnormal_value = .true.
174 PURE FUNCTION angle(a, b)
RESULT(angle_ab)
176 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: a, b
177 REAL(kind=
dp) :: angle_ab
179 REAL(kind=
dp) :: length_of_a, length_of_b
180 REAL(kind=
dp),
DIMENSION(SIZE(a, 1)) :: a_norm, b_norm
182 length_of_a = sqrt(dot_product(a, a))
183 length_of_b = sqrt(dot_product(b, b))
185 IF ((length_of_a > eps_geo) .AND. (length_of_b > eps_geo))
THEN
186 a_norm(:) = a(:)/length_of_a
187 b_norm(:) = b(:)/length_of_b
188 angle_ab = acos(min(max(dot_product(a_norm, b_norm), -1.0_dp), 1.0_dp))
207 INTEGER,
INTENT(IN) :: n, k
208 REAL(kind=
dp) :: n_over_k
210 IF ((k >= 0) .AND. (k <= n))
THEN
232 REAL(kind=
dp),
INTENT(IN) :: z
233 INTEGER,
INTENT(IN) :: k
234 REAL(kind=
dp) :: z_over_k
241 z_over_k = z_over_k*(z - i + 1)/real(i,
dp)
258 INTEGER,
INTENT(IN) :: n
259 INTEGER,
DIMENSION(:),
INTENT(IN) :: k
263 REAL(kind=
dp) :: denom
265 IF (all(k >= 0) .AND. sum(k) == n)
THEN
268 denom = denom*
fac(k(i))
289 REAL(kind=
dp),
INTENT(IN) :: phi
290 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: a
291 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(OUT) :: rotmat
293 REAL(kind=
dp) :: cosp, cost, length_of_a, sinp
294 REAL(kind=
dp),
DIMENSION(3) :: d
296 length_of_a = sqrt(a(1)*a(1) + a(2)*a(2) + a(3)*a(3))
298 IF (length_of_a > eps_geo)
THEN
300 d(:) = a(:)/length_of_a
306 rotmat(1, 1) = d(1)*d(1)*cost + cosp
307 rotmat(1, 2) = d(1)*d(2)*cost - d(3)*sinp
308 rotmat(1, 3) = d(1)*d(3)*cost + d(2)*sinp
309 rotmat(2, 1) = d(2)*d(1)*cost + d(3)*sinp
310 rotmat(2, 2) = d(2)*d(2)*cost + cosp
311 rotmat(2, 3) = d(2)*d(3)*cost - d(1)*sinp
312 rotmat(3, 1) = d(3)*d(1)*cost - d(2)*sinp
313 rotmat(3, 2) = d(3)*d(2)*cost + d(1)*sinp
314 rotmat(3, 3) = d(3)*d(3)*cost + cosp
316 CALL unit_matrix(rotmat)
329 PURE FUNCTION det_3x3_1(a)
RESULT(det_a)
330 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(IN) :: a
331 REAL(kind=
dp) :: det_a
333 det_a = a(1, 1)*(a(2, 2)*a(3, 3) - a(2, 3)*a(3, 2)) + &
334 a(1, 2)*(a(2, 3)*a(3, 1) - a(2, 1)*a(3, 3)) + &
335 a(1, 3)*(a(2, 1)*a(3, 2) - a(2, 2)*a(3, 1))
337 END FUNCTION det_3x3_1
349 PURE FUNCTION det_3x3_2(a1, a2, a3)
RESULT(det_a)
350 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: a1, a2, a3
351 REAL(kind=
dp) :: det_a
353 det_a = a1(1)*(a2(2)*a3(3) - a3(2)*a2(3)) + &
354 a2(1)*(a3(2)*a1(3) - a1(2)*a3(3)) + &
355 a3(1)*(a1(2)*a2(3) - a2(2)*a1(3))
357 END FUNCTION det_3x3_2
376 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: a
377 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: eigval
378 LOGICAL,
INTENT(IN),
OPTIONAL :: dac
380 CHARACTER(len=*),
PARAMETER :: routinen =
'diamat_all'
382 INTEGER :: handle, info, liwork, lwork, n, nb
383 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: iwork
384 INTEGER,
EXTERNAL :: ilaenv
385 LOGICAL :: divide_and_conquer
386 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: work
388 EXTERNAL dsyev, dsyevd
390 CALL timeset(routinen, handle)
396 IF (
SIZE(a, 2) /= n)
THEN
397 cpabort(
"Check the size of matrix a (parameter #1)")
401 IF (
SIZE(eigval) /= n)
THEN
402 cpabort(
"The dimension of vector eigval is too small")
407 IF (
PRESENT(dac))
THEN
408 divide_and_conquer = dac
410 divide_and_conquer = .false.
415 IF (divide_and_conquer)
THEN
416 lwork = 2*n**2 + 6*n + 1
419 nb = ilaenv(1,
"DSYTRD",
"U", n, -1, -1, -1)
425 ALLOCATE (work(lwork))
426 IF (divide_and_conquer)
THEN
427 ALLOCATE (iwork(liwork))
433 IF (divide_and_conquer)
THEN
434 CALL dsyevd(
"V",
"U", n, a, n, eigval, work, lwork, iwork, liwork, info)
436 CALL dsyev(
"V",
"U", n, a, n, eigval, work, lwork, info)
440 IF (divide_and_conquer)
THEN
441 cpabort(
"The matrix diagonalization with dsyevd failed")
443 cpabort(
"The matrix diagonalization with dsyev failed")
451 IF (divide_and_conquer)
THEN
455 CALL timestop(handle)
472 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: ab, bc, cd
473 REAL(kind=
dp) :: dihedral_angle_abcd
475 REAL(kind=
dp) :: det_abcd
476 REAL(kind=
dp),
DIMENSION(3) :: abc, bcd
483 det_abcd =
det_3x3(abc, bcd, -bc)
484 dihedral_angle_abcd = sign(1.0_dp, det_abcd)*
angle(abc, bcd)
497 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: a
499 DIMENSION(MIN(SIZE(a, 1), SIZE(a, 2))) :: a_diag
503 n = min(
SIZE(a, 1),
SIZE(a, 2))
521 REAL(kind=
dp),
DIMENSION(3, 3),
INTENT(IN) :: a
522 REAL(kind=
dp),
DIMENSION(3, 3) :: a_inv
524 REAL(kind=
dp) :: det_a
528 a_inv(1, 1) = (a(2, 2)*a(3, 3) - a(3, 2)*a(2, 3))*det_a
529 a_inv(2, 1) = (a(2, 3)*a(3, 1) - a(3, 3)*a(2, 1))*det_a
530 a_inv(3, 1) = (a(2, 1)*a(3, 2) - a(3, 1)*a(2, 2))*det_a
532 a_inv(1, 2) = (a(1, 3)*a(3, 2) - a(3, 3)*a(1, 2))*det_a
533 a_inv(2, 2) = (a(1, 1)*a(3, 3) - a(3, 1)*a(1, 3))*det_a
534 a_inv(3, 2) = (a(1, 2)*a(3, 1) - a(3, 2)*a(1, 1))*det_a
536 a_inv(1, 3) = (a(1, 2)*a(2, 3) - a(2, 2)*a(1, 3))*det_a
537 a_inv(2, 3) = (a(1, 3)*a(2, 1) - a(2, 3)*a(1, 1))*det_a
538 a_inv(3, 3) = (a(1, 1)*a(2, 2) - a(2, 1)*a(1, 2))*det_a
548 REAL(kind=
dp),
INTENT(INOUT) :: a(:, :)
549 INTEGER,
INTENT(OUT) :: info
552 INTEGER,
ALLOCATABLE :: ipiv(:)
553 REAL(kind=
dp),
ALLOCATABLE :: work(:)
558 ALLOCATE (work(lwork))
562 CALL dgetrf(n, n, a, n, ipiv, info)
564 CALL dgetri(n, a, n, ipiv, work, lwork, info)
566 DEALLOCATE (ipiv, work)
579 REAL(kind=
dp),
INTENT(INOUT) :: a(:, :)
580 CHARACTER(LEN=1),
INTENT(IN),
OPTIONAL :: cholesky_triangle
582 CHARACTER(LEN=*),
PARAMETER :: routinen =
'invmat_symm'
584 CHARACTER(LEN=1) :: my_triangle
585 INTEGER :: handle, i, info, n
587 CALL timeset(routinen, handle)
592 IF (
PRESENT(cholesky_triangle))
THEN
593 my_triangle = cholesky_triangle
599 IF (.NOT.
PRESENT(cholesky_triangle))
THEN
600 CALL dpotrf(my_triangle, n, a, n, info)
602 cpabort(
"DPOTRF failed")
607 CALL dpotri(my_triangle, n, a, n, info)
609 cpabort(
"Matrix inversion failed")
612 IF (my_triangle ==
"U")
THEN
614 a(i + 1:n, i) = a(i, i + 1:n)
618 a(i, i + 1:n) = a(i + 1:n, i)
622 CALL timestop(handle)
648 SUBROUTINE invert_matrix_d(a, a_inverse, eval_error, option, improve)
649 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: a
650 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: a_inverse
651 REAL(kind=
dp),
INTENT(OUT) :: eval_error
652 CHARACTER(LEN=1),
INTENT(IN),
OPTIONAL :: option
653 LOGICAL,
INTENT(IN),
OPTIONAL :: improve
655 CHARACTER(LEN=1) :: norm, trans
656 CHARACTER(LEN=default_string_length) :: message
657 INTEGER :: info, iter, n
658 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ipiv, iwork
659 LOGICAL :: do_improve
660 REAL(kind=
dp) :: a_norm, old_eval_error, r_cond
661 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: berr, ferr, work
662 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: a_lu, b
663 REAL(kind=
dp),
EXTERNAL :: dlange
665 EXTERNAL dgecon, dgerfs, dgetrf, dgetrs
668 IF (
PRESENT(option))
THEN
674 IF (
PRESENT(improve))
THEN
685 cpabort(
"Matrix to be inverted of zero size")
688 IF (n /=
SIZE(a, 2))
THEN
689 cpabort(
"Check the array bounds of parameter #1")
692 IF ((n /=
SIZE(a_inverse, 1)) .OR. &
693 (n /=
SIZE(a_inverse, 2)))
THEN
694 cpabort(
"Check the array bounds of parameter #2")
698 ALLOCATE (a_lu(n, n))
712 a_lu(1:n, 1:n) = a(1:n, 1:n)
715 CALL dgetrf(n, n, a_lu, n, ipiv, info)
718 cpabort(
"The LU factorization in dgetrf failed")
723 IF (trans ==
"N")
THEN
729 a_norm = dlange(norm, n, n, a, n, work)
733 CALL dgecon(norm, n, a_lu, n, a_norm, r_cond, work, iwork, info)
736 cpabort(
"The computation of the condition number in dgecon failed")
739 IF (r_cond < epsilon(0.0_dp))
THEN
740 WRITE (message,
"(A,ES10.3)")
"R_COND =", r_cond
741 CALL cp_abort(__location__, &
742 "Bad condition number "//trim(message)//
" (smaller than the machine "// &
743 "working precision)")
748 CALL unit_matrix(a_inverse)
750 CALL dgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
753 cpabort(
"Solving the system of linear equations in dgetrs failed")
764 CALL dgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
768 cpabort(
"Improving the computed solution in dgerfs failed")
771 old_eval_error = eval_error
772 eval_error = maxval(ferr)
774 IF (abs(eval_error - old_eval_error) <= epsilon(1.0_dp))
EXIT
788 END SUBROUTINE invert_matrix_d
810 SUBROUTINE invert_matrix_z(a, a_inverse, eval_error, option)
811 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(IN) :: a
812 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(OUT) :: a_inverse
813 REAL(kind=
dp),
INTENT(OUT) :: eval_error
814 CHARACTER(LEN=1),
INTENT(IN),
OPTIONAL :: option
816 CHARACTER(LEN=1) :: norm, trans
817 CHARACTER(LEN=default_string_length) :: message
818 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:) :: work
819 COMPLEX(KIND=dp),
ALLOCATABLE,
DIMENSION(:, :) :: a_lu, b
820 INTEGER :: info, iter, n
821 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ipiv
822 REAL(kind=
dp) :: a_norm, old_eval_error, r_cond
823 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: berr, ferr, rwork
824 REAL(kind=
dp),
EXTERNAL :: zlange
826 EXTERNAL zgecon, zgerfs, zgetrf, zgetrs
829 IF (
PRESENT(option))
THEN
840 cpabort(
"Matrix to be inverted of zero size")
843 IF (n /=
SIZE(a, 2))
THEN
844 cpabort(
"Check the array bounds of parameter #1")
847 IF ((n /=
SIZE(a_inverse, 1)) .OR. &
848 (n /=
SIZE(a_inverse, 2)))
THEN
849 cpabort(
"Check the array bounds of parameter #2")
853 ALLOCATE (a_lu(n, n))
863 ALLOCATE (rwork(2*n))
867 a_lu(1:n, 1:n) = a(1:n, 1:n)
870 CALL zgetrf(n, n, a_lu, n, ipiv, info)
873 cpabort(
"The LU factorization in dgetrf failed")
878 IF (trans ==
"N")
THEN
884 a_norm = zlange(norm, n, n, a, n, work)
888 CALL zgecon(norm, n, a_lu, n, a_norm, r_cond, work, rwork, info)
891 cpabort(
"The computation of the condition number in dgecon failed")
894 IF (r_cond < epsilon(0.0_dp))
THEN
895 WRITE (message,
"(A,ES10.3)")
"R_COND =", r_cond
896 CALL cp_abort(__location__, &
897 "Bad condition number "//trim(message)//
" (smaller than the machine "// &
898 "working precision)")
903 CALL unit_matrix(a_inverse)
905 CALL zgetrs(trans, n, n, a_lu, n, ipiv, a_inverse, n, info)
908 cpabort(
"Solving the system of linear equations in dgetrs failed")
918 CALL zgerfs(trans, n, n, a, n, a_lu, n, ipiv, b, n, a_inverse, n, ferr, berr, &
922 cpabort(
"Improving the computed solution in dgerfs failed")
925 old_eval_error = eval_error
926 eval_error = maxval(ferr)
928 IF (abs(eval_error - old_eval_error) <= epsilon(1.0_dp))
EXIT
941 END SUBROUTINE invert_matrix_z
955 REAL(kind=
dp),
DIMENSION(:, :) :: a, a_pinverse
956 REAL(kind=
dp),
INTENT(IN) :: rskip
957 REAL(kind=
dp),
INTENT(OUT),
OPTIONAL :: determinant
958 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT), &
959 OPTIONAL,
POINTER :: sval
961 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_pseudo_inverse_svd'
963 INTEGER :: handle, i, info, lwork, n
964 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: iwork
965 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: sig, work
966 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: sig_plus, temp_mat, u, vt
968 CALL timeset(routinen, handle)
971 ALLOCATE (u(n, n), vt(n, n), sig(n), sig_plus(n, n), iwork(8*n), work(1), temp_mat(n, n))
978 IF (
PRESENT(determinant)) determinant = 1.0_dp
982 CALL dgesdd(
'A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
983 lwork, iwork(1), info)
986 cpabort(
"ERROR in DGESDD: Could not retrieve work array sizes")
990 ALLOCATE (work(lwork))
993 CALL dgesdd(
'A', n, n, a(1, 1), n, sig(1), u(1, 1), n, vt(1, 1), n, work(1), &
994 lwork, iwork(1), info)
997 cpabort(
"SVD failed")
1000 IF (
PRESENT(sval))
THEN
1001 cpassert(.NOT.
ASSOCIATED(sval))
1008 IF (sig(i) > rskip*maxval(sig))
THEN
1009 IF (
PRESENT(determinant)) &
1010 determinant = determinant*sig(i)
1011 sig_plus(i, i) = 1._dp/sig(i)
1013 sig_plus(i, i) = 0.0_dp
1018 CALL dgemm(
"N",
"T", n, n, n, 1._dp, sig_plus, n, u, n, 0._dp, temp_mat, n)
1019 CALL dgemm(
"T",
"N", n, n, n, 1._dp, vt, n, temp_mat, n, 0._dp, a_pinverse, n)
1021 DEALLOCATE (u, vt, sig, iwork, work, sig_plus, temp_mat)
1023 CALL timestop(handle)
1037 REAL(kind=
dp),
DIMENSION(:, :) :: a, a_pinverse
1038 REAL(kind=
dp),
INTENT(IN) :: rskip
1040 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_pseudo_inverse_diag'
1042 INTEGER :: handle, i, info, lwork, n
1043 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: eig, work
1044 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: dinv, temp_mat
1046 CALL timeset(routinen, handle)
1050 ALLOCATE (dinv(n, n), eig(n), work(1), temp_mat(n, n))
1058 CALL dsyev(
'V',
'U', n, a, n, eig(1), work(1), lwork, info)
1060 cpabort(
"ERROR in DSYEV: Could not retrieve work array sizes")
1062 lwork = int(work(1))
1064 ALLOCATE (work(lwork))
1068 CALL dsyev(
'V',
'U', n, a, n, eig(1), work(1), lwork, info)
1071 cpabort(
"Matrix diagonalization failed")
1076 IF (eig(i) > rskip*maxval(eig))
THEN
1077 dinv(i, i) = 1.0_dp/eig(i)
1084 CALL dgemm(
"N",
"T", n, n, n, 1._dp, dinv, n, a, n, 0._dp, temp_mat, n)
1085 CALL dgemm(
"N",
"N", n, n, n, 1._dp, a, n, temp_mat, n, 0._dp, a_pinverse, n)
1087 DEALLOCATE (eig, work, dinv, temp_mat)
1089 CALL timestop(handle)
1105 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: a, b
1106 REAL(kind=
dp),
DIMENSION(3) :: a_mirror
1108 REAL(kind=
dp) :: length_of_b, scapro
1109 REAL(kind=
dp),
DIMENSION(3) :: d
1111 length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1113 IF (length_of_b > eps_geo)
THEN
1115 d(:) = b(:)/length_of_b
1118 scapro = a(1)*d(1) + a(2)*d(2) + a(3)*d(3)
1120 a_mirror(:) = a(:) - 2.0_dp*scapro*d(:)
1124 a_mirror(:) = 0.0_dp
1143 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: a
1144 REAL(kind=
dp),
INTENT(IN) :: phi
1145 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: b
1146 REAL(kind=
dp),
DIMENSION(3) :: a_rot
1148 REAL(kind=
dp) :: length_of_b
1149 REAL(kind=
dp),
DIMENSION(3, 3) :: rotmat
1151 length_of_b = sqrt(b(1)*b(1) + b(2)*b(2) + b(3)*b(3))
1152 IF (length_of_b > eps_geo)
THEN
1158 a_rot(:) = matmul(rotmat, a)
1176 PURE SUBROUTINE set_diag_scalar_d(a, b)
1177 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: a
1178 REAL(kind=
dp),
INTENT(IN) :: b
1182 n = min(
SIZE(a, 1),
SIZE(a, 2))
1187 END SUBROUTINE set_diag_scalar_d
1194 PURE SUBROUTINE set_diag_scalar_z(a, b)
1195 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(INOUT) :: a
1196 COMPLEX(KIND=dp),
INTENT(IN) :: b
1200 n = min(
SIZE(a, 1),
SIZE(a, 2))
1205 END SUBROUTINE set_diag_scalar_z
1216 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: a
1217 CHARACTER(LEN=*),
INTENT(IN) :: option
1221 n = min(
SIZE(a, 1),
SIZE(a, 2))
1223 IF (option ==
"lower_to_upper")
THEN
1225 a(i, i + 1:n) = a(i + 1:n, i)
1227 ELSE IF (option ==
"upper_to_lower")
THEN
1229 a(i + 1:n, i) = a(i, i + 1:n)
1231 ELSE IF (option ==
"anti_lower_to_upper")
THEN
1233 a(i, i + 1:n) = -a(i + 1:n, i)
1235 ELSE IF (option ==
"anti_upper_to_lower")
THEN
1237 a(i + 1:n, i) = -a(i, i + 1:n)
1240 cpabort(
"Invalid option <"//trim(option)//
"> was specified for parameter #2")
1252 PURE SUBROUTINE unit_matrix_d(a)
1253 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: a
1256 CALL set_diag(a, 1.0_dp)
1258 END SUBROUTINE unit_matrix_d
1264 PURE SUBROUTINE unit_matrix_z(a)
1265 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(INOUT) :: a
1267 a(:, :) = (0.0_dp, 0.0_dp)
1268 CALL set_diag(a, (1.0_dp, 0.0_dp))
1270 END SUBROUTINE unit_matrix_z
1282 REAL(kind=
dp),
DIMENSION(3),
INTENT(IN) :: a, b
1283 REAL(kind=
dp),
DIMENSION(3) :: c
1285 c(1) = a(2)*b(3) - a(3)*b(2)
1286 c(2) = a(3)*b(1) - a(1)*b(3)
1287 c(3) = a(1)*b(2) - a(2)*b(1)
1299 INTEGER,
INTENT(IN) :: a, b
1302 INTEGER :: aa, ab, l, rem, s
1334 INTEGER,
INTENT(IN) :: a, b
1344 lcm = abs((a/tmp)*b)
1358 INTEGER,
PARAMETER :: maxit = 100
1359 REAL(
dp),
PARAMETER :: eps = epsilon(0.0_dp), &
1360 fpmin = tiny(0.0_dp)
1363 REAL(
dp) :: fact, prev, sum1, term
1365 IF (x <= 0._dp)
THEN
1366 cpabort(
"Invalid argument")
1371 ELSE IF (x <= -log(eps))
THEN
1375 fact = fact*x/real(k,
dp)
1376 term = fact/real(k,
dp)
1378 IF (term < eps*sum1)
EXIT
1380 ei = sum1 + log(x) +
euler
1386 term = term*real(k,
dp)/x
1387 IF (term < eps)
EXIT
1388 IF (term < prev)
THEN
1395 ei = exp(x)*(1._dp + sum1)/x
1412 INTEGER,
INTENT(IN) :: n
1413 REAL(
dp),
INTENT(IN) :: x
1416 INTEGER,
PARAMETER :: maxit = 100
1417 REAL(
dp),
PARAMETER :: eps = 6.e-14_dp,
euler = 0.5772156649015328606065120_dp, &
1418 fpmin = tiny(0.0_dp)
1420 INTEGER :: i, ii, nm1
1421 REAL(
dp) :: a, b, c, d, del, fact, h, psi
1425 IF (n .LT. 0 .OR. x .LT. 0.0_dp .OR. (x .EQ. 0.0_dp .AND. (n .EQ. 0 .OR. n .EQ. 1)))
THEN
1426 cpabort(
"Invalid argument")
1427 ELSE IF (n .EQ. 0)
THEN
1429 ELSE IF (x .EQ. 0.0_dp)
THEN
1431 ELSE IF (x .GT. 1.0_dp)
THEN
1439 d = 1.0_dp/(a*d + b)
1443 IF (abs(del - 1.0_dp) .LT. eps)
THEN
1448 cpabort(
"continued fraction failed in expint")
1450 IF (nm1 .NE. 0)
THEN
1458 IF (i .NE. nm1)
THEN
1459 del = -fact/(i - nm1)
1463 psi = psi + 1.0_dp/ii
1465 del = fact*(-log(x) + psi)
1468 IF (abs(del) .LT. abs(
expint)*eps)
RETURN
1470 cpabort(
"series failed in expint")
1487 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: a
1488 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: d
1489 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: v
1496 CALL diag(n, a, d, v)
1499 CALL eigsrt(n, d, v)
1516 INTEGER,
INTENT(IN) :: n
1517 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: a
1518 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: d
1519 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: v
1521 REAL(kind=
dp),
PARAMETER :: a_eps = 1.0e-10_dp, d_eps = 1.0e-3_dp
1523 INTEGER :: i, ip, iq
1524 REAL(kind=
dp) :: a_max, apq, c, d_min, dip, diq, g, h, s, &
1525 t, tau, theta, tresh
1526 REAL(kind=
dp),
DIMENSION(n) :: b, z
1530 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1540 d_min = max(d_eps, minval(abs(b)))
1541 IF (a_max < a_eps*d_min)
RETURN
1542 tresh = merge(a_max, 0.0_dp, (i < 4))
1549 g = 100.0_dp*abs(apq)
1550 IF (tresh < abs(apq))
THEN
1552 IF ((abs(h) + g) .NE. abs(h))
THEN
1553 theta = 0.5_dp*h/apq
1554 t = 1.0_dp/(abs(theta) + sqrt(1.0_dp + theta**2))
1555 IF (theta < 0.0_dp) t = -t
1559 c = 1.0_dp/sqrt(1.0_dp + t**2)
1561 tau = s/(1.0_dp + c)
1568 CALL jrotate(a(1:ip - 1, ip), a(1:ip - 1, iq), s, tau)
1569 CALL jrotate(a(ip, ip + 1:iq - 1), a(ip + 1:iq - 1, iq), s, tau)
1570 CALL jrotate(a(ip, iq + 1:n), a(iq, iq + 1:n), s, tau)
1571 CALL jrotate(v(:, ip), v(:, iq), s, tau)
1572 ELSE IF ((4 < i) .AND. &
1573 ((abs(dip) + g) == abs(dip)) .AND. &
1574 ((abs(diq) + g) == abs(diq)))
THEN
1582 a_max = max(a_max, maxval(abs(a(ip, ip + 1:n))))
1585 WRITE (*,
'(/,T2,A,/)')
'Too many iterations in jacobi'
1599 PURE SUBROUTINE jrotate(a, b, ss, tt)
1600 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: a, b
1601 REAL(kind=
dp),
INTENT(IN) :: ss, tt
1603 REAL(kind=
dp) :: u, v
1609 b = b*(u + ss*v) + a*v
1611 END SUBROUTINE jrotate
1623 SUBROUTINE eigsrt(n, d, v)
1624 INTEGER,
INTENT(IN) :: n
1625 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: d
1626 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: v
1631 j = sum(minloc(d(i:n))) + i - 1
1633 CALL swap(d(i), d(j))
1634 CALL swap(v(:, i), v(:, j))
1638 END SUBROUTINE eigsrt
1648 ELEMENTAL SUBROUTINE swap_scalar(a, b)
1649 REAL(kind=
dp),
INTENT(INOUT) :: a, b
1657 END SUBROUTINE swap_scalar
1667 SUBROUTINE swap_vector(a, b)
1668 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: a, b
1675 IF (n /=
SIZE(b))
THEN
1676 cpabort(
"Check the array bounds of the parameters")
1685 END SUBROUTINE swap_vector
1700 REAL(
dp),
INTENT(in) :: eps, omg
1701 REAL(
dp),
INTENT(out) :: r_cutoff
1703 CHARACTER(LEN=*),
PARAMETER :: routinen =
'erfc_cutoff'
1705 REAL(
dp),
PARAMETER :: abstol = 1e-10_dp, soltol = 1e-16_dp
1706 REAL(
dp) :: r0, f0, fprime0, delta_r
1707 INTEGER :: iter, handle
1708 INTEGER,
PARAMETER :: itermax = 1000
1710 CALL timeset(routinen, handle)
1713 r0 = sqrt(-log(eps*omg*10**2))/omg
1714 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1716 DO iter = 1, itermax
1717 delta_r = f0/fprime0
1719 CALL eval_transc_func(r0, eps, omg, f0, fprime0)
1720 IF (abs(delta_r) .LT. abstol .OR. abs(f0) .LT. soltol)
EXIT
1722 cpassert(iter <= itermax)
1725 CALL timestop(handle)
1735 ELEMENTAL SUBROUTINE eval_transc_func(r, eps, omega, fn, df)
1736 REAL(
dp),
INTENT(in) :: r, eps, omega
1737 REAL(
dp),
INTENT(out) :: fn, df
1739 fn = erfc(omega*r) - r*eps
1740 df = -omega*2*exp(-omega**2*r**2)/sqrt(
pi) - eps
1741 END SUBROUTINE eval_transc_func
1753 COMPLEX(KIND=dp),
DIMENSION(:, :),
INTENT(INOUT) :: matrix, eigenvectors
1754 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: eigenvalues
1756 COMPLEX(KIND=dp),
DIMENSION(:),
ALLOCATABLE :: work
1757 INTEGER :: info, liwork, lrwork, lwork, n
1758 INTEGER,
DIMENSION(:),
ALLOCATABLE :: iwork
1759 REAL(kind=
dp),
DIMENSION(:),
ALLOCATABLE :: rwork
1761 eigenvectors(:, :) = matrix(:, :)
1764 ALLOCATE (iwork(1), rwork(1), work(1))
1771 CALL zheevd(
'V',
'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1773 lwork = ceiling(real(work(1), kind=
dp))
1774 lrwork = ceiling(real(rwork(1), kind=
dp))
1777 DEALLOCATE (iwork, rwork, work)
1779 ALLOCATE (iwork(liwork))
1781 ALLOCATE (rwork(lrwork))
1783 ALLOCATE (work(lwork))
1784 work(:) = cmplx(0.0_dp, 0.0_dp, kind=
dp)
1787 CALL zheevd(
'V',
'U', n, eigenvectors, n, eigenvalues, work, lwork, rwork, lrwork, iwork, liwork, info)
1790 cpabort(
"Diagonalisation of a complex matrix failed")
1791 DEALLOCATE (iwork, rwork, work)
real(kind=dp) function det_3x3(a)
...
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.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public euler
real(kind=dp), dimension(0:maxfac), parameter, public fac
Collection of simple mathematical functions and subroutines.
subroutine, public get_pseudo_inverse_svd(a, a_pinverse, rskip, determinant, sval)
returns the pseudoinverse of a real, square matrix using singular value decomposition
subroutine, public jacobi(a, d, v)
Jacobi matrix diagonalization. The eigenvalues are returned in vector d and the eigenvectors are retu...
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.
elemental integer function, public lcm(a, b)
computes the least common multiplier of two numbers
subroutine, public get_pseudo_inverse_diag(a, a_pinverse, rskip)
returns the pseudoinverse of a real, symmetric and positive definite matrix using diagonalization.
subroutine, public invmat_symm(a, cholesky_triangle)
returns inverse of real symmetric, positive definite matrix
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)*....
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.
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...
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.
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...
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
subroutine, public symmetrize_matrix(a, option)
Symmetrize the matrix a.
subroutine, public complex_diag(matrix, eigenvectors, eigenvalues)
Diagonalizes a local complex Hermitian matrix using LAPACK. Based on cp_cfm_heevd.
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
subroutine, public invmat(a, info)
returns inverse of matrix using the lapack routines DGETRF and DGETRI
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...
subroutine, public diag(n, a, d, v)
Diagonalize matrix a. The eigenvalues are returned in vector d and the eigenvectors are returned in m...
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...
subroutine, public erfc_cutoff(eps, omg, r_cutoff)
compute a truncation radius for the shortrange operator
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...
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,...
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,...
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) ...
pure real(kind=dp) function, dimension(3), public vector_product(a, b)
Calculation of the vector product c = a x b.
pure real(kind=dp) function, public multinomial(n, k)
Calculates the multinomial coefficients.