21 #include "../base/base_uses.f90"
25 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'ai_contraction'
29 PUBLIC :: contraction, decontraction, block_add, force_trace
32 MODULE PROCEDURE contraction_ab, contraction_abc
35 INTERFACE decontraction
36 MODULE PROCEDURE decontraction_ab
40 MODULE PROCEDURE force_trace_ab
44 MODULE PROCEDURE block_add_ab
71 SUBROUTINE contraction_ab(sab, qab, ca, na, ma, cb, nb, mb, fscale, trans)
73 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN) :: sab
74 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: qab
75 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN), &
77 INTEGER,
INTENT(IN),
OPTIONAL :: na, ma
78 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN), &
80 INTEGER,
INTENT(IN),
OPTIONAL :: nb, mb
81 REAL(KIND=
dp),
INTENT(IN),
OPTIONAL :: fscale
82 LOGICAL,
INTENT(IN),
OPTIONAL :: trans
84 INTEGER :: lda, ldb, ldq, lds, ldw, mal, mbl, nal, &
88 REAL(KIND=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: work
92 IF (
PRESENT(trans))
THEN
99 IF (
PRESENT(fscale))
THEN
106 IF (
PRESENT(ca))
THEN
107 IF (
PRESENT(na))
THEN
112 IF (
PRESENT(ma))
THEN
119 IF (
PRESENT(cb))
THEN
120 IF (
PRESENT(nb))
THEN
125 IF (
PRESENT(mb))
THEN
136 IF (
PRESENT(ca) .AND.
PRESENT(cb))
THEN
138 ALLOCATE (work(nal, mbl))
141 work(1:nal, 1:mbl) = matmul(sab(1:nal, 1:nbl), cb(1:nbl, 1:mbl))
144 qab(1:mbl, 1:mal) = fs*matmul(transpose(work(1:nal, 1:mbl)), ca(1:nal, 1:mal))
147 qab(1:mal, 1:mbl) = fs*matmul(transpose(ca(1:nal, 1:mal)), work(1:nal, 1:mbl))
150 ELSE IF (
PRESENT(ca))
THEN
151 IF (
PRESENT(nb))
THEN
158 qab(1:nbl, 1:mal) = fs*matmul(transpose(sab(1:nal, 1:nbl)), ca(1:nal, 1:mal))
161 qab(1:mal, 1:nbl) = fs*matmul(transpose(ca(1:nal, 1:mal)), sab(1:nal, 1:nbl))
163 ELSE IF (
PRESENT(cb))
THEN
164 IF (
PRESENT(na))
THEN
171 qab(1:nal, 1:mbl) = fs*matmul(sab(1:nal, 1:nbl), cb(1:nbl, 1:mbl))
174 qab(1:mbl, 1:nal) = fs*matmul(transpose(cb(1:nbl, 1:mbl)), transpose(sab(1:nal, 1:nbl)))
178 cpabort(
"Copy of arrays is not covered here")
181 END SUBROUTINE contraction_ab
200 SUBROUTINE contraction_abc(sabc, qabc, ca, na, ma, cb, nb, mb, cc, nc, mc)
202 REAL(KIND=
dp),
DIMENSION(:, :, :),
INTENT(IN) :: sabc
203 REAL(KIND=
dp),
DIMENSION(:, :, :),
INTENT(INOUT) :: qabc
204 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN), &
206 INTEGER,
INTENT(IN),
OPTIONAL :: na, ma
207 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN), &
209 INTEGER,
INTENT(IN),
OPTIONAL :: nb, mb
210 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN), &
212 INTEGER,
INTENT(IN),
OPTIONAL :: nc, mc
214 INTEGER :: lda, ldb, ldc, mal, mbl, mcl, nal, nbl, &
216 REAL(KIND=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: work1, work2, work3, work4
220 IF (
PRESENT(ca))
THEN
221 IF (
PRESENT(na))
THEN
226 IF (
PRESENT(ma))
THEN
233 IF (
PRESENT(cb))
THEN
234 IF (
PRESENT(nb))
THEN
239 IF (
PRESENT(mb))
THEN
246 IF (
PRESENT(cc))
THEN
247 IF (
PRESENT(nc))
THEN
252 IF (
PRESENT(mc))
THEN
260 IF (
PRESENT(ca) .AND.
PRESENT(cb) .AND.
PRESENT(cc))
THEN
262 ALLOCATE (work1(nal, nbl, ncl))
264 work1(1:nal, 1:nbl, 1:ncl) = sabc(1:nal, 1:nbl, 1:ncl)
266 ALLOCATE (work2(nbl, ncl, mal))
267 CALL dgemm(
"T",
"N", nbl*ncl, mal, nal, 1.0_dp, work1(1, 1, 1), nal, ca(1, 1), lda, &
268 0.0_dp, work2(1, 1, 1), nbl*ncl)
270 ALLOCATE (work3(ncl, mal, mbl))
271 CALL dgemm(
"T",
"N", ncl*mal, mbl, nbl, 1.0_dp, work2(1, 1, 1), nbl, cb(1, 1), ldb, &
272 0.0_dp, work3(1, 1, 1), ncl*mal)
274 ALLOCATE (work4(mal, mbl, mcl))
275 CALL dgemm(
"T",
"N", mal*mbl, mcl, ncl, 1.0_dp, work3(1, 1, 1), ncl, cc(1, 1), ldc, &
276 0.0_dp, work4(1, 1, 1), mal*mbl)
278 qabc(1:mal, 1:mbl, 1:mcl) = work4(1:mal, 1:mbl, 1:mcl)
280 DEALLOCATE (work1, work2, work3, work4)
282 ELSE IF (
PRESENT(ca) .AND.
PRESENT(cb))
THEN
283 cpabort(
"Not implemented")
284 ELSE IF (
PRESENT(ca) .AND.
PRESENT(cc))
THEN
285 cpabort(
"Not implemented")
286 ELSE IF (
PRESENT(cb) .AND.
PRESENT(cc))
THEN
287 cpabort(
"Not implemented")
288 ELSE IF (
PRESENT(ca))
THEN
289 cpabort(
"Not implemented")
290 ELSE IF (
PRESENT(cb))
THEN
291 cpabort(
"Not implemented")
292 ELSE IF (
PRESENT(cc))
THEN
293 cpabort(
"Not implemented")
296 cpabort(
"Copy of arrays is not covered here")
299 END SUBROUTINE contraction_abc
316 SUBROUTINE decontraction_ab(sab, qab, ca, na, ma, cb, nb, mb, trans)
318 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN) :: sab
319 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: qab
320 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN) :: ca
321 INTEGER,
INTENT(IN) :: na, ma
322 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN) :: cb
323 INTEGER,
INTENT(IN) :: nb, mb
324 LOGICAL,
INTENT(IN),
OPTIONAL :: trans
326 INTEGER :: lda, ldb, ldq, lds, ldw
328 REAL(KIND=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: work
331 IF (
PRESENT(trans))
THEN
342 ALLOCATE (work(na, mb))
347 work(1:na, 1:mb) = matmul(ca(1:na, 1:ma), transpose(sab(1:mb, 1:ma)))
350 work(1:na, 1:mb) = matmul(ca(1:na, 1:ma), sab(1:ma, 1:mb))
353 qab(1:na, 1:nb) = matmul(work(1:na, 1:mb), transpose(cb(1:nb, 1:mb)))
357 END SUBROUTINE decontraction_ab
370 SUBROUTINE force_trace_ab(force, sab, pab, na, nb, m, trans)
372 REAL(KIND=
dp),
DIMENSION(:),
INTENT(INOUT) :: force
373 REAL(KIND=
dp),
DIMENSION(:, :, :),
INTENT(IN) :: sab
374 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(IN) :: pab
375 INTEGER,
INTENT(IN) :: na, nb, m
376 LOGICAL,
INTENT(IN),
OPTIONAL :: trans
381 cpassert(m <=
SIZE(sab, 3))
382 cpassert(m <=
SIZE(force, 1))
385 IF (
PRESENT(trans))
THEN
393 force(i) = sum(sab(1:nb, 1:na, i)*pab(1:nb, 1:na))
395 force(i) = sum(sab(1:na, 1:nb, i)*pab(1:na, 1:nb))
399 END SUBROUTINE force_trace_ab
415 SUBROUTINE block_add_ab(dir, sab, na, nb, qab, ia, ib, trans)
417 CHARACTER(LEN=*),
INTENT(IN) :: dir
418 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: sab
419 INTEGER,
INTENT(IN) :: na, nb
420 REAL(KIND=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: qab
421 INTEGER,
INTENT(IN) :: ia, ib
422 LOGICAL,
INTENT(IN),
OPTIONAL :: trans
427 IF (
PRESENT(trans))
THEN
433 IF (dir ==
"IN" .OR. dir ==
"in")
THEN
438 qab(ib:jb, ia:ja) = qab(ib:jb, ia:ja) + sab(1:nb, 1:na)
440 qab(ia:ja, ib:jb) = qab(ia:ja, ib:jb) + sab(1:na, 1:nb)
442 ELSEIF (dir ==
"OUT" .OR. dir ==
"out")
THEN
447 sab(1:nb, 1:na) = sab(1:nb, 1:na) + qab(ib:jb, ia:ja)
449 sab(1:na, 1:nb) = sab(1:na, 1:nb) + qab(ia:ja, ib:jb)
455 END SUBROUTINE block_add_ab
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.
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Defines the basic variable types.
integer, parameter, public dp