Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 : ! **************************************************************************************************
8 : !> \brief routines to contract density matrix blocks with the four-center
9 : !> integrals to yield the Kohn-Sham matrix. The specialized routines
10 : !> are about 1.2-2.0 as fast as the default one.
11 : !> \par History
12 : !> 10.2009 created [Joost VandeVondele]
13 : !> \author Joost VandeVondele
14 : ! **************************************************************************************************
15 : MODULE hfx_contract_block
16 : USE kinds, ONLY: dp
17 : #include "../base/base_uses.f90"
18 :
19 : IMPLICIT NONE
20 : PRIVATE
21 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hfx_contract_block'
22 : PUBLIC :: contract_block
23 : CONTAINS
24 :
25 : #! Keep the original dimensions and limit for further specialization.
26 : #:set block_sizes = [1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 15]
27 : #:set max_product = 20
28 : #:set block_args = ['kbd', 'kbc', 'kad', 'kac', 'pbd', 'pbc', 'pad', 'pac', 'prim', 'scale']
29 :
30 : #:def call_block(fixed=())
31 : #:set name = 'block_' + ('_'.join(map(str, fixed)) if fixed else 'default')
32 : #:set runtime_dims = ['m' + axis + '_max' for axis in 'abcd'[len(fixed):]]
33 : CALL ${name}$ (${', '.join(runtime_dims + block_args)}$)
34 : #:enddef
35 :
36 : #! The same loop body generates the default kernel and every specialization.
37 : #! Fixed dimensions remain literal constants in the generated Fortran code.
38 : #:def block(fixed=())
39 : #:set name = 'block_' + ('_'.join(map(str, fixed)) if fixed else 'default')
40 : #:set runtime_dims = ['m' + axis + '_max' for axis in 'abcd'[len(fixed):]]
41 : #:set a, b, c, d = list(map(str, fixed)) + runtime_dims
42 : ! **************************************************************************************************
43 : !> \brief Contract density matrix blocks with four-center integrals.
44 : #:for dim in runtime_dims
45 : !> \param ${dim}$ Block dimension
46 : #:endfor
47 : #:for arg in block_args
48 : !> \param ${arg}$ ...
49 : #:endfor
50 : ! **************************************************************************************************
51 : SUBROUTINE ${name}$ (${', '.join(runtime_dims + block_args)}$)
52 : #:if runtime_dims
53 : INTEGER :: ${', '.join(runtime_dims)}$
54 : #:endif
55 : REAL(KIND=dp) :: kbd(${b}$*${d}$), kbc(${b}$*${c}$), kad(${a}$*${d}$), &
56 : kac(${a}$*${c}$), pbd(${b}$*${d}$), pbc(${b}$*${c}$), pad(${a}$*${d}$), &
57 : pac(${a}$*${c}$), prim(${a}$*${b}$*${c}$*${d}$), scale
58 :
59 : INTEGER :: ma, mb, mc, md, p_index
60 : REAL(KIND=dp) :: ks_bc, ks_bd, p_bc, p_bd, tmp
61 :
62 : kbd(1:${b}$*${d}$) = 0.0_dp
63 : kbc(1:${b}$*${c}$) = 0.0_dp
64 : kad(1:${a}$*${d}$) = 0.0_dp
65 : kac(1:${a}$*${c}$) = 0.0_dp
66 : p_index = 0
67 : DO md = 1, ${d}$
68 : DO mc = 1, ${c}$
69 : DO mb = 1, ${b}$
70 : ks_bd = 0.0_dp
71 : ks_bc = 0.0_dp
72 : p_bd = pbd((md - 1)*${b}$+mb)
73 : p_bc = pbc((mc - 1)*${b}$+mb)
74 : DO ma = 1, ${a}$
75 : p_index = p_index + 1
76 : tmp = scale*prim(p_index)
77 : ks_bc = ks_bc + tmp*pad((md - 1)*${a}$+ma)
78 : ks_bd = ks_bd + tmp*pac((mc - 1)*${a}$+ma)
79 : kad((md - 1)*${a}$+ma) = kad((md - 1)*${a}$+ma) - tmp*p_bc
80 : kac((mc - 1)*${a}$+ma) = kac((mc - 1)*${a}$+ma) - tmp*p_bd
81 : END DO
82 : kbd((md - 1)*${b}$+mb) = kbd((md - 1)*${b}$+mb) - ks_bd
83 : kbc((mc - 1)*${b}$+mb) = kbc((mc - 1)*${b}$+mb) - ks_bc
84 : END DO
85 : END DO
86 : END DO
87 : END SUBROUTINE ${name}$
88 : #:enddef
89 :
90 : ! **************************************************************************************************
91 : !> \brief ...
92 : !> \param ma_max ...
93 : !> \param mb_max ...
94 : !> \param mc_max ...
95 : !> \param md_max ...
96 : !> \param kbd ...
97 : !> \param kbc ...
98 : !> \param kad ...
99 : !> \param kac ...
100 : !> \param pbd ...
101 : !> \param pbc ...
102 : !> \param pad ...
103 : !> \param pac ...
104 : !> \param prim ...
105 : !> \param scale ...
106 : ! **************************************************************************************************
107 87104725 : SUBROUTINE contract_block(ma_max, mb_max, mc_max, md_max, kbd, kbc, kad, kac, pbd, pbc, pad, pac, prim, scale)
108 : INTEGER :: ma_max, mb_max, mc_max, md_max
109 : REAL(KIND=dp) :: kbd(mb_max*md_max), kbc(mb_max*mc_max), &
110 : kad(ma_max*md_max), kac(ma_max*mc_max), pbd(mb_max*md_max), &
111 : pbc(mb_max*mc_max), pad(ma_max*md_max), pac(ma_max*mc_max), &
112 : prim(ma_max*mb_max*mc_max*md_max), scale
113 :
114 : #if !defined (__LIBINT)
115 : MARK_USED(ma_max)
116 : MARK_USED(mb_max)
117 : MARK_USED(mc_max)
118 : MARK_USED(md_max)
119 : MARK_USED(kbd)
120 : MARK_USED(kbc)
121 : MARK_USED(kad)
122 : MARK_USED(kac)
123 : MARK_USED(pbd)
124 : MARK_USED(pbc)
125 : MARK_USED(pad)
126 : MARK_USED(pac)
127 : MARK_USED(prim)
128 : MARK_USED(scale)
129 : CPABORT("libint not compiled in")
130 : #else
131 49895169 : SELECT CASE (ma_max)
132 : #:for a in block_sizes
133 : CASE (${a}$)
134 99210415 : SELECT CASE (mb_max)
135 : #:for b in block_sizes
136 : CASE (${b}$)
137 : #:if a*b < max_product
138 65211629 : SELECT CASE (mc_max)
139 : #:for c in block_sizes
140 : CASE (${c}$)
141 : #:if a*b*c < max_product
142 86844475 : SELECT CASE (md_max)
143 : #:for d in block_sizes
144 : CASE (${d}$)
145 : #:if a*b*c*d < max_product
146 65222684 : $:call_block((a, b, c, d))
147 : #:else
148 11198296 : $:call_block((a, b, c))
149 : #:endif
150 : #:endfor
151 : CASE DEFAULT
152 76420980 : $:call_block()
153 : END SELECT
154 : #:else
155 8592190 : $:call_block((a, b))
156 : #:endif
157 : #:endfor
158 : CASE DEFAULT
159 85287896 : $:call_block()
160 : END SELECT
161 : #:else
162 1471838 : $:call_block((a, b))
163 : #:endif
164 : #:endfor
165 : CASE DEFAULT
166 86807710 : $:call_block()
167 : END SELECT
168 : #:endfor
169 : CASE DEFAULT
170 87104725 : $:call_block()
171 : END SELECT
172 : #endif
173 87104725 : END SUBROUTINE contract_block
174 :
175 : #if defined(__LIBINT)
176 661090672 : $:block()
177 :
178 : #:for a in block_sizes
179 : #:for b in block_sizes
180 : #:if a*b < max_product
181 : #:for c in block_sizes
182 : #:if a*b*c < max_product
183 : #:for d in block_sizes
184 : #:if a*b*c*d < max_product
185 952534444 : $:block((a, b, c, d))
186 : #:endif
187 : #:endfor
188 : #:if a*b*c*block_sizes[-1] >= max_product
189 950608353 : $:block((a, b, c))
190 : #:endif
191 : #:endif
192 : #:endfor
193 : #:endif
194 : #:if a*b*block_sizes[-1] >= max_product
195 2277945521 : $:block((a, b))
196 : #:endif
197 : #:endfor
198 : #:endfor
199 : #endif
200 : END MODULE hfx_contract_block
|
