LIBINT  2.1.0-stable
OSVRR_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_osvrrxsxs_h_
20 #define _libint2_src_lib_libint_osvrrxsxs_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 unit_b, bool vectorize> struct OSVRR_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  const LIBINT2_REALTYPE* src4);
37  };
38 
46  template <int La, int Lc,
47  bool unit_b,
48  bool vectorize> struct OSVRR_xs_xs<0,La,Lc,unit_b,vectorize> {
49 
50  static void compute(const Libint_t* inteval,
51  LIBINT2_REALTYPE* target,
52  const LIBINT2_REALTYPE* src0,
53  const LIBINT2_REALTYPE* src1,
54  const LIBINT2_REALTYPE* src2,
55  const LIBINT2_REALTYPE* src3,
56  const LIBINT2_REALTYPE* src4) {
57 
58  // works for (ds|ps) and higher
59  assert(not (La < 2 || Lc < 1));
60 
61  const unsigned int veclen = vectorize ? inteval->veclen : 1;
62 
63  const unsigned int Nc = INT_NCART(Lc);
64  const unsigned int NcV = Nc * veclen;
65 
66  int ax, ay, az;
67  FOR_CART(ax, ay, az, La)
68 
69  int a[3]; a[0] = ax; a[1] = ay; a[2] = az;
70 
71  enum XYZ {x=0, y=1, z=2};
72  // Build along x, if possible
73  XYZ xyz = z;
74  if (ay != 0) xyz = y;
75  if (ax != 0) xyz = x;
76  --a[xyz];
77 
78  // redirect
79  const LIBINT2_REALTYPE *PA, *WP;
80  switch(xyz) {
81  case x:
82 #if LIBINT2_DEFINED(eri,PA_x)
83  if (not unit_b) PA = inteval->PA_x;
84 #endif
85  WP = inteval->WP_x;
86  break;
87  case y:
88 #if LIBINT2_DEFINED(eri,PA_y)
89  if (not unit_b) PA = inteval->PA_y;
90 #endif
91  WP = inteval->WP_y;
92  break;
93  case z:
94 #if LIBINT2_DEFINED(eri,PA_z)
95  if (not unit_b) PA = inteval->PA_z;
96 #endif
97  WP = inteval->WP_z;
98  break;
99  }
100 
101  const unsigned int iam1 = INT_CARTINDEX(La-1,a[0],a[1]);
102  const unsigned int am10c0_offset = iam1 * NcV;
103  const LIBINT2_REALTYPE* src0_ptr = unit_b ? 0 : src0 + am10c0_offset;
104  const LIBINT2_REALTYPE* src1_ptr = src1 + am10c0_offset;
105 
106  // if a-2_xyz exists, include (a-2_xyz 0 | c 0)
107  if (a[xyz] > 0) {
108  --a[xyz];
109  const unsigned int iam2 = INT_CARTINDEX(La-2,a[0],a[1]);
110  const unsigned int am20c0_offset = iam2 * NcV;
111  ++a[xyz];
112  const LIBINT2_REALTYPE* src2_ptr = src2 + am20c0_offset;
113  const LIBINT2_REALTYPE* src3_ptr = src3 + am20c0_offset;
114  const LIBINT2_REALTYPE axyz = (LIBINT2_REALTYPE)a[xyz];
115 
116  unsigned int cv = 0;
117  for(unsigned int c = 0; c < Nc; ++c) {
118  for(unsigned int v=0; v<veclen; ++v, ++cv) {
119  LIBINT2_REALTYPE value = WP[v] * src1_ptr[cv] + axyz * inteval->oo2z[v] * (src2_ptr[cv] - inteval->roz[v] * src3_ptr[cv]);
120  if (not unit_b) value += PA[v] * src0_ptr[cv];
121  target[cv] = value;
122  }
123  }
124 #if LIBINT2_FLOP_COUNT
125  inteval->nflops[0] += (unit_b ? 6 : 8) * NcV;
126 #endif
127 
128  }
129  else {
130  unsigned int cv = 0;
131  for(unsigned int c = 0; c < Nc; ++c) {
132  for(unsigned int v=0; v<veclen; ++v, ++cv) {
133  LIBINT2_REALTYPE value = WP[v] * src1_ptr[cv];
134  if (not unit_b) value += PA[v] * src0_ptr[cv];
135  target[cv] = value;
136  }
137  }
138 #if LIBINT2_FLOP_COUNT
139  inteval->nflops[0] += (unit_b ? 1 : 3) * NcV;
140 #endif
141  }
142 
143  {
144  const unsigned int Ncm1 = INT_NCART(Lc-1);
145  const unsigned int Ncm1V = Ncm1 * veclen;
146  const unsigned int am10cm10_offset = iam1 * Ncm1V;
147  const LIBINT2_REALTYPE* src4_ptr = src4 + am10cm10_offset;
148 
149  // loop over c-1 shell and include (a-1_xyz 0 | c-1_xyz 0) to (a 0 | c 0)
150  int cx, cy, cz;
151  FOR_CART(cx, cy, cz, Lc-1)
152 
153  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
154  ++c[xyz];
155 
156  const unsigned int cc = INT_CARTINDEX(Lc,c[0],c[1]);
157  const unsigned int cc_offset = cc * veclen;
158  LIBINT2_REALTYPE* tptr = target + cc_offset;
159  const LIBINT2_REALTYPE cxyz = (LIBINT2_REALTYPE)c[xyz];
160  for(unsigned int v=0; v<veclen; ++v) {
161  tptr[v] += cxyz * inteval->oo2ze[v] * src4_ptr[v];
162  }
163 #if LIBINT2_FLOP_COUNT
164  inteval->nflops[0] += 3 * veclen;
165 #endif
166  src4_ptr += veclen;
167 
168  END_FOR_CART
169  }
170 
171  target += NcV;
172 
173  END_FOR_CART
174 
176  //inteval->nflops[0] = inteval->nflops[0] + 222 * 1 * 1 * veclen;
177 
178  }
179 
180  };
181 
182  // Ahlrichs' extension of OS VRR
183  template <int part, int La, int Lc, bool vectorize> struct OSAVRR_xs_xs {
184  static void compute(const Libint_t* inteval,
185  LIBINT2_REALTYPE* target,
186  const LIBINT2_REALTYPE* src1,
187  const LIBINT2_REALTYPE* src4);
188  };
189 
194  template <int La, int Lc, bool vectorize> struct OSAVRR_xs_xs<0,La,Lc,vectorize> {
195 
196  static void compute(const Libint_t* inteval,
197  LIBINT2_REALTYPE* target,
198  const LIBINT2_REALTYPE* src1,
199  const LIBINT2_REALTYPE* src4) {
200 
201  // works for (ps|ps) and higher
202  assert(not (La < 1 || Lc < 1));
203 
204  const unsigned int veclen = vectorize ? inteval->veclen : 1;
205 
206  const unsigned int Nc = INT_NCART(Lc);
207  const unsigned int NcV = Nc * veclen;
208 
209  int ax, ay, az;
210  FOR_CART(ax, ay, az, La)
211 
212  int a[3]; a[0] = ax; a[1] = ay; a[2] = az;
213 
214  enum XYZ {x=0, y=1, z=2};
215  // Build along x, if possible
216  XYZ xyz = z;
217  if (ay != 0) xyz = y;
218  if (ax != 0) xyz = x;
219  --a[xyz];
220 
221  // redirect
222  const LIBINT2_REALTYPE *WP;
223  switch(xyz) {
224  case x:
225  WP = inteval->WP_x;
226  break;
227  case y:
228  WP = inteval->WP_y;
229  break;
230  case z:
231  WP = inteval->WP_z;
232  break;
233  }
234 
235  const unsigned int iam1 = INT_CARTINDEX(La-1,a[0],a[1]);
236  const unsigned int am10c0_offset = iam1 * NcV;
237  const LIBINT2_REALTYPE* src1_ptr = src1 + am10c0_offset;
238 
239  {
240  unsigned int cv = 0;
241  for(unsigned int c = 0; c < Nc; ++c) {
242  for(unsigned int v=0; v<veclen; ++v, ++cv) {
243  target[cv] = WP[v] * src1_ptr[cv];
244  }
245  }
246 #if LIBINT2_FLOP_COUNT
247  inteval->nflops[0] += NcV;
248 #endif
249  }
250 
251  {
252  const unsigned int Ncm1 = INT_NCART(Lc-1);
253  const unsigned int Ncm1V = Ncm1 * veclen;
254  const unsigned int am10cm10_offset = iam1 * Ncm1V;
255  const LIBINT2_REALTYPE* src4_ptr = src4 + am10cm10_offset;
256 
257  // loop over c-1 shell and include (a-1_xyz 0 | c-1_xyz 0) to (a 0 | c 0)
258  int cx, cy, cz;
259  FOR_CART(cx, cy, cz, Lc-1)
260 
261  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
262  ++c[xyz];
263 
264  const unsigned int cc = INT_CARTINDEX(Lc,c[0],c[1]);
265  const unsigned int cc_offset = cc * veclen;
266  LIBINT2_REALTYPE* tptr = target + cc_offset;
267  const LIBINT2_REALTYPE cxyz = (LIBINT2_REALTYPE)c[xyz];
268  for(unsigned int v=0; v<veclen; ++v) {
269  tptr[v] += cxyz * inteval->oo2ze[v] * src4_ptr[v];
270  }
271 #if LIBINT2_FLOP_COUNT
272  inteval->nflops[0] += 3 * veclen;
273 #endif
274  src4_ptr += veclen;
275 
276  END_FOR_CART
277  }
278 
279  target += NcV;
280 
281  END_FOR_CART
282 
284  //inteval->nflops[0] = inteval->nflops[0] + 222 * 1 * 1 * veclen;
285 
286  }
287 
288  };
289 
290 };
291 
292 #endif // header guard
static void compute(const Libint_t *inteval, LIBINT2_REALTYPE *target, const LIBINT2_REALTYPE *src1, const LIBINT2_REALTYPE *src4)
Definition: OSVRR_xs_xs.h:196
static void compute(const Libint_t *inteval, LIBINT2_REALTYPE *target, const LIBINT2_REALTYPE *src0, const LIBINT2_REALTYPE *src1, const LIBINT2_REALTYPE *src2, const LIBINT2_REALTYPE *src3, const LIBINT2_REALTYPE *src4)
Definition: OSVRR_xs_xs.h:50
Definition: OSVRR_xs_xs.h:183
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
Definition: test_eri_rys.cc:46
Definition: OSVRR_xs_xs.h:29