LIBINT  2.1.0-stable
ITR_xs_xs.h
1 /*
2  * This file is a part of Libint.
3  * Copyright (C) 2004-2014 Edward F. Valeev
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU Library General Public License, version 2,
7  * as published by the Free Software Foundation.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU Library General Public License
15  * along with this program. If not, see http://www.gnu.org/licenses/.
16  *
17  */
18 
19 #ifndef _libint2_src_lib_libint_itrxsxs_h_
20 #define _libint2_src_lib_libint_itrxsxs_h_
21 
22 #include <cstdlib>
23 #include <libint2.h>
24 #include <util_types.h>
25 #include <libint2/cgshell_ordering.h>
26 
27 namespace libint2 {
28 
29  template <int part, int La, int Lc, bool InBra, bool vectorize> struct ITR_xs_xs {
30  static void compute(const Libint_t* inteval,
31  LIBINT2_REALTYPE* target,
32  const LIBINT2_REALTYPE* src0,
33  const LIBINT2_REALTYPE* src1,
34  const LIBINT2_REALTYPE* src2,
35  const LIBINT2_REALTYPE* src3);
36  };
37 
44  template <int La, int Lc, bool InBra, bool vectorize> struct ITR_xs_xs<0,La,Lc,InBra,vectorize> {
45 
46  static void compute(const Libint_t* inteval,
47  LIBINT2_REALTYPE* target,
48  const LIBINT2_REALTYPE* src0,
49  const LIBINT2_REALTYPE* src1,
50  const LIBINT2_REALTYPE* src2,
51  const LIBINT2_REALTYPE* src3) {
52 
53  // works for (ds|ps) and higher
54  if (La < 2 || Lc < 1)
55  abort();
56 
57  const unsigned int veclen = vectorize ? inteval->veclen : 1;
58 
59  const unsigned int Nc = INT_NCART(Lc);
60  const unsigned int NcV = Nc * veclen;
61 
62  int ax, ay, az;
63  FOR_CART(ax, ay, az, La)
64 
65  int a[3]; a[0] = ax; a[1] = ay; a[2] = az;
66 
67  enum XYZ {x=0, y=1, z=2};
68  // Build along x, if possible
69  XYZ xyz = z;
70  if (ay != 0) xyz = y;
71  if (ax != 0) xyz = x;
72  --a[xyz];
73 
74  // redirect
75  const LIBINT2_REALTYPE *pfac0;
76  switch(xyz) {
77  case x:
78  pfac0 = InBra ?
79 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_x)
80  inteval->TwoPRepITR_pfac0_0_0_x
81 #else
82  0
83 #endif
84  :
85 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x)
86  inteval->TwoPRepITR_pfac0_0_1_x;
87 #else
88  0
89 #endif
90  ;
91  break;
92  case y:
93  pfac0 = InBra ?
94 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_y)
95  inteval->TwoPRepITR_pfac0_0_0_y
96 #else
97  0
98 #endif
99  :
100 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_y)
101  inteval->TwoPRepITR_pfac0_0_1_y;
102 #else
103  0
104 #endif
105  ;
106  break;
107  case z:
108  pfac0 = InBra ?
109 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_0_z)
110  inteval->TwoPRepITR_pfac0_0_0_z
111 #else
112  0
113 #endif
114  :
115 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_z)
116  inteval->TwoPRepITR_pfac0_0_1_z;
117 #else
118  0
119 #endif
120  ;
121  break;
122  }
123 
124  const unsigned int iam1 = INT_CARTINDEX(La-1,a[0],a[1]);
125  const unsigned int am10c0_offset = iam1 * NcV;
126  const LIBINT2_REALTYPE* src0_ptr = src0 + am10c0_offset;
127 
128  // if a-2_xyz exists, include (a-2_xyz 0 | c 0)
129  if (a[xyz] > 0) {
130  --a[xyz];
131  const unsigned int iam2 = INT_CARTINDEX(La-2,a[0],a[1]);
132  const unsigned int am20c0_offset = iam2 * NcV;
133  ++a[xyz];
134  const LIBINT2_REALTYPE* src2_ptr = src2 + am20c0_offset;
135  const LIBINT2_REALTYPE axyz = (LIBINT2_REALTYPE)a[xyz];
136 
137  unsigned int cv = 0;
138  for(unsigned int c = 0; c < Nc; ++c) {
139  for(unsigned int v=0; v<veclen; ++v, ++cv) {
140  target[cv] = pfac0[v] * src0_ptr[cv] + axyz * inteval->oo2z[v] * src2_ptr[cv];
141  }
142  }
143 #if LIBINT2_FLOP_COUNT
144  inteval->nflops[0] += 4 * NcV;
145 #endif
146 
147  }
148  else {
149  unsigned int cv = 0;
150  for(unsigned int c = 0; c < Nc; ++c) {
151  for(unsigned int v=0; v<veclen; ++v, ++cv) {
152  target[cv] = pfac0[v] * src0_ptr[cv];
153  }
154  }
155 #if LIBINT2_FLOP_COUNT
156  inteval->nflops[0] += 1 * NcV;
157 #endif
158  }
159 
160  {
161  const unsigned int Ncp1 = INT_NCART(Lc+1);
162  const unsigned int Ncp1V = Ncp1 * veclen;
163  const unsigned int am10cp10_offset = iam1 * Ncp1V;
164  const LIBINT2_REALTYPE* src1_ptr = src1 + am10cp10_offset;
165 
166  // loop over c shell and include (a-1_xyz 0 | c+1_xyz 0) to (a 0 | c 0)
167  int cx, cy, cz;
168  int cv = 0;
169  FOR_CART(cx, cy, cz, Lc)
170 
171  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
172  ++c[xyz];
173 
174  const unsigned int cp1 = INT_CARTINDEX(Lc+1,c[0],c[1]);
175  const unsigned int cp1_offset = cp1 * veclen;
176  const LIBINT2_REALTYPE* sptr = src1_ptr + cp1_offset;
177  const LIBINT2_REALTYPE cxyz = (LIBINT2_REALTYPE)c[xyz];
178  for(unsigned int v=0; v<veclen; ++v, ++cv) {
179  target[cv] -= inteval->eoz[v] * sptr[v];
180  }
181 #if LIBINT2_FLOP_COUNT
182  inteval->nflops[0] += 2 * veclen;
183 #endif
184 
185  END_FOR_CART
186  }
187 
188 
189  {
190  const unsigned int Ncm1 = INT_NCART(Lc-1);
191  const unsigned int Ncm1V = Ncm1 * veclen;
192  const unsigned int am10cm10_offset = iam1 * Ncm1V;
193  const LIBINT2_REALTYPE* src3_ptr = src3 + am10cm10_offset;
194 
195  // loop over c-1 shell and include (a-1_xyz 0 | c-1_xyz 0) to (a 0 | c 0)
196  int cx, cy, cz;
197  FOR_CART(cx, cy, cz, Lc-1)
198 
199  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
200  ++c[xyz];
201 
202  const unsigned int cc = INT_CARTINDEX(Lc,c[0],c[1]);
203  const unsigned int cc_offset = cc * veclen;
204  LIBINT2_REALTYPE* tptr = target + cc_offset;
205  const LIBINT2_REALTYPE cxyz = (LIBINT2_REALTYPE)c[xyz];
206  for(unsigned int v=0; v<veclen; ++v) {
207  tptr[v] += cxyz * inteval->oo2z[v] * src3_ptr[v];
208  }
209 #if LIBINT2_FLOP_COUNT
210  inteval->nflops[0] += 3 * veclen;
211 #endif
212  src3_ptr += veclen;
213 
214  END_FOR_CART
215  }
216 
217  target += NcV;
218 
219  END_FOR_CART
220 
221  }
222 
223  };
224 
225 #if 0
226 
232  template <int La, int Lc, bool InBra, bool vectorize> struct ITR_xs_xs<1,La,Lc,InBra,vectorize> {
233 
234  static void compute(const Libint_t* inteval,
235  LIBINT2_REALTYPE* target,
236  const LIBINT2_REALTYPE* src0,
237  const LIBINT2_REALTYPE* src1,
238  const LIBINT2_REALTYPE* src2,
239  const LIBINT2_REALTYPE* src3) {
240 
241  // works for (ps|ds) and higher
242  if (La < 1 || Lc < 2)
243  abort();
244 
245  const unsigned int veclen = vectorize ? inteval->veclen : 1;
246 
247  const unsigned int Na = INT_NCART(La);
248  const unsigned int NaV = Na * veclen;
249 
250  int cx, cy, cz;
251  FOR_CART(cx, cy, cz, Lc)
252 
253  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
254 
255  enum XYZ {x=0, y=1, z=2};
256  // Build along x, if possible
257  XYZ xyz = z;
258  if (cy != 0) xyz = y;
259  if (cx != 0) xyz = x;
260  --c[xyz];
261 
262  // redirect
263  const LIBINT2_REALTYPE *pfac0;
264  switch(xyz) {
265  case x:
266  pfac0 = InBra ?
267 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_x)
268  inteval->TwoPRepITR_pfac0_1_0_x
269 #else
270  0
271 #endif
272  :
273 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_0_1_x)
274  inteval->TwoPRepITR_pfac0_0_1_x;
275 #else
276  0
277 #endif
278  ;
279  break;
280  case y:
281  pfac0 = InBra ?
282 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_y)
283  inteval->TwoPRepITR_pfac0_1_0_y
284 #else
285  0
286 #endif
287  :
288 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_y)
289  inteval->TwoPRepITR_pfac0_1_1_y;
290 #else
291  0
292 #endif
293  ;
294  break;
295  case z:
296  pfac0 = InBra ?
297 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_0_z)
298  inteval->TwoPRepITR_pfac0_1_0_z
299 #else
300  0
301 #endif
302  :
303 #if LIBINT2_DEFINED(eri,TwoPRepITR_pfac0_1_1_z)
304  inteval->TwoPRepITR_pfac0_1_1_z;
305 #else
306  0
307 #endif
308  ;
309  break;
310  }
311 
312  const unsigned int icm1 = INT_CARTINDEX(Lc-1,c[0],c[1]);
313  const unsigned int a0cm10_offset = icm1 * NaV;
314  const LIBINT2_REALTYPE* src0_ptr = src0 + a0cm10_offset;
315 
316  // if c-2_xyz exists, include (a 0 | c-2_xyz 0)
317  if (c[xyz] > 0) {
318  --c[xyz];
319  const unsigned int icm2 = INT_CARTINDEX(Lc-2,c[0],c[1]);
320  const unsigned int a0cm20_offset = icm2 * NaV;
321  ++c[xyz];
322  const LIBINT2_REALTYPE* src2_ptr = src2 + am20c0_offset;
323  const LIBINT2_REALTYPE axyz = (LIBINT2_REALTYPE)c[xyz];
324 
325  unsigned int cv = 0;
326  for(unsigned int c = 0; c < Nc; ++c) {
327  for(unsigned int v=0; v<veclen; ++v, ++cv) {
328  target[cv] = pfac0[v] * src0_ptr[cv] + axyz * inteval->oo2z[v] * src2_ptr[cv];
329  }
330  }
331 #if LIBINT2_FLOP_COUNT
332  inteval->nflops[0] += 4 * NcV;
333 #endif
334 
335  }
336  else {
337  unsigned int cv = 0;
338  for(unsigned int c = 0; c < Nc; ++c) {
339  for(unsigned int v=0; v<veclen; ++v, ++cv) {
340  target[cv] = pfac0[v] * src0_ptr[cv];
341  }
342  }
343 #if LIBINT2_FLOP_COUNT
344  inteval->nflops[0] += 1 * NcV;
345 #endif
346  }
347 
348  {
349  const unsigned int Ncp1 = INT_NCART(Lc+1);
350  const unsigned int Ncp1V = Ncp1 * veclen;
351  const unsigned int am10cp10_offset = iam1 * Ncp1V;
352  const LIBINT2_REALTYPE* src1_ptr = src1 + am10cp10_offset;
353 
354  // loop over c shell and include (c-1_xyz 0 | c+1_xyz 0) to (c 0 | c 0)
355  int cx, cy, cz;
356  int cv = 0;
357  FOR_CART(cx, cy, cz, Lc)
358 
359  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
360  ++c[xyz];
361 
362  const unsigned int cp1 = INT_CARTINDEX(Lc+1,c[0],c[1]);
363  const unsigned int cp1_offset = cp1 * veclen;
364  const LIBINT2_REALTYPE* sptr = src1_ptr + cp1_offset;
365  const LIBINT2_REALTYPE cxyz = (LIBINT2_REALTYPE)c[xyz];
366  for(unsigned int v=0; v<veclen; ++v, ++cv) {
367  target[cv] -= inteval->eoz[v] * sptr[v];
368  }
369 #if LIBINT2_FLOP_COUNT
370  inteval->nflops[0] += 2 * veclen;
371 #endif
372 
373  END_FOR_CART
374  }
375 
376 
377  {
378  const unsigned int Ncm1 = INT_NCART(Lc-1);
379  const unsigned int Ncm1V = Ncm1 * veclen;
380  const unsigned int am10cm10_offset = iam1 * Ncm1V;
381  const LIBINT2_REALTYPE* src3_ptr = src3 + am10cm10_offset;
382 
383  // loop over c-1 shell and include (c-1_xyz 0 | c-1_xyz 0) to (c 0 | c 0)
384  int cx, cy, cz;
385  FOR_CART(cx, cy, cz, Lc-1)
386 
387  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
388  ++c[xyz];
389 
390  const unsigned int cc = INT_CARTINDEX(Lc,c[0],c[1]);
391  const unsigned int cc_offset = cc * veclen;
392  LIBINT2_REALTYPE* tptr = target + cc_offset;
393  const LIBINT2_REALTYPE cxyz = (LIBINT2_REALTYPE)c[xyz];
394  for(unsigned int v=0; v<veclen; ++v) {
395  tptr[v] += cxyz * inteval->oo2z[v] * src3_ptr[v];
396  }
397 #if LIBINT2_FLOP_COUNT
398  inteval->nflops[0] += 3 * veclen;
399 #endif
400  src3_ptr += veclen;
401 
402  END_FOR_CART
403  }
404 
405  target += NcV;
406 
407  END_FOR_CART
408 
409  }
410 
411  };
412 #endif
413 
414 };
415 
416 #endif // header guard
417 
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
Definition: test_eri_rys.cc:46
Definition: ITR_xs_xs.h:29