Line data Source code
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 : ! **************************************************************************************************
18 : MODULE 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 :
49 : CONTAINS
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 21577442 : 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 21577442 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
89 :
90 : ! Should output matrix be transposed?
91 :
92 21577442 : IF (PRESENT(trans)) THEN
93 21576632 : my_trans = trans
94 : ELSE
95 : my_trans = .FALSE.
96 : END IF
97 :
98 : ! Scaling of output matrix
99 21577442 : IF (PRESENT(fscale)) THEN
100 21575606 : fs = fscale
101 : ELSE
102 : fs = 1.0_dp
103 : END IF
104 :
105 : ! Active matrix size
106 21577442 : IF (PRESENT(ca)) THEN
107 21577442 : IF (PRESENT(na)) THEN
108 21577442 : nal = na
109 : ELSE
110 0 : nal = SIZE(ca, 1)
111 : END IF
112 21577442 : IF (PRESENT(ma)) THEN
113 21577442 : mal = ma
114 : ELSE
115 0 : mal = SIZE(ca, 2)
116 : END IF
117 21577442 : lda = SIZE(ca, 1)
118 : END IF
119 21577442 : IF (PRESENT(cb)) THEN
120 21577442 : IF (PRESENT(nb)) THEN
121 21577442 : nbl = nb
122 : ELSE
123 0 : nbl = SIZE(cb, 1)
124 : END IF
125 21577442 : IF (PRESENT(mb)) THEN
126 21577442 : mbl = mb
127 : ELSE
128 0 : mbl = SIZE(cb, 2)
129 : END IF
130 21577442 : ldb = SIZE(cb, 1)
131 : END IF
132 :
133 21577442 : lds = SIZE(sab, 1)
134 21577442 : ldq = SIZE(qab, 1)
135 :
136 21577442 : IF (PRESENT(ca) .AND. PRESENT(cb)) THEN
137 : ! Full transform
138 86309768 : ALLOCATE (work(nal, mbl))
139 21577442 : 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 20447688344 : work(1:nal, 1:mbl) = MATMUL(sab(1:nal, 1:nbl), cb(1:nbl, 1:mbl))
142 21577442 : 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 2274620424 : 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 5651400763 : qab(1:mal, 1:mbl) = fs*MATMUL(TRANSPOSE(ca(1:nal, 1:mal)), work(1:nal, 1:mbl))
148 : END IF
149 21577442 : DEALLOCATE (work)
150 0 : ELSE IF (PRESENT(ca)) THEN
151 0 : IF (PRESENT(nb)) THEN
152 0 : nbl = nb
153 : ELSE
154 0 : nbl = SIZE(sab, 2)
155 : END IF
156 0 : 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 0 : 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 0 : qab(1:mal, 1:nbl) = fs*MATMUL(TRANSPOSE(ca(1:nal, 1:mal)), sab(1:nal, 1:nbl))
162 : END IF
163 0 : ELSE IF (PRESENT(cb)) THEN
164 0 : IF (PRESENT(na)) THEN
165 0 : nal = na
166 : ELSE
167 0 : nal = SIZE(sab, 1)
168 : END IF
169 0 : 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 0 : 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 0 : 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 0 : CPABORT("Copy of arrays is not covered here")
179 : END IF
180 :
181 21577442 : 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 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: work1, work2, work3, work4
217 :
218 : ! Active matrix size
219 :
220 0 : IF (PRESENT(ca)) THEN
221 0 : IF (PRESENT(na)) THEN
222 0 : nal = na
223 : ELSE
224 0 : nal = SIZE(ca, 1)
225 : END IF
226 0 : IF (PRESENT(ma)) THEN
227 0 : mal = ma
228 : ELSE
229 0 : mal = SIZE(ca, 2)
230 : END IF
231 0 : lda = SIZE(ca, 1)
232 : END IF
233 0 : IF (PRESENT(cb)) THEN
234 0 : IF (PRESENT(nb)) THEN
235 0 : nbl = nb
236 : ELSE
237 0 : nbl = SIZE(cb, 1)
238 : END IF
239 0 : IF (PRESENT(mb)) THEN
240 0 : mbl = mb
241 : ELSE
242 0 : mbl = SIZE(cb, 2)
243 : END IF
244 0 : ldb = SIZE(cb, 1)
245 : END IF
246 0 : IF (PRESENT(cc)) THEN
247 0 : IF (PRESENT(nc)) THEN
248 0 : ncl = nc
249 : ELSE
250 0 : ncl = SIZE(cc, 1)
251 : END IF
252 0 : IF (PRESENT(mc)) THEN
253 0 : mcl = mc
254 : ELSE
255 0 : mcl = SIZE(cc, 2)
256 : END IF
257 0 : ldc = SIZE(cc, 1)
258 : END IF
259 :
260 0 : IF (PRESENT(ca) .AND. PRESENT(cb) .AND. PRESENT(cc)) THEN
261 : ! Full transform
262 0 : ALLOCATE (work1(nal, nbl, ncl))
263 : ! make sure that we have contiguous memory, needed for transpose algorithm
264 0 : work1(1:nal, 1:nbl, 1:ncl) = sabc(1:nal, 1:nbl, 1:ncl)
265 : !
266 0 : 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.0_dp, work2(1, 1, 1), nbl*ncl)
269 : !
270 0 : 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.0_dp, work3(1, 1, 1), ncl*mal)
273 : !
274 0 : 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.0_dp, work4(1, 1, 1), mal*mbl)
277 : !
278 0 : qabc(1:mal, 1:mbl, 1:mcl) = work4(1:mal, 1:mbl, 1:mcl)
279 : !
280 0 : DEALLOCATE (work1, work2, work3, work4)
281 : !
282 0 : ELSE IF (PRESENT(ca) .AND. PRESENT(cb)) THEN
283 0 : CPABORT("Not implemented")
284 0 : ELSE IF (PRESENT(ca) .AND. PRESENT(cc)) THEN
285 0 : CPABORT("Not implemented")
286 0 : ELSE IF (PRESENT(cb) .AND. PRESENT(cc)) THEN
287 0 : CPABORT("Not implemented")
288 0 : ELSE IF (PRESENT(ca)) THEN
289 0 : CPABORT("Not implemented")
290 0 : ELSE IF (PRESENT(cb)) THEN
291 0 : CPABORT("Not implemented")
292 0 : ELSE IF (PRESENT(cc)) THEN
293 0 : CPABORT("Not implemented")
294 : ELSE
295 : ! Copy of arrays is not covered here
296 0 : CPABORT("Copy of arrays is not covered here")
297 : END IF
298 :
299 0 : 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 1704713 : 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 1704713 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
329 :
330 : ! Should input matrix be transposed?
331 1704713 : IF (PRESENT(trans)) THEN
332 1704713 : my_trans = trans
333 : ELSE
334 : my_trans = .FALSE.
335 : END IF
336 :
337 1704713 : lds = SIZE(sab, 1)
338 1704713 : ldq = SIZE(qab, 1)
339 1704713 : lda = SIZE(ca, 1)
340 1704713 : ldb = SIZE(cb, 1)
341 :
342 6818852 : ALLOCATE (work(na, mb))
343 1704713 : ldw = na
344 :
345 1704713 : 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 655918707 : 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 605641340 : 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 3638728355 : qab(1:na, 1:nb) = MATMUL(work(1:na, 1:mb), TRANSPOSE(cb(1:nb, 1:mb)))
354 :
355 1704713 : DEALLOCATE (work)
356 :
357 1704713 : 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 1704713 : 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 1704713 : CPASSERT(m <= SIZE(SAB, 3))
382 1704713 : CPASSERT(m <= SIZE(force, 1))
383 :
384 : ! are matrices transposed?
385 1704713 : IF (PRESENT(trans)) THEN
386 0 : my_trans = trans
387 : ELSE
388 : my_trans = .FALSE.
389 : END IF
390 :
391 6818852 : DO i = 1, m
392 6818852 : IF (my_trans) THEN
393 0 : force(i) = SUM(sab(1:nb, 1:na, i)*pab(1:nb, 1:na))
394 : ELSE
395 1094926710 : force(i) = SUM(sab(1:na, 1:nb, i)*pab(1:na, 1:nb))
396 : END IF
397 : END DO
398 :
399 1704713 : 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 24050706 : 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 24050706 : IF (PRESENT(trans)) THEN
428 24049896 : my_trans = trans
429 : ELSE
430 : my_trans = .FALSE.
431 : END IF
432 :
433 24050706 : IF (dir == "IN" .OR. dir == "in") THEN
434 : ! QAB(block) <= SAB
435 22345993 : ja = ia + na - 1
436 22345993 : jb = ib + nb - 1
437 22345993 : IF (my_trans) THEN
438 193087649 : qab(ib:jb, ia:ja) = qab(ib:jb, ia:ja) + sab(1:nb, 1:na)
439 : ELSE
440 411360533 : qab(ia:ja, ib:jb) = qab(ia:ja, ib:jb) + sab(1:na, 1:nb)
441 : END IF
442 1704713 : ELSEIF (dir == "OUT" .OR. dir == "out") THEN
443 : ! SAB <= QAB(block)
444 1704713 : ja = ia + na - 1
445 1704713 : jb = ib + nb - 1
446 1704713 : IF (my_trans) THEN
447 31895236 : sab(1:nb, 1:na) = sab(1:nb, 1:na) + qab(ib:jb, ia:ja)
448 : ELSE
449 27153533 : sab(1:na, 1:nb) = sab(1:na, 1:nb) + qab(ia:ja, ib:jb)
450 : END IF
451 : ELSE
452 0 : CPABORT("")
453 : END IF
454 :
455 24050706 : END SUBROUTINE block_add_ab
456 : ! **************************************************************************************************
457 :
458 : END MODULE ai_contraction
|
