19 USE libxsmm,
ONLY: libxsmm_prefetch_none, &
20 libxsmm_blasint_kind, &
28 #include "../base/base_uses.f90"
34 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ai_contraction_sphi'
51 SUBROUTINE ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
53 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: abint
54 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: sab, sphi_a, sphi_b
55 INTEGER,
INTENT(IN) :: ncoa, ncob, nsgfa, nsgfb
57 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ab_contract'
59 INTEGER :: handle, m1, m2, msphia, msphib, nn
60 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: cpp
62 CALL timeset(routinen, handle)
64 msphia =
SIZE(sphi_a, 1)
65 msphib =
SIZE(sphi_b, 1)
72 ALLOCATE (cpp(nsgfa, m2))
74 CALL dgemm(
"T",
"N", nsgfa, m2, ncoa, 1._dp, sphi_a, msphia, sab, m1, 0.0_dp, cpp, nsgfa)
75 CALL dgemm(
"N",
"N", nsgfa, nsgfb, ncob, 1._dp, cpp, nsgfa, sphi_b, msphib, 0.0_dp, &
99 SUBROUTINE abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, &
102 REAL(kind=
dp),
DIMENSION(:, :, :) :: abcint, sabc
103 REAL(kind=
dp),
DIMENSION(:, :) :: sphi_a, sphi_b, sphi_c
104 INTEGER,
INTENT(IN) :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
106 CHARACTER(LEN=*),
PARAMETER :: routinen =
'abc_contract'
108 INTEGER :: handle, i, m1, m2, m3, msphia, msphib, &
110 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: cpc, cpp
112 CALL timeset(routinen, handle)
114 cpassert(
SIZE(abcint, 1) == nsgfa)
115 cpassert(
SIZE(abcint, 2) == nsgfb)
117 msphia =
SIZE(sphi_a, 1)
118 msphib =
SIZE(sphi_b, 1)
119 msphic =
SIZE(sphi_c, 1)
125 ALLOCATE (cpp(nsgfa, m2, m3), cpc(nsgfa, nsgfb, m3))
127 CALL dgemm(
"T",
"N", nsgfa, m2*m3, ncoa, 1._dp, sphi_a, msphia, sabc, m1, 0.0_dp, cpp, nsgfa)
129 CALL dgemm(
"N",
"N", nsgfa, nsgfb, ncob, 1._dp, cpp(:, :, i), nsgfa, sphi_b, msphib, &
130 0.0_dp, cpc(:, :, i), nsgfa)
132 CALL dgemm(
"N",
"N", nsgfa*nsgfb, nsgfc, ncoc, 1._dp, cpc, nsgfa*nsgfb, sphi_c, msphic, 0.0_dp, &
135 DEALLOCATE (cpp, cpc)
137 CALL timestop(handle)
159 SUBROUTINE abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, ncob, &
160 ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd)
162 REAL(kind=
dp),
DIMENSION(:, :, :, :), &
163 INTENT(INOUT) :: abcdint
164 REAL(kind=
dp),
DIMENSION(:, :, :, :),
INTENT(IN) :: sabcd
165 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: sphi_a, sphi_b, sphi_c, sphi_d
166 INTEGER,
INTENT(IN) :: ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, &
169 CHARACTER(LEN=*),
PARAMETER :: routinen =
'abcd_contract'
171 INTEGER :: handle, isgfc, isgfd, m1, m2, m3, m4, &
172 msphia, msphib, msphic, msphid
173 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: temp_cccc, work_cpcc
174 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: temp_cpcc, work_cppc
175 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: cpcc, cppc, cppp
177 CALL timeset(routinen, handle)
179 msphia =
SIZE(sphi_a, 1)
180 msphib =
SIZE(sphi_b, 1)
181 msphic =
SIZE(sphi_c, 1)
182 msphid =
SIZE(sphi_d, 1)
189 ALLOCATE (cppp(nsgfa, m2, m3, m4), cppc(nsgfa, m2, m3, nsgfd), &
190 cpcc(nsgfa, m2, nsgfc, nsgfd))
192 ALLOCATE (work_cppc(nsgfa, m2, m3), temp_cpcc(nsgfa, m2, nsgfc))
196 ALLOCATE (work_cpcc(nsgfa, m2), temp_cccc(nsgfa, nsgfb))
200 CALL dgemm(
"T",
"N", nsgfa, m2*m3*m4, ncoa, 1._dp, sphi_a, msphia, sabcd, m1, &
202 CALL dgemm(
"N",
"N", nsgfa*m2*m3, nsgfd, ncod, 1._dp, cppp, nsgfa*m2*m3, &
203 sphi_d, msphid, 0.0_dp, cppc, nsgfa*m2*m3)
206 work_cppc(:, :, :) = cppc(:, :, :, isgfd)
207 CALL dgemm(
"N",
"N", nsgfa*m2, nsgfc, ncoc, 1._dp, work_cppc, nsgfa*m2, &
208 sphi_c, msphic, 0.0_dp, temp_cpcc, nsgfa*m2)
209 cpcc(:, :, :, isgfd) = temp_cpcc(:, :, :)
211 work_cpcc(:, :) = cpcc(:, :, isgfc, isgfd)
212 CALL dgemm(
"N",
"N", nsgfa, nsgfb, ncob, 1._dp, work_cpcc, nsgfa, sphi_b, &
213 msphib, 0.0_dp, temp_cccc, nsgfa)
214 abcdint(:, :, isgfc, isgfd) = temp_cccc(:, :)
218 DEALLOCATE (cpcc, cppc, cppp)
219 DEALLOCATE (work_cpcc, work_cppc, temp_cpcc, temp_cccc)
221 CALL timestop(handle)
246 nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer)
248 REAL(
dp),
DIMENSION(*) :: abcint, sabc
249 REAL(
dp),
DIMENSION(:, :) :: tsphi_a, sphi_b, sphi_c
250 INTEGER,
INTENT(IN) :: ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc
251 REAL(
dp),
DIMENSION(nsgfa*ncob) :: cpp_buffer
252 REAL(
dp),
DIMENSION(nsgfa*nsgfb*ncoc) :: ccp_buffer
254 CHARACTER(LEN=*),
PARAMETER :: routinen =
'libxsmm_abc_contract'
257 LOGICAL :: libxsmm_kernels_available
259 INTEGER(libxsmm_blasint_kind) :: m, n, k
260 TYPE(libxsmm_dmmfunction) :: xmm1, xmm2
263 CALL timeset(routinen, handle)
264 libxsmm_kernels_available = .false.
271 m = nsgfa; n = ncob; k = ncoa
272 CALL libxsmm_dispatch(xmm1, m, n, k, beta=0.0_dp, prefetch=libxsmm_prefetch_none)
273 m = nsgfa; n = nsgfb; k = ncob
274 CALL libxsmm_dispatch(xmm2, m, n, k, beta=0.0_dp, prefetch=libxsmm_prefetch_none)
276 libxsmm_kernels_available = libxsmm_available(xmm1) .AND. libxsmm_available(xmm2)
278 IF (libxsmm_kernels_available)
THEN
281 CALL libxsmm_dmmcall(xmm1, tsphi_a, sabc((i - 1)*ncoa*ncob + 1), cpp_buffer)
282 CALL libxsmm_dmmcall(xmm2, cpp_buffer, sphi_b, ccp_buffer((i - 1)*nsgfa*nsgfb + 1))
285 cpwarn(
"libxsmm available, but kernels are not, fallback to dgemm")
289 IF (.NOT. libxsmm_kernels_available)
THEN
293 CALL dgemm(
"N",
"N", nsgfa, ncob, ncoa, 1.0_dp, tsphi_a, nsgfa, sabc((i - 1)*ncoa*ncob + 1), &
294 ncoa, 0.0_dp, cpp_buffer, nsgfa)
295 CALL dgemm(
"N",
"N", nsgfa, nsgfb, ncob, 1.0_dp, cpp_buffer, nsgfa, sphi_b, ncob, 0.0_dp, &
296 ccp_buffer((i - 1)*nsgfa*nsgfb + 1), nsgfa)
301 CALL dgemm(
"N",
"N", nsgfa*nsgfb, nsgfc, ncoc, 1.0_dp, ccp_buffer, nsgfa*nsgfb, sphi_c, ncoc, &
302 0.0_dp, abcint, nsgfa*nsgfb)
304 CALL timestop(handle)
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.
Contraction of integrals over primitive Cartesian Gaussians based on the contraction matrix sphi whic...
subroutine, public abc_contract(abcint, sabc, sphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc)
contract three-center overlap integrals (a,b,c) and transfer to spherical Gaussians
subroutine, public abcd_contract(abcdint, sabcd, sphi_a, sphi_b, sphi_c, sphi_d, ncoa, ncob, ncoc, ncod, nsgfa, nsgfb, nsgfc, nsgfd)
contract four-center overlap integrals (a,b,c,d) and transfer to spherical Gaussians
subroutine, public libxsmm_abc_contract(abcint, sabc, tsphi_a, sphi_b, sphi_c, ncoa, ncob, ncoc, nsgfa, nsgfb, nsgfc, cpp_buffer, ccp_buffer)
3-center contraction routine from primitive cartesain Gaussians to spherical Gaussian functions....
subroutine, public ab_contract(abint, sab, sphi_a, sphi_b, ncoa, ncob, nsgfa, nsgfb)
contract overlap integrals (a,b) and transfer to spherical Gaussians
Defines the basic variable types.
integer, parameter, public dp