(git:e68414f)
Loading...
Searching...
No Matches
ai_contraction.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Set of routines to:
10!> Contract integrals over primitive Gaussians
11!> Decontract (density) matrices
12!> Trace matrices to get forces
13!> Block copy and add matrices
14!> \par History
15!> Replace dgemm by MATMUL: Massive speedups in openMP loops (JGH, 12.2019)
16!> \author JGH (01.07.2014)
17! **************************************************************************************************
18MODULE ai_contraction
19
20 USE kinds, ONLY: dp
21#include "../base/base_uses.f90"
22
23 IMPLICIT NONE
24
25 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ai_contraction'
26
27 PRIVATE
28
29 PUBLIC :: contraction, decontraction, block_add, force_trace
30
31 INTERFACE contraction
32 MODULE PROCEDURE contraction_ab, contraction_abc
33 END INTERFACE
34
35 INTERFACE decontraction
36 MODULE PROCEDURE decontraction_ab
37 END INTERFACE
38
39 INTERFACE force_trace
40 MODULE PROCEDURE force_trace_ab
41 END INTERFACE
42
43 INTERFACE block_add
44 MODULE PROCEDURE block_add_ab
45 END INTERFACE
46
47! **************************************************************************************************
48
49CONTAINS
50
51! **************************************************************************************************
52!> \brief Applying the contraction coefficients to a set of two-center primitive
53!> integrals
54!> QAB <- CA(T) * SAB * CB
55!> QAB is optionally scaled with "fscale"
56!> Variable "trans" requests the output to be QAB(T)
57!> If only one of the transformation matrix is given, only a half
58!> transformation is done
59!> Active dimensions are: QAB(ma,mb), SAB(na,nb)
60!> \param sab Input matrix, dimension(:,:)
61!> \param qab Output matrix, dimension(:,:)
62!> \param ca Left transformation matrix, optional
63!> \param na First dimension of ca, optional
64!> \param ma Second dimension of ca, optional
65!> \param cb Right transformation matrix, optional
66!> \param nb First dimension of cb, optional
67!> \param mb Second dimension of cb, optional
68!> \param fscale Optional scaling of output
69!> \param trans Optional transposition of output
70! **************************************************************************************************
71 SUBROUTINE contraction_ab(sab, qab, ca, na, ma, cb, nb, mb, fscale, trans)
72
73 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: sab
74 REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: qab
75 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
76 OPTIONAL :: ca
77 INTEGER, INTENT(IN), OPTIONAL :: na, ma
78 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
79 OPTIONAL :: cb
80 INTEGER, INTENT(IN), OPTIONAL :: nb, mb
81 REAL(KIND=dp), INTENT(IN), OPTIONAL :: fscale
82 LOGICAL, INTENT(IN), OPTIONAL :: trans
83
84 INTEGER :: lda, ldb, ldq, lds, ldw, mal, mbl, nal, &
85 nbl
86 LOGICAL :: my_trans
87 REAL(KIND=dp) :: fs
88 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
89
90! Should output matrix be transposed?
91
92 IF (PRESENT(trans)) THEN
93 my_trans = trans
94 ELSE
95 my_trans = .false.
96 END IF
97
98 ! Scaling of output matrix
99 IF (PRESENT(fscale)) THEN
100 fs = fscale
101 ELSE
102 fs = 1.0_dp
103 END IF
104
105 ! Active matrix size
106 IF (PRESENT(ca)) THEN
107 IF (PRESENT(na)) THEN
108 nal = na
109 ELSE
110 nal = SIZE(ca, 1)
111 END IF
112 IF (PRESENT(ma)) THEN
113 mal = ma
114 ELSE
115 mal = SIZE(ca, 2)
116 END IF
117 lda = SIZE(ca, 1)
118 END IF
119 IF (PRESENT(cb)) THEN
120 IF (PRESENT(nb)) THEN
121 nbl = nb
122 ELSE
123 nbl = SIZE(cb, 1)
124 END IF
125 IF (PRESENT(mb)) THEN
126 mbl = mb
127 ELSE
128 mbl = SIZE(cb, 2)
129 END IF
130 ldb = SIZE(cb, 1)
131 END IF
132
133 lds = SIZE(sab, 1)
134 ldq = SIZE(qab, 1)
135
136 IF (PRESENT(ca) .AND. PRESENT(cb)) THEN
137 ! Full transform
138 ALLOCATE (work(nal, mbl))
139 ldw = nal
140!dg CALL dgemm("N", "N", nal, mbl, nbl, 1.0_dp, sab(1, 1), lds, cb(1, 1), ldb, 0.0_dp, work(1, 1), ldw)
141 work(1:nal, 1:mbl) = matmul(sab(1:nal, 1:nbl), cb(1:nbl, 1:mbl))
142 IF (my_trans) THEN
143!dg CALL dgemm("T", "N", mbl, mal, nal, fs, work(1, 1), ldw, ca(1, 1), lda, 0.0_dp, qab(1, 1), ldq)
144 qab(1:mbl, 1:mal) = fs*matmul(transpose(work(1:nal, 1:mbl)), ca(1:nal, 1:mal))
145 ELSE
146!dg CALL dgemm("T", "N", mal, mbl, nal, fs, ca(1, 1), lda, work(1, 1), ldw, 0.0_dp, qab(1, 1), ldq)
147 qab(1:mal, 1:mbl) = fs*matmul(transpose(ca(1:nal, 1:mal)), work(1:nal, 1:mbl))
148 END IF
149 DEALLOCATE (work)
150 ELSE IF (PRESENT(ca)) THEN
151 IF (PRESENT(nb)) THEN
152 nbl = nb
153 ELSE
154 nbl = SIZE(sab, 2)
155 END IF
156 IF (my_trans) THEN
157!dg CALL dgemm("T", "N", nbl, mal, nal, fs, sab(1, 1), lds, ca(1, 1), lda, 0.0_dp, qab(1, 1), ldq)
158 qab(1:nbl, 1:mal) = fs*matmul(transpose(sab(1:nal, 1:nbl)), ca(1:nal, 1:mal))
159 ELSE
160!dg CALL dgemm("T", "N", mal, nbl, nal, fs, ca(1, 1), lda, sab(1, 1), lds, 0.0_dp, qab(1, 1), ldq)
161 qab(1:mal, 1:nbl) = fs*matmul(transpose(ca(1:nal, 1:mal)), sab(1:nal, 1:nbl))
162 END IF
163 ELSE IF (PRESENT(cb)) THEN
164 IF (PRESENT(na)) THEN
165 nal = na
166 ELSE
167 nal = SIZE(sab, 1)
168 END IF
169 IF (my_trans) THEN
170!dg CALL dgemm("N", "N", nal, mbl, nbl, fs, sab(1, 1), lds, cb(1, 1), ldb, 0.0_dp, qab, ldq)
171 qab(1:nal, 1:mbl) = fs*matmul(sab(1:nal, 1:nbl), cb(1:nbl, 1:mbl))
172 ELSE
173!dg CALL dgemm("T", "T", mbl, nal, nbl, fs, cb(1, 1), ldb, sab(1, 1), lds, 0.0_dp, qab, ldq)
174 qab(1:mbl, 1:nal) = fs*matmul(transpose(cb(1:nbl, 1:mbl)), transpose(sab(1:nal, 1:nbl)))
175 END IF
176 ELSE
177 ! Copy of arrays is not covered here
178 cpabort("Copy of arrays is not covered here")
179 END IF
180
181 END SUBROUTINE contraction_ab
182
183! **************************************************************************************************
184!> \brief Applying the contraction coefficients to a tripple set integrals
185!> QABC <- CA(T) * SABC * CB * CC
186!> If only one or two of the transformation matrices are given, only a
187!> part transformation is done
188!> \param sabc Input matrix, dimension(:,:)
189!> \param qabc Output matrix, dimension(:,:)
190!> \param ca Transformation matrix (index 1), optional
191!> \param na First dimension of ca, optional
192!> \param ma Second dimension of ca, optional
193!> \param cb Transformation matrix (index 2), optional
194!> \param nb First dimension of cb, optional
195!> \param mb Second dimension of cb, optional
196!> \param cc Transformation matrix (index 3), optional
197!> \param nc First dimension of cc, optional
198!> \param mc Second dimension of cc, optional
199! **************************************************************************************************
200 SUBROUTINE contraction_abc(sabc, qabc, ca, na, ma, cb, nb, mb, cc, nc, mc)
201
202 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: sabc
203 REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: qabc
204 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
205 OPTIONAL :: ca
206 INTEGER, INTENT(IN), OPTIONAL :: na, ma
207 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
208 OPTIONAL :: cb
209 INTEGER, INTENT(IN), OPTIONAL :: nb, mb
210 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
211 OPTIONAL :: cc
212 INTEGER, INTENT(IN), OPTIONAL :: nc, mc
213
214 INTEGER :: lda, ldb, ldc, mal, mbl, mcl, nal, nbl, &
215 ncl
216 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work1, work2, work3, work4
217
218! Active matrix size
219
220 IF (PRESENT(ca)) THEN
221 IF (PRESENT(na)) THEN
222 nal = na
223 ELSE
224 nal = SIZE(ca, 1)
225 END IF
226 IF (PRESENT(ma)) THEN
227 mal = ma
228 ELSE
229 mal = SIZE(ca, 2)
230 END IF
231 lda = SIZE(ca, 1)
232 END IF
233 IF (PRESENT(cb)) THEN
234 IF (PRESENT(nb)) THEN
235 nbl = nb
236 ELSE
237 nbl = SIZE(cb, 1)
238 END IF
239 IF (PRESENT(mb)) THEN
240 mbl = mb
241 ELSE
242 mbl = SIZE(cb, 2)
243 END IF
244 ldb = SIZE(cb, 1)
245 END IF
246 IF (PRESENT(cc)) THEN
247 IF (PRESENT(nc)) THEN
248 ncl = nc
249 ELSE
250 ncl = SIZE(cc, 1)
251 END IF
252 IF (PRESENT(mc)) THEN
253 mcl = mc
254 ELSE
255 mcl = SIZE(cc, 2)
256 END IF
257 ldc = SIZE(cc, 1)
258 END IF
259
260 IF (PRESENT(ca) .AND. PRESENT(cb) .AND. PRESENT(cc)) THEN
261 ! Full transform
262 ALLOCATE (work1(nal, nbl, ncl))
263 ! make sure that we have contiguous memory, needed for transpose algorithm
264 work1(1:nal, 1:nbl, 1:ncl) = sabc(1:nal, 1:nbl, 1:ncl)
265 !
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)
269 !
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)
273 !
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)
277 !
278 qabc(1:mal, 1:mbl, 1:mcl) = work4(1:mal, 1:mbl, 1:mcl)
279 !
280 DEALLOCATE (work1, work2, work3, work4)
281 !
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")
294 ELSE
295 ! Copy of arrays is not covered here
296 cpabort("Copy of arrays is not covered here")
297 END IF
298
299 END SUBROUTINE contraction_abc
300
301! **************************************************************************************************
302!> \brief Applying the de-contraction coefficients to a matrix
303!> QAB <- CA * SAB * CB(T)
304!> Variable "trans" requests the input matrix to be SAB(T)
305!> Active dimensions are: QAB(na,nb), SAB(ma,mb)
306!> \param sab Input matrix, dimension(:,:)
307!> \param qab Output matrix, dimension(:,:)
308!> \param ca Left transformation matrix
309!> \param na First dimension of ca
310!> \param ma Second dimension of ca
311!> \param cb Right transformation matrix
312!> \param nb First dimension of cb
313!> \param mb Second dimension of cb
314!> \param trans Optional transposition of input matrix
315! **************************************************************************************************
316 SUBROUTINE decontraction_ab(sab, qab, ca, na, ma, cb, nb, mb, trans)
317
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
325
326 INTEGER :: lda, ldb, ldq, lds, ldw
327 LOGICAL :: my_trans
328 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
329
330 ! Should input matrix be transposed?
331 IF (PRESENT(trans)) THEN
332 my_trans = trans
333 ELSE
334 my_trans = .false.
335 END IF
336
337 lds = SIZE(sab, 1)
338 ldq = SIZE(qab, 1)
339 lda = SIZE(ca, 1)
340 ldb = SIZE(cb, 1)
341
342 ALLOCATE (work(na, mb))
343 ldw = na
344
345 IF (my_trans) THEN
346!dg CALL dgemm("N", "T", na, mb, ma, 1.0_dp, ca, lda, sab, lds, 0.0_dp, work, ldw)
347 work(1:na, 1:mb) = matmul(ca(1:na, 1:ma), transpose(sab(1:mb, 1:ma)))
348 ELSE
349!dg CALL dgemm("N", "N", na, mb, ma, 1.0_dp, ca, lda, sab, lds, 0.0_dp, work, ldw)
350 work(1:na, 1:mb) = matmul(ca(1:na, 1:ma), sab(1:ma, 1:mb))
351 END IF
352!dg CALL dgemm("N", "T", na, nb, mb, 1.0_dp, work, ldw, cb, ldb, 0.0_dp, qab, ldq)
353 qab(1:na, 1:nb) = matmul(work(1:na, 1:mb), transpose(cb(1:nb, 1:mb)))
354
355 DEALLOCATE (work)
356
357 END SUBROUTINE decontraction_ab
358
359! **************************************************************************************************
360!> \brief Routine to trace a series of matrices with another matrix
361!> Calculate forces of type f(:) = Trace(Pab*Sab(:))
362!> \param force Vector to hold output forces
363!> \param sab Input vector of matrices, dimension (:,:,:)
364!> \param pab Input matrix
365!> \param na Active first dimension
366!> \param nb Active second dimension
367!> \param m Number of matrices to be traced
368!> \param trans Matrices are transposed (Sab and Pab)
369! **************************************************************************************************
370 SUBROUTINE force_trace_ab(force, sab, pab, na, nb, m, trans)
371
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
377
378 INTEGER :: i
379 LOGICAL :: my_trans
380
381 cpassert(m <= SIZE(sab, 3))
382 cpassert(m <= SIZE(force, 1))
383
384 ! are matrices transposed?
385 IF (PRESENT(trans)) THEN
386 my_trans = trans
387 ELSE
388 my_trans = .false.
389 END IF
390
391 DO i = 1, m
392 IF (my_trans) THEN
393 force(i) = sum(sab(1:nb, 1:na, i)*pab(1:nb, 1:na))
394 ELSE
395 force(i) = sum(sab(1:na, 1:nb, i)*pab(1:na, 1:nb))
396 END IF
397 END DO
398
399 END SUBROUTINE force_trace_ab
400
401! **************************************************************************************************
402!> \brief Copy a block out of a matrix and add it to another matrix
403!> SAB = SAB + QAB or QAB = QAB + SAB
404!> QAB(ia:,ib:) and SAB(1:,1:)
405!> \param dir "IN" and "OUT" defines direction of copy
406!> \param sab Matrix input for "IN", output for "OUT"
407!> \param na first dimension of matrix to copy
408!> \param nb second dimension of matrix to copy
409!> \param qab Matrix output for "IN", input for "OUT"
410!> Use subblock of this matrix
411!> \param ia Starting index in qab first dimension
412!> \param ib Starting index in qab second dimension
413!> \param trans Matrices (qab and sab) are transposed
414! **************************************************************************************************
415 SUBROUTINE block_add_ab(dir, sab, na, nb, qab, ia, ib, trans)
416
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
423
424 INTEGER :: ja, jb
425 LOGICAL :: my_trans
426
427 IF (PRESENT(trans)) THEN
428 my_trans = trans
429 ELSE
430 my_trans = .false.
431 END IF
432
433 IF (dir == "IN" .OR. dir == "in") THEN
434 ! QAB(block) <= SAB
435 ja = ia + na - 1
436 jb = ib + nb - 1
437 IF (my_trans) THEN
438 qab(ib:jb, ia:ja) = qab(ib:jb, ia:ja) + sab(1:nb, 1:na)
439 ELSE
440 qab(ia:ja, ib:jb) = qab(ia:ja, ib:jb) + sab(1:na, 1:nb)
441 END IF
442 ELSEIF (dir == "OUT" .OR. dir == "out") THEN
443 ! SAB <= QAB(block)
444 ja = ia + na - 1
445 jb = ib + nb - 1
446 IF (my_trans) THEN
447 sab(1:nb, 1:na) = sab(1:nb, 1:na) + qab(ib:jb, ia:ja)
448 ELSE
449 sab(1:na, 1:nb) = sab(1:na, 1:nb) + qab(ia:ja, ib:jb)
450 END IF
451 ELSE
452 cpabort("")
453 END IF
454
455 END SUBROUTINE block_add_ab
456! **************************************************************************************************
457
458END MODULE ai_contraction
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:214
ai_contraction::block_add
Definition ai_contraction.F:43
ai_contraction::contraction
Definition ai_contraction.F:31
ai_contraction::decontraction
Definition ai_contraction.F:35
ai_contraction::force_trace
Definition ai_contraction.F:39
ai_contraction
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Definition ai_contraction.F:18
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34