(git:58e3e09)
ai_contraction.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ! **************************************************************************************************
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  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 
458 END MODULE ai_contraction
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.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34