Line data Source code
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: BSD-3-Clause */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : #include "grid_dgemm_prepare_pab.h"
9 :
10 : #include <assert.h>
11 : #include <stdbool.h>
12 :
13 : #include "../common/grid_common.h"
14 : #include "../common/grid_constants.h"
15 : #include "grid_dgemm_utils.h"
16 :
17 : struct pab_computation_struct_ {
18 : int offset[2];
19 : int lmax[2];
20 : int lmin[2];
21 : double zeta[2];
22 : tensor *pab;
23 : tensor *pab_prep;
24 : int dir1, dir2;
25 : };
26 :
27 : // *****************************************************************************
28 44389 : static void grid_prepare_pab_AB(struct pab_computation_struct_ *const tp) {
29 117327 : for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
30 195459 : for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
31 300574 : for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
32 437913 : for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
33 259860 : for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
34 562143 : lza <= tp->lmax[0] - lxa - lya; lza++) {
35 302283 : for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
36 655153 : lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
37 352870 : const int ico = tp->offset[0] + coset(lxa, lya, lza);
38 352870 : const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
39 352870 : idx2(tp->pab_prep[0], coset(lxb, lyb, lzb),
40 352870 : coset(lxa, lya, lza)) = idx2(tp->pab[0], jco, ico);
41 : }
42 : }
43 : }
44 : }
45 : }
46 : }
47 44389 : }
48 :
49 : // *****************************************************************************
50 0 : static void grid_prepare_pab_DADB(struct pab_computation_struct_ *const tp) {
51 : // create a new pab_prep so that mapping pab_prep with pgf_a pgf_b
52 : // is equivalent to mapping pab with 0.5 * (nabla pgf_a) . (nabla pgf_b)
53 : // (ddx pgf_a ) (ddx pgf_b) = (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x})*(lbx
54 : // pgf_{b-1x} - 2*tp->zeta[1]*pgf_{b+1x})
55 :
56 0 : for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
57 0 : for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
58 0 : for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
59 0 : for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
60 0 : for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
61 0 : lza <= tp->lmax[0] - lxa - lya; lza++) {
62 0 : for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
63 0 : lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
64 0 : const int ico = tp->offset[0] + coset(lxa, lya, lza);
65 0 : const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
66 0 : const double pab = idx2(tp->pab[0], jco, ico);
67 0 : int ico_l, jco_l;
68 : // x (all safe if lxa = 0, as the spurious added terms have zero
69 : // prefactor)
70 :
71 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
72 0 : jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
73 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += 0.5 * lxa * lxb * pab;
74 :
75 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
76 0 : jco_l = coset((lxb + 1), lyb, lzb);
77 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * tp->zeta[1] * pab;
78 :
79 0 : ico_l = coset((lxa + 1), lya, lza);
80 0 : jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
81 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= tp->zeta[0] * lxb * pab;
82 :
83 0 : ico_l = coset((lxa + 1), lya, lza);
84 0 : jco_l = coset((lxb + 1), lyb, lzb);
85 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
86 0 : 2.0 * tp->zeta[0] * tp->zeta[1] * pab;
87 :
88 : // y
89 :
90 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
91 0 : jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
92 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += 0.5 * lya * lyb * pab;
93 :
94 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
95 0 : jco_l = coset(lxb, (lyb + 1), lzb);
96 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * tp->zeta[1] * pab;
97 :
98 0 : ico_l = coset(lxa, (lya + 1), lza);
99 0 : jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
100 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= tp->zeta[0] * lyb * pab;
101 :
102 0 : ico_l = coset(lxa, (lya + 1), lza);
103 0 : jco_l = coset(lxb, (lyb + 1), lzb);
104 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
105 0 : 2.0 * tp->zeta[0] * tp->zeta[1] * pab;
106 :
107 : // z
108 :
109 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
110 0 : jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
111 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += 0.5 * lza * lzb * pab;
112 :
113 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
114 0 : jco_l = coset(lxb, lyb, (lzb + 1));
115 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * tp->zeta[1] * pab;
116 :
117 0 : ico_l = coset(lxa, lya, (lza + 1));
118 0 : jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
119 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= tp->zeta[0] * lzb * pab;
120 :
121 0 : ico_l = coset(lxa, lya, (lza + 1));
122 0 : jco_l = coset(lxb, lyb, (lzb + 1));
123 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
124 0 : 2.0 * tp->zeta[0] * tp->zeta[1] * pab;
125 : }
126 : }
127 : }
128 : }
129 : }
130 : }
131 0 : }
132 :
133 : // *****************************************************************************
134 0 : static void grid_prepare_pab_ADBmDAB(struct pab_computation_struct_ *const tp) {
135 : // create a new pab_prep so that mapping pab_prep with pgf_a pgf_b
136 : // is equivalent to mapping pab with
137 : // pgf_a (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) pgf_b
138 : // ( pgf_a ) (ddx pgf_b) - (ddx pgf_a)( pgf_b ) =
139 : // pgf_a *(lbx pgf_{b-1x} - 2*tp->zeta[1]*pgf_{b+1x}) -
140 : // (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_b
141 :
142 0 : for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
143 0 : for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
144 0 : for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
145 0 : for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
146 0 : for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
147 0 : lza <= tp->lmax[0] - lxa - lya; lza++) {
148 0 : for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
149 0 : lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
150 0 : const int ico = tp->offset[0] + coset(lxa, lya, lza);
151 0 : const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
152 0 : const double pab = idx2(tp->pab[0], jco, ico);
153 0 : int ico_l, jco_l;
154 :
155 : // ! this element of pab results in 4 elements of pab_prep
156 0 : switch (tp->dir1) {
157 0 : case 'X': { // x
158 0 : ico_l = coset(lxa, lya, lza);
159 0 : jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
160 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
161 :
162 : // ico_l = coset(lxa, lya, lza);
163 0 : jco_l = coset((lxb + 1), lyb, lzb);
164 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
165 :
166 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
167 0 : jco_l = coset(lxb, lyb, lzb);
168 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * pab;
169 :
170 0 : ico_l = coset(lxa + 1, lya, lza);
171 : // jco_l = coset(lxb, lyb, lzb);
172 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += 2.0 * tp->zeta[0] * pab;
173 0 : } break;
174 0 : case 'Y': { // y
175 0 : ico_l = coset(lxa, lya, lza);
176 0 : jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
177 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
178 :
179 : // ico_l = coset(lxa, lya, lza);
180 0 : jco_l = coset(lxb, lyb + 1, lzb);
181 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
182 :
183 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
184 0 : jco_l = coset(lxb, lyb, lzb);
185 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * pab;
186 :
187 0 : ico_l = coset(lxa, lya + 1, lza);
188 : // jco_l = coset(lxb, lyb, lzb);
189 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += 2.0 * tp->zeta[0] * pab;
190 0 : } break;
191 0 : case 'Z': { // z
192 0 : ico_l = coset(lxa, lya, lza);
193 0 : jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
194 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
195 :
196 : // ico_l = coset(lxa, lya, lza);
197 0 : jco_l = coset(lxb, lyb, lzb + 1);
198 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
199 :
200 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
201 0 : jco_l = coset(lxb, lyb, lzb);
202 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * pab;
203 :
204 0 : ico_l = coset(lxa, lya, lza + 1);
205 : // jco_l = coset(lxb, lyb, lzb);
206 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += 2.0 * tp->zeta[0] * pab;
207 0 : } break;
208 : default:
209 : break;
210 : }
211 : }
212 : }
213 : }
214 : }
215 : }
216 : }
217 0 : }
218 : // *****************************************************************************
219 : static void
220 0 : grid_prepare_pab_ARDBmDARB(struct pab_computation_struct_ *const tp) {
221 : // create a new pab_prep so that mapping pab_prep with pgf_a pgf_b
222 : // is equivalent to mapping pab with
223 : // pgf_a (r-Rb)_{ir} (nabla_{idir} pgf_b) - (nabla_{idir} pgf_a) (r-Rb)_{ir}
224 : // pgf_b ( pgf_a )(r-Rb)_{ir} (ddx pgf_b) - (ddx pgf_a) (r-Rb)_{ir} ( pgf_b )
225 : // =
226 : // pgf_a *(lbx pgf_{b-1x+1ir} -
227 : // 2*tp->zeta[1]*pgf_{b+1x+1ir}) -
228 : // (lax pgf_{a-1x} - 2*zeta*pgf_{a+1x}) pgf_{b+1ir}
229 :
230 0 : assert(1 <= tp->dir1 && tp->dir1 <= 3);
231 0 : assert(1 <= tp->dir2 && tp->dir2 <= 3);
232 :
233 0 : for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
234 0 : for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
235 0 : for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
236 0 : for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
237 0 : for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
238 0 : lza <= tp->lmax[0] - lxa - lya; lza++) {
239 0 : for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
240 0 : lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
241 0 : const int ico = tp->offset[0] + coset(lxa, lya, lza);
242 0 : const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
243 0 : const double pab = idx2(tp->pab[0], jco, ico);
244 :
245 0 : int ico_l, jco_l;
246 :
247 : // this element of pab results in 4 elements of pab_prep
248 0 : switch (tp->dir1) {
249 0 : case 'X': {
250 0 : switch (tp->dir2) {
251 0 : case 'X': {
252 0 : ico_l = coset(lxa, lya, lza);
253 0 : jco_l = coset(lxb, lyb, lzb);
254 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
255 :
256 0 : ico_l = coset(lxa, lya, lza);
257 0 : jco_l = coset((lxb + 2), lyb, lzb);
258 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
259 0 : 2.0 * tp->zeta[1] * pab;
260 :
261 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
262 0 : jco_l = coset((lxb + 1), lyb, lzb);
263 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * pab;
264 :
265 0 : ico_l = coset((lxa + 1), lya, lza);
266 0 : jco_l = coset((lxb + 1), lyb, lzb);
267 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
268 0 : 2.0 * tp->zeta[0] * pab;
269 0 : } break;
270 0 : case 'Y': {
271 0 : ico_l = coset(lxa, lya, lza);
272 0 : jco_l = coset(imax(lxb - 1, 0), (lyb + 1), lzb);
273 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
274 :
275 0 : ico_l = coset(lxa, lya, lza);
276 0 : jco_l = coset((lxb + 1), (lyb + 1), lzb);
277 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
278 0 : 2.0 * tp->zeta[1] * pab;
279 :
280 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
281 0 : jco_l = coset(lxb, (lyb + 1), lzb);
282 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * pab;
283 :
284 0 : ico_l = coset((lxa + 1), lya, lza);
285 0 : jco_l = coset(lxb, (lyb + 1), lzb);
286 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
287 0 : 2.0 * tp->zeta[0] * pab;
288 0 : } break;
289 0 : case 'Z': {
290 0 : ico_l = coset(lxa, lya, lza);
291 0 : jco_l = coset(imax(lxb - 1, 0), lyb, (lzb + 1));
292 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
293 :
294 0 : ico_l = coset(lxa, lya, lza);
295 0 : jco_l = coset((lxb + 1), lyb, (lzb + 1));
296 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
297 0 : 2.0 * tp->zeta[1] * pab;
298 :
299 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
300 0 : jco_l = coset(lxb, lyb, (lzb + 1));
301 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lxa * pab;
302 :
303 0 : ico_l = coset((lxa + 1), lya, lza);
304 0 : jco_l = coset(lxb, lyb, (lzb + 1));
305 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
306 0 : 2.0 * tp->zeta[0] * pab;
307 0 : } break;
308 : default:
309 : break;
310 : }
311 : } break;
312 0 : case 'Y': {
313 0 : switch (tp->dir2) {
314 0 : case 'X': {
315 0 : ico_l = coset(lxa, lya, lza);
316 0 : jco_l = coset((lxb + 1), imax(lyb - 1, 0), lzb);
317 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
318 :
319 0 : ico_l = coset(lxa, lya, lza);
320 0 : jco_l = coset((lxb + 1), (lyb + 1), lzb);
321 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
322 0 : 2.0 * tp->zeta[1] * pab;
323 :
324 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
325 0 : jco_l = coset((lxb + 1), lyb, lzb);
326 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * pab;
327 :
328 0 : ico_l = coset(lxa, (lya + 1), lza);
329 0 : jco_l = coset((lxb + 1), lyb, lzb);
330 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
331 0 : 2.0 * tp->zeta[0] * pab;
332 0 : } break;
333 0 : case 'Y': {
334 0 : ico_l = coset(lxa, lya, lza);
335 0 : jco_l = coset(lxb, lyb, lzb);
336 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
337 :
338 0 : ico_l = coset(lxa, lya, lza);
339 0 : jco_l = coset(lxb, (lyb + 2), lzb);
340 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
341 0 : 2.0 * tp->zeta[1] * pab;
342 :
343 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
344 0 : jco_l = coset(lxb, (lyb + 1), lzb);
345 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * pab;
346 :
347 0 : ico_l = coset(lxa, (lya + 1), lza);
348 0 : jco_l = coset(lxb, (lyb + 1), lzb);
349 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
350 0 : 2.0 * tp->zeta[0] * pab;
351 0 : } break;
352 0 : case 'Z': {
353 0 : ico_l = coset(lxa, lya, lza);
354 0 : jco_l = coset(lxb, imax(lyb - 1, 0), (lzb + 1));
355 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
356 :
357 0 : ico_l = coset(lxa, lya, lza);
358 0 : jco_l = coset(lxb, (lyb + 1), (lzb + 1));
359 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
360 0 : 2.0 * tp->zeta[1] * pab;
361 :
362 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
363 0 : jco_l = coset(lxb, lyb, (lzb + 1));
364 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lya * pab;
365 :
366 0 : ico_l = coset(lxa, (lya + 1), lza);
367 0 : jco_l = coset(lxb, lyb, (lzb + 1));
368 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
369 0 : 2.0 * tp->zeta[0] * pab;
370 0 : } break;
371 : default:
372 : break;
373 : }
374 : } break;
375 0 : case 'Z': {
376 0 : switch (tp->dir2) {
377 0 : case 'X': {
378 0 : ico_l = coset(lxa, lya, lza);
379 0 : jco_l = coset((lxb + 1), lyb, imax(lzb - 1, 0));
380 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
381 :
382 0 : ico_l = coset(lxa, lya, lza);
383 0 : jco_l = coset((lxb + 1), lyb, (lzb + 1));
384 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
385 0 : 2.0 * tp->zeta[1] * pab;
386 :
387 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
388 0 : jco_l = coset((lxb + 1), lyb, lzb);
389 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * pab;
390 :
391 0 : ico_l = coset(lxa, lya, (lza + 1));
392 0 : jco_l = coset((lxb + 1), lyb, lzb);
393 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
394 0 : 2.0 * tp->zeta[0] * pab;
395 0 : } break;
396 0 : case 'Y': {
397 0 : ico_l = coset(lxa, lya, lza);
398 0 : jco_l = coset(lxb, (lyb + 1), imax(lzb - 1, 0));
399 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
400 :
401 0 : ico_l = coset(lxa, lya, lza);
402 0 : jco_l = coset(lxb, (lyb + 1), (lzb + 1));
403 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
404 0 : -2.0 * tp->zeta[1] * pab;
405 :
406 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
407 0 : jco_l = coset(lxb, (lyb + 1), lzb);
408 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * pab;
409 :
410 0 : ico_l = coset(lxa, lya, (lza + 1));
411 0 : jco_l = coset(lxb, (lyb + 1), lzb);
412 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
413 0 : 2.0 * tp->zeta[0] * pab;
414 0 : } break;
415 0 : case 'Z': {
416 0 : ico_l = coset(lxa, lya, lza);
417 0 : jco_l = coset(lxb, lyb, lzb);
418 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
419 :
420 0 : ico_l = coset(lxa, lya, lza);
421 0 : jco_l = coset(lxb, lyb, (lzb + 2));
422 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
423 0 : 2.0 * tp->zeta[1] * pab;
424 :
425 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
426 0 : jco_l = coset(lxb, lyb, (lzb + 1));
427 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= lza * pab;
428 :
429 0 : ico_l = coset(lxa, lya, (lza + 1));
430 0 : jco_l = coset(lxb, lyb, (lzb + 1));
431 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
432 0 : 2.0 * tp->zeta[0] * pab;
433 0 : } break;
434 : default:
435 : break;
436 : }
437 : } break;
438 : default:
439 : break;
440 : }
441 : }
442 : }
443 : }
444 : }
445 : }
446 : }
447 0 : }
448 : // *****************************************************************************
449 0 : static void grid_prepare_pab_DABpADB(struct pab_computation_struct_ *const tp) {
450 : // create a new pab_prep so that mapping pab_prep with pgf_a pgf_b
451 : // is equivalent to mapping pab with
452 : // pgf_a (nabla_{idir} pgf_b) + (nabla_{idir} pgf_a) pgf_b
453 : // ( pgf_a ) (ddx pgf_b) + (ddx pgf_a)( pgf_b ) =
454 : // pgf_a *(lbx pgf_{b-1x} + 2*tp->zeta[1]*pgf_{b+1x}) +
455 : // (lax pgf_{a-1x} + 2*zeta*pgf_{a+1x}) pgf_b
456 0 : for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
457 0 : for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
458 0 : for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
459 0 : for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
460 0 : for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
461 0 : lza <= tp->lmax[0] - lxa - lya; lza++) {
462 0 : for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
463 0 : lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
464 0 : const int ico = tp->offset[0] + coset(lxa, lya, lza);
465 0 : const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
466 0 : const double pab = idx2(tp->pab[0], jco, ico);
467 :
468 0 : int ico_l, jco_l;
469 :
470 : // this element of pab results in 4 elements of pab_prep
471 :
472 0 : switch (tp->dir1) {
473 0 : case 'X': {
474 0 : ico_l = coset(lxa, lya, lza);
475 0 : jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
476 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lxb * pab;
477 :
478 0 : ico_l = coset(lxa, lya, lza);
479 0 : jco_l = coset((lxb + 1), lyb, lzb);
480 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
481 :
482 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
483 0 : jco_l = coset(lxb, lyb, lzb);
484 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lxa * pab;
485 :
486 0 : ico_l = coset((lxa + 1), lya, lza);
487 0 : jco_l = coset(lxb, lyb, lzb);
488 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[0] * pab;
489 0 : } break;
490 0 : case 'Y': { // y
491 0 : ico_l = coset(lxa, lya, lza);
492 0 : jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
493 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lyb * pab;
494 :
495 0 : ico_l = coset(lxa, lya, lza);
496 0 : jco_l = coset(lxb, (lyb + 1), lzb);
497 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
498 :
499 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
500 0 : jco_l = coset(lxb, lyb, lzb);
501 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lya * pab;
502 :
503 0 : ico_l = coset(lxa, (lya + 1), lza);
504 0 : jco_l = coset(lxb, lyb, lzb);
505 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[0] * pab;
506 0 : } break;
507 0 : case 'Z': { // z
508 0 : ico_l = coset(lxa, lya, lza);
509 0 : jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
510 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lzb * pab;
511 :
512 0 : ico_l = coset(lxa, lya, lza);
513 0 : jco_l = coset(lxb, lyb, (lzb + 1));
514 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[1] * pab;
515 :
516 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
517 0 : jco_l = coset(lxb, lyb, lzb);
518 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lza * pab;
519 :
520 0 : ico_l = coset(lxa, lya, (lza + 1));
521 0 : jco_l = coset(lxb, lyb, lzb);
522 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -= 2.0 * tp->zeta[0] * pab;
523 0 : break;
524 : }
525 : default:
526 : break;
527 : }
528 : }
529 : }
530 : }
531 : }
532 : }
533 : }
534 0 : }
535 : // *****************************************************************************
536 0 : static void grid_prepare_pab_Di(struct pab_computation_struct_ *const tp) {
537 : // create a new pab_local so that mapping pab_local with pgf_a pgf_b
538 : // is equivalent to mapping pab with
539 : // d_{ider1} pgf_a d_{ider1} pgf_b
540 : // dx pgf_a dx pgf_b =
541 : // (lax pgf_{a-1x})*(lbx pgf_{b-1x}) -
542 : // 2*tp->zeta[1]*lax*pgf_{a-1x}*pgf_{b+1x} -
543 : // lbx pgf_{b-1x}*2*zeta*pgf_{a+1x}+ 4*zeta*zetab*pgf_{a+1x}pgf_{b+1x}
544 :
545 0 : for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
546 0 : for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
547 0 : for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
548 0 : for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
549 0 : for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
550 0 : lza <= tp->lmax[0] - lxa - lya; lza++) {
551 0 : for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
552 0 : lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
553 0 : const int ico = tp->offset[0] + coset(lxa, lya, lza);
554 0 : const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
555 0 : const double pab = idx2(tp->pab[0], jco, ico);
556 :
557 0 : int ico_l, jco_l;
558 : // this element of pab results in 12 elements of pab_prep
559 :
560 0 : switch (tp->dir1) {
561 0 : case 'X': {
562 : // x (all safe if lxa = 0, as the spurious added terms have
563 : // zero prefactor)
564 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
565 0 : jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
566 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lxa * lxb * pab;
567 :
568 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza);
569 0 : jco_l = coset((lxb + 1), lyb, lzb);
570 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
571 0 : 2.0 * lxa * tp->zeta[1] * pab;
572 :
573 0 : ico_l = coset((lxa + 1), lya, lza);
574 0 : jco_l = coset(imax(lxb - 1, 0), lyb, lzb);
575 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
576 0 : 2.0 * tp->zeta[0] * lxb * pab;
577 :
578 0 : ico_l = coset((lxa + 1), lya, lza);
579 0 : jco_l = coset((lxb + 1), lyb, lzb);
580 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
581 0 : 4.0 * tp->zeta[0] * tp->zeta[1] * pab;
582 0 : } break;
583 0 : case 'Y': {
584 : // y
585 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
586 0 : jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
587 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lya * lyb * pab;
588 :
589 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza);
590 0 : jco_l = coset(lxb, (lyb + 1), lzb);
591 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
592 0 : 2.0 * lya * tp->zeta[1] * pab;
593 :
594 0 : ico_l = coset(lxa, (lya + 1), lza);
595 0 : jco_l = coset(lxb, imax(lyb - 1, 0), lzb);
596 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
597 0 : 2.0 * tp->zeta[0] * lyb * pab;
598 :
599 0 : ico_l = coset(lxa, (lya + 1), lza);
600 0 : jco_l = coset(lxb, (lyb + 1), lzb);
601 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
602 0 : 4.0 * tp->zeta[0] * tp->zeta[1] * pab;
603 0 : } break;
604 0 : case 'Z': {
605 : // z
606 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
607 0 : jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
608 0 : idx2(tp->pab_prep[0], jco_l, ico_l) += lza * lzb * pab;
609 :
610 0 : ico_l = coset(lxa, lya, imax(lza - 1, 0));
611 0 : jco_l = coset(lxb, lyb, (lzb + 1));
612 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
613 0 : 2.0 * lza * tp->zeta[1] * pab;
614 :
615 0 : ico_l = coset(lxa, lya, (lza + 1));
616 0 : jco_l = coset(lxb, lyb, imax(lzb - 1, 0));
617 0 : idx2(tp->pab_prep[0], jco_l, ico_l) -=
618 0 : 2.0 * tp->zeta[0] * lzb * pab;
619 :
620 0 : ico_l = coset(lxa, lya, (lza + 1));
621 0 : jco_l = coset(lxb, lyb, (lzb + 1));
622 0 : idx2(tp->pab_prep[0], jco_l, ico_l) +=
623 0 : 4.0 * tp->zeta[0] * tp->zeta[1] * pab;
624 0 : } break;
625 : default:
626 : break;
627 : }
628 : }
629 : }
630 : }
631 : }
632 : }
633 : }
634 0 : }
635 :
636 : // *****************************************************************************
637 0 : static void oneterm_dijdij(const int idir, const double func_a, const int ico_l,
638 : const int lx, const int ly, const int lz,
639 : const double zet, tensor *const pab_prep) {
640 0 : int jco_l;
641 :
642 0 : switch (idir) {
643 0 : case 'X': {
644 0 : const int l1 = lx;
645 0 : const int l2 = ly;
646 0 : jco_l = coset(imax(lx - 1, 0), imax(ly - 1, 0), lz);
647 0 : idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
648 :
649 0 : jco_l = coset(lx + 1, imax(ly - 1, 0), lz);
650 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
651 :
652 0 : jco_l = coset(imax(lx - 1, 0), ly + 1, lz);
653 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
654 :
655 0 : jco_l = coset(lx + 1, ly + 1, lz);
656 0 : idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
657 0 : } break;
658 0 : case 'Y': {
659 0 : const int l1 = ly;
660 0 : const int l2 = lz;
661 0 : jco_l = coset(lx, imax(ly - 1, 0), imax(lz - 1, 0));
662 0 : idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
663 :
664 0 : jco_l = coset(lx, ly + 1, imax(lz - 1, 0));
665 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
666 :
667 0 : jco_l = coset(lx, imax(ly - 1, 0), lz + 1);
668 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
669 :
670 0 : jco_l = coset(lx, ly + 1, lz + 1);
671 0 : idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
672 0 : } break;
673 0 : case 'Z': {
674 0 : const int l1 = lz;
675 0 : const int l2 = lx;
676 0 : jco_l = coset(imax(lx - 1, 0), ly, imax(lz - 1, 0));
677 0 : idx2(pab_prep[0], jco_l, ico_l) += l1 * l2 * func_a;
678 :
679 0 : jco_l = coset(imax(lx - 1, 0), ly, lz + 1);
680 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l2 * func_a;
681 :
682 0 : jco_l = coset(lx + 1, ly, imax(lz - 1, 0));
683 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * l1 * func_a;
684 :
685 0 : jco_l = coset(lx + 1, ly, lz + 1);
686 0 : idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
687 0 : } break;
688 : default:
689 : break;
690 : }
691 0 : }
692 : // *****************************************************************************
693 0 : static void grid_prepare_pab_DiDj(struct pab_computation_struct_ *const tp) {
694 : // create a new pab_local so that mapping pab_local with pgf_a pgf_b
695 : // is equivalent to mapping pab with
696 : // d_{ider1} pgf_a d_{ider1} pgf_b
697 :
698 0 : for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
699 0 : for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
700 0 : for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
701 0 : for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
702 0 : for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
703 0 : lza <= tp->lmax[0] - lxa - lya; lza++) {
704 0 : for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
705 0 : lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
706 0 : const int ico = tp->offset[0] + coset(lxa, lya, lza);
707 0 : const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
708 0 : const double pab = idx2(tp->pab[0], jco, ico);
709 :
710 0 : int ico_l;
711 0 : double func_a;
712 :
713 : // this element of pab results in 12 elements of pab_local
714 :
715 0 : if ((tp->dir1 == 'X' && tp->dir2 == 'Y') ||
716 : (tp->dir1 == 'Y' && tp->dir2 == 'X')) {
717 : // xy
718 0 : ico_l = coset(imax(lxa - 1, 0), imax(lya - 1, 0), lza);
719 0 : func_a = lxa * lya * pab;
720 0 : oneterm_dijdij('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
721 : tp->pab_prep);
722 :
723 0 : ico_l = coset(lxa + 1, imax(lya - 1, 0), lza);
724 0 : func_a = -2.0 * tp->zeta[0] * lya * pab;
725 0 : oneterm_dijdij('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
726 : tp->pab_prep);
727 :
728 0 : ico_l = coset(imax(lxa - 1, 0), lya + 1, lza);
729 0 : func_a = -2.0 * tp->zeta[0] * lxa * pab;
730 0 : oneterm_dijdij('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
731 : tp->pab_prep);
732 :
733 0 : ico_l = coset(lxa + 1, lya + 1, lza);
734 0 : func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
735 0 : oneterm_dijdij('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
736 : tp->pab_prep);
737 0 : } else if ((tp->dir1 == 'Y' && tp->dir2 == 'Z') ||
738 : (tp->dir1 == 'Z' && tp->dir2 == 'Y')) {
739 : // yz
740 0 : ico_l = coset(lxa, imax(lya - 1, 0), imax(lza - 1, 0));
741 0 : func_a = lya * lza * pab;
742 0 : oneterm_dijdij('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
743 : tp->pab_prep);
744 :
745 0 : ico_l = coset(lxa, lya + 1, imax(lza - 1, 0));
746 0 : func_a = -2.0 * tp->zeta[0] * lza * pab;
747 0 : oneterm_dijdij('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
748 : tp->pab_prep);
749 :
750 0 : ico_l = coset(lxa, imax(lya - 1, 0), lza + 1);
751 0 : func_a = -2.0 * tp->zeta[0] * lya * pab;
752 0 : oneterm_dijdij('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
753 : tp->pab_prep);
754 :
755 0 : ico_l = coset(lxa, lya + 1, lza + 1);
756 0 : func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
757 0 : oneterm_dijdij(2, func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
758 : tp->pab_prep);
759 0 : } else if ((tp->dir1 == 'Z' && tp->dir2 == 'X') ||
760 : (tp->dir1 == 'X' && tp->dir2 == 'Z')) {
761 : // zx
762 0 : ico_l = coset(imax(lxa - 1, 0), lya, imax(lza - 1, 0));
763 0 : func_a = lza * lxa * pab;
764 0 : oneterm_dijdij('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
765 : tp->pab_prep);
766 :
767 0 : ico_l = coset(imax(lxa - 1, 0), lya, lza + 1);
768 0 : func_a = -2.0 * tp->zeta[0] * lxa * pab;
769 0 : oneterm_dijdij('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
770 : tp->pab_prep);
771 :
772 0 : ico_l = coset(lxa + 1, lya, imax(lza - 1, 0));
773 0 : func_a = -2.0 * tp->zeta[0] * lza * pab;
774 0 : oneterm_dijdij('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
775 : tp->pab_prep);
776 :
777 0 : ico_l = coset(lxa + 1, lya, lza + 1);
778 0 : func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
779 0 : oneterm_dijdij('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
780 : tp->pab_prep);
781 : }
782 : }
783 : }
784 : }
785 : }
786 : }
787 : }
788 0 : }
789 :
790 : // *****************************************************************************
791 0 : static void oneterm_diidii(const int idir, const double func_a, const int ico_l,
792 : const int lx, const int ly, const int lz,
793 : const double zet, tensor *const pab_prep) {
794 0 : int jco_l;
795 :
796 0 : switch (idir) {
797 0 : case 'X': {
798 0 : const int l1 = lx;
799 0 : jco_l = coset(imax(lx - 2, 0), ly, lz);
800 0 : idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
801 :
802 0 : jco_l = coset(lx, ly, lz);
803 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
804 :
805 0 : jco_l = coset(lx + 2, ly, lz);
806 0 : idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
807 0 : } break;
808 0 : case 'Y': {
809 0 : const int l1 = ly;
810 0 : jco_l = coset(lx, imax(ly - 2, 0), lz);
811 0 : idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
812 :
813 0 : jco_l = coset(lx, ly, lz);
814 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
815 :
816 0 : jco_l = coset(lx, ly + 2, lz);
817 0 : idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
818 0 : } break;
819 0 : case 'Z': {
820 0 : const int l1 = lz;
821 0 : jco_l = coset(lx, ly, imax(lz - 2, 0));
822 0 : idx2(pab_prep[0], jco_l, ico_l) += l1 * (l1 - 1) * func_a;
823 :
824 0 : jco_l = coset(lx, ly, lz);
825 0 : idx2(pab_prep[0], jco_l, ico_l) -= 2.0 * zet * (2 * l1 + 1) * func_a;
826 :
827 0 : jco_l = coset(lx, ly, lz + 2);
828 0 : idx2(pab_prep[0], jco_l, ico_l) += 4.0 * zet * zet * func_a;
829 0 : } break;
830 : default:
831 0 : printf("Wrong value for ider: should be 1, 2, or 3\n");
832 0 : abort();
833 0 : break;
834 : }
835 0 : }
836 :
837 : // *****************************************************************************
838 0 : static void grid_prepare_pab_Di2(struct pab_computation_struct_ *const tp) {
839 : // create a new pab_local so that mapping pab_local with pgf_a pgf_b
840 : // is equivalent to mapping pab with
841 : // dd_{ider1} pgf_a dd_{ider1} pgf_b
842 :
843 0 : for (int lxa = 0; lxa <= tp->lmax[0]; lxa++) {
844 0 : for (int lxb = 0; lxb <= tp->lmax[1]; lxb++) {
845 0 : for (int lya = 0; lya <= tp->lmax[0] - lxa; lya++) {
846 0 : for (int lyb = 0; lyb <= tp->lmax[1] - lxb; lyb++) {
847 0 : for (int lza = imax(tp->lmin[0] - lxa - lya, 0);
848 0 : lza <= tp->lmax[0] - lxa - lya; lza++) {
849 0 : for (int lzb = imax(tp->lmin[1] - lxb - lyb, 0);
850 0 : lzb <= tp->lmax[1] - lxb - lyb; lzb++) {
851 0 : const int ico = tp->offset[0] + coset(lxa, lya, lza);
852 0 : const int jco = tp->offset[1] + coset(lxb, lyb, lzb);
853 0 : const double pab = idx2(tp->pab[0], jco, ico);
854 :
855 0 : int ico_l;
856 0 : double func_a;
857 :
858 : // this element of pab results in 9 elements of pab_local
859 0 : switch (tp->dir1) {
860 0 : case 'X': {
861 : // x
862 0 : ico_l = coset(imax(lxa - 2, 0), lya, lza);
863 0 : func_a = lxa * (lxa - 1) * pab;
864 0 : oneterm_diidii('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
865 : tp->pab_prep);
866 :
867 0 : ico_l = coset(lxa, lya, lza);
868 0 : func_a = -2.0 * tp->zeta[0] * (2 * lxa + 1) * pab;
869 0 : oneterm_diidii('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
870 : tp->pab_prep);
871 :
872 0 : ico_l = coset(lxa + 2, lya, lza);
873 0 : func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
874 0 : oneterm_diidii('X', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
875 : tp->pab_prep);
876 0 : } break;
877 0 : case 'Y': {
878 : // y
879 0 : ico_l = coset(lxa, imax(lya - 2, 0), lza);
880 0 : func_a = lya * (lya - 1) * pab;
881 0 : oneterm_diidii('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
882 : tp->pab_prep);
883 :
884 0 : ico_l = coset(lxa, lya, lza);
885 0 : func_a = -2.0 * tp->zeta[0] * (2 * lya + 1) * pab;
886 0 : oneterm_diidii('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
887 : tp->pab_prep);
888 :
889 0 : ico_l = coset(lxa, lya + 2, lza);
890 0 : func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
891 0 : oneterm_diidii('Y', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
892 : tp->pab_prep);
893 0 : } break;
894 0 : case 'Z': {
895 : // z
896 0 : ico_l = coset(lxa, lya, imax(lza - 2, 0));
897 0 : func_a = lza * (lza - 1) * pab;
898 0 : oneterm_diidii('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
899 : tp->pab_prep);
900 :
901 0 : ico_l = coset(lxa, lya, lza);
902 0 : func_a = -2.0 * tp->zeta[0] * (2 * lza + 1) * pab;
903 0 : oneterm_diidii('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
904 : tp->pab_prep);
905 :
906 0 : ico_l = coset(lxa, lya, lza + 2);
907 0 : func_a = 4.0 * tp->zeta[0] * tp->zeta[0] * pab;
908 0 : oneterm_diidii('Z', func_a, ico_l, lxb, lyb, lzb, tp->zeta[1],
909 : tp->pab_prep);
910 0 : } break;
911 : default:
912 0 : printf("Wrong value for ider: should be 'X', 'Y' or 'Z'.\n");
913 0 : abort();
914 0 : break;
915 : }
916 : }
917 : }
918 : }
919 : }
920 : }
921 : }
922 0 : }
923 :
924 : // *****************************************************************************
925 632 : void grid_prepare_get_ldiffs_dgemm(const enum grid_func func,
926 : int *const lmin_diff, int *const lmax_diff) {
927 632 : switch (func) {
928 632 : case GRID_FUNC_AB:
929 632 : lmax_diff[0] = 0;
930 632 : lmin_diff[0] = 0;
931 632 : lmax_diff[1] = 0;
932 632 : lmin_diff[1] = 0;
933 632 : break;
934 0 : case GRID_FUNC_DADB:
935 : case GRID_FUNC_ADBmDAB_X:
936 : case GRID_FUNC_ADBmDAB_Y:
937 : case GRID_FUNC_ADBmDAB_Z:
938 : case GRID_FUNC_DABpADB_X:
939 : case GRID_FUNC_DABpADB_Y:
940 : case GRID_FUNC_DABpADB_Z:
941 0 : lmax_diff[0] = 1;
942 0 : lmin_diff[0] = -1;
943 0 : lmax_diff[1] = 1;
944 0 : lmin_diff[1] = -1;
945 0 : break;
946 0 : case GRID_FUNC_ARDBmDARB_XX:
947 : case GRID_FUNC_ARDBmDARB_XY:
948 : case GRID_FUNC_ARDBmDARB_XZ:
949 : case GRID_FUNC_ARDBmDARB_YX:
950 : case GRID_FUNC_ARDBmDARB_YY:
951 : case GRID_FUNC_ARDBmDARB_YZ:
952 : case GRID_FUNC_ARDBmDARB_ZX:
953 : case GRID_FUNC_ARDBmDARB_ZY:
954 : case GRID_FUNC_ARDBmDARB_ZZ:
955 0 : lmax_diff[0] = 1; // TODO: mistake???, then we could merge la and lb.
956 0 : lmin_diff[0] = -1;
957 0 : lmax_diff[1] = 2;
958 0 : lmin_diff[1] = -1;
959 0 : break;
960 0 : case GRID_FUNC_DX:
961 : case GRID_FUNC_DY:
962 : case GRID_FUNC_DZ:
963 0 : lmax_diff[0] = 1;
964 0 : lmin_diff[0] = -1;
965 0 : lmax_diff[1] = 1;
966 0 : lmin_diff[1] = -1;
967 0 : break;
968 0 : case GRID_FUNC_DXDY:
969 : case GRID_FUNC_DYDZ:
970 : case GRID_FUNC_DZDX:
971 : case GRID_FUNC_DXDX:
972 : case GRID_FUNC_DYDY:
973 : case GRID_FUNC_DZDZ:
974 0 : lmax_diff[0] = 2;
975 0 : lmin_diff[0] = -2;
976 0 : lmax_diff[1] = 2;
977 0 : lmin_diff[1] = -2;
978 0 : break;
979 : default:
980 0 : printf("Unkown ga-gb function");
981 0 : abort();
982 : }
983 632 : }
984 :
985 : // *****************************************************************************
986 44389 : void grid_prepare_pab_dgemm(const enum grid_func func, const int *const offset,
987 : const int *const lmin, const int *const lmax,
988 : const double *const zeta, tensor *const pab,
989 : tensor *const pab_prep) {
990 44389 : struct pab_computation_struct_ tmp;
991 :
992 44389 : tmp.offset[0] = offset[0];
993 44389 : tmp.offset[1] = offset[1];
994 :
995 44389 : tmp.lmin[0] = lmin[0];
996 44389 : tmp.lmin[1] = lmin[1];
997 :
998 44389 : tmp.lmax[0] = lmax[0];
999 44389 : tmp.lmax[1] = lmax[1];
1000 :
1001 44389 : tmp.pab = pab;
1002 44389 : tmp.pab_prep = pab_prep;
1003 :
1004 44389 : tmp.zeta[0] = zeta[0];
1005 44389 : tmp.zeta[1] = zeta[1];
1006 44389 : memset(pab_prep->data, 0, pab_prep->alloc_size_ * sizeof(double));
1007 :
1008 44389 : switch (func) {
1009 44389 : case GRID_FUNC_AB:
1010 44389 : grid_prepare_pab_AB(&tmp);
1011 44389 : break;
1012 0 : case GRID_FUNC_DADB:
1013 0 : grid_prepare_pab_DADB(&tmp);
1014 0 : break;
1015 0 : case GRID_FUNC_ADBmDAB_X:
1016 0 : tmp.dir1 = 'X';
1017 0 : grid_prepare_pab_ADBmDAB(&tmp);
1018 0 : break;
1019 0 : case GRID_FUNC_ADBmDAB_Y:
1020 0 : tmp.dir1 = 'Y';
1021 0 : grid_prepare_pab_ADBmDAB(&tmp);
1022 0 : break;
1023 0 : case GRID_FUNC_ADBmDAB_Z:
1024 0 : tmp.dir1 = 'Z';
1025 0 : grid_prepare_pab_ADBmDAB(&tmp);
1026 0 : break;
1027 0 : case GRID_FUNC_ARDBmDARB_XX:
1028 0 : tmp.dir1 = 'X';
1029 0 : tmp.dir2 = 'X';
1030 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1031 0 : break;
1032 0 : case GRID_FUNC_ARDBmDARB_XY:
1033 0 : tmp.dir1 = 'X';
1034 0 : tmp.dir2 = 'Y';
1035 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1036 0 : break;
1037 0 : case GRID_FUNC_ARDBmDARB_XZ:
1038 0 : tmp.dir1 = 'X';
1039 0 : tmp.dir2 = 'Z';
1040 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1041 0 : break;
1042 0 : case GRID_FUNC_ARDBmDARB_YX:
1043 0 : tmp.dir1 = 'Y';
1044 0 : tmp.dir2 = 'X';
1045 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1046 0 : break;
1047 0 : case GRID_FUNC_ARDBmDARB_YY:
1048 0 : tmp.dir1 = 'Y';
1049 0 : tmp.dir2 = 'Y';
1050 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1051 0 : break;
1052 0 : case GRID_FUNC_ARDBmDARB_YZ:
1053 0 : tmp.dir1 = 'Y';
1054 0 : tmp.dir2 = 'Z';
1055 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1056 0 : break;
1057 0 : case GRID_FUNC_ARDBmDARB_ZX:
1058 0 : tmp.dir1 = 'Z';
1059 0 : tmp.dir2 = 'X';
1060 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1061 0 : break;
1062 0 : case GRID_FUNC_ARDBmDARB_ZY:
1063 0 : tmp.dir1 = 'Z';
1064 0 : tmp.dir2 = 'Y';
1065 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1066 0 : break;
1067 0 : case GRID_FUNC_ARDBmDARB_ZZ:
1068 0 : tmp.dir1 = 'Z';
1069 0 : tmp.dir2 = 'Z';
1070 0 : grid_prepare_pab_ARDBmDARB(&tmp);
1071 0 : break;
1072 0 : case GRID_FUNC_DABpADB_X:
1073 0 : tmp.dir1 = 'X';
1074 0 : grid_prepare_pab_DABpADB(&tmp);
1075 0 : break;
1076 0 : case GRID_FUNC_DABpADB_Y:
1077 0 : tmp.dir1 = 'Y';
1078 0 : grid_prepare_pab_DABpADB(&tmp);
1079 0 : break;
1080 0 : case GRID_FUNC_DABpADB_Z:
1081 0 : tmp.dir1 = 'Z';
1082 0 : grid_prepare_pab_DABpADB(&tmp);
1083 0 : break;
1084 0 : case GRID_FUNC_DX:
1085 0 : tmp.dir1 = 'X';
1086 0 : grid_prepare_pab_Di(&tmp);
1087 0 : break;
1088 0 : case GRID_FUNC_DY:
1089 0 : tmp.dir1 = 'Y';
1090 0 : grid_prepare_pab_Di(&tmp);
1091 0 : break;
1092 0 : case GRID_FUNC_DZ:
1093 0 : tmp.dir1 = 'Z';
1094 0 : grid_prepare_pab_Di(&tmp);
1095 0 : break;
1096 0 : case GRID_FUNC_DXDY:
1097 0 : tmp.dir1 = 'X';
1098 0 : tmp.dir2 = 'Y';
1099 0 : grid_prepare_pab_DiDj(&tmp);
1100 0 : break;
1101 0 : case GRID_FUNC_DYDZ:
1102 0 : tmp.dir1 = 'Y';
1103 0 : tmp.dir2 = 'Z';
1104 0 : grid_prepare_pab_DiDj(&tmp);
1105 0 : break;
1106 0 : case GRID_FUNC_DZDX:
1107 0 : tmp.dir1 = 'Z';
1108 0 : tmp.dir2 = 'X';
1109 0 : grid_prepare_pab_DiDj(&tmp);
1110 0 : break;
1111 0 : case GRID_FUNC_DXDX:
1112 0 : tmp.dir1 = 'X';
1113 0 : grid_prepare_pab_Di2(&tmp);
1114 0 : break;
1115 0 : case GRID_FUNC_DYDY:
1116 0 : tmp.dir1 = 'Y';
1117 0 : grid_prepare_pab_Di2(&tmp);
1118 0 : break;
1119 0 : case GRID_FUNC_DZDZ:
1120 0 : tmp.dir1 = 'Z';
1121 0 : grid_prepare_pab_Di2(&tmp);
1122 0 : break;
1123 : default:
1124 0 : printf("Unkown ga-gb function");
1125 0 : abort();
1126 : }
1127 44389 : }
|
