LIBINT  2.1.0-stable
VRR_r12kg12_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_vrrr12kg12xsxs_h_
20 #define _libint2_src_lib_libint_vrrr12kg12xsxs_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, int K, bool vectorize> struct VRR_r12kg12_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  const LIBINT2_REALTYPE* src5);
38  };
39 
50  template <int La, int Lc, int K, bool vectorize> struct VRR_r12kg12_xs_xs<0,La,Lc,K,vectorize> {
51 
52  static void compute(const Libint_t* inteval,
53  LIBINT2_REALTYPE* target,
54  const LIBINT2_REALTYPE* src0,
55  const LIBINT2_REALTYPE* src1,
56  const LIBINT2_REALTYPE* src2,
57  const LIBINT2_REALTYPE* src3,
58  const LIBINT2_REALTYPE* src4,
59  const LIBINT2_REALTYPE* src5) {
60 
61  // works for (ds|ps) and higher, K >= 0
62  if ((La < 2 && Lc < 1) || K < 0)
63  abort();
64 
65  const unsigned int veclen = vectorize ? inteval->veclen : 1;
66 
67  const unsigned int Nc = INT_NCART(Lc);
68  const unsigned int NcV = Nc * veclen;
69  const unsigned int Ncm1 = INT_NCART(Lc-1);
70  const unsigned int Ncm1V = Ncm1 * veclen;
71  const unsigned int Ncp1 = INT_NCART(Lc+1);
72  const unsigned int Ncp1V = Ncp1 * veclen;
73 
74  int ax, ay, az;
75  FOR_CART(ax, ay, az, La)
76 
77  int a[3]; a[0] = ax; a[1] = ay; a[2] = az;
78 
79  enum XYZ {x=0, y=1, z=2};
80  // Build along x, if possible
81  XYZ xyz = z;
82  if (ay != 0) xyz = y;
83  if (ax != 0) xyz = x;
84  --a[xyz];
85 
86  // redirect
87  const LIBINT2_REALTYPE *pfac0_0, *pfac4_0;
88  switch(xyz) {
89  case x:
90  pfac0_0 = inteval->R12kG12_pfac0_0_x;
91  pfac4_0 = inteval->R12kG12_pfac4_0_x;
92  break;
93  case y:
94  pfac0_0 = inteval->R12kG12_pfac0_0_y;
95  pfac4_0 = inteval->R12kG12_pfac4_0_y;
96  break;
97  case z:
98  pfac0_0 = inteval->R12kG12_pfac0_0_z;
99  pfac4_0 = inteval->R12kG12_pfac4_0_z;
100  break;
101  }
102 
103  const unsigned int iam1 = INT_CARTINDEX(La-1,a[0],a[1]);
104  const unsigned int am10c0_offset = iam1 * NcV;
105  const LIBINT2_REALTYPE* src0_ptr = src0 + am10c0_offset;
106 
107  // if a-2_xyz exists, include (a-2_xyz 0 | c 0)
108  if (a[xyz] > 0) {
109  --a[xyz];
110  const unsigned int iam2 = INT_CARTINDEX(La-2,a[0],a[1]);
111  const unsigned int am20c0_offset = iam2 * NcV;
112  ++a[xyz];
113  const LIBINT2_REALTYPE* src1_ptr = src1 + 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  target[cv] = pfac0_0[v] * src0_ptr[cv] + axyz * inteval->R12kG12_pfac1_0[v] * src1_ptr[cv];
120  }
121  }
122 #if LIBINT2_FLOP_COUNT
123  inteval->nflops[0] += 4 * NcV;
124 #endif
125 
126  }
127  else {
128  unsigned int cv = 0;
129  for(unsigned int c = 0; c < Nc; ++c) {
130  for(unsigned int v=0; v<veclen; ++v, ++cv) {
131  target[cv] = pfac0_0[v] * src0_ptr[cv];
132  }
133  }
134 #if LIBINT2_FLOP_COUNT
135  inteval->nflops[0] += NcV;
136 #endif
137  }
138 
139  {
140  const unsigned int am10cm10_offset = iam1 * Ncm1V;
141  const LIBINT2_REALTYPE* src2_ptr = src2 + am10cm10_offset;
142 
143  // loop over c-1 shell and include (a-1_xyz 0 | c-1_xyz 0) to (a 0 | c 0)
144  int cx, cy, cz;
145  FOR_CART(cx, cy, cz, Lc-1)
146 
147  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
148  ++c[xyz];
149 
150  const unsigned int cc = INT_CARTINDEX(Lc,c[0],c[1]);
151  const unsigned int cc_offset = cc * veclen;
152  LIBINT2_REALTYPE* tptr = target + cc_offset;
153  const LIBINT2_REALTYPE cxyz = (LIBINT2_REALTYPE)c[xyz];
154  for(unsigned int v=0; v<veclen; ++v) {
155  tptr[v] += cxyz * inteval->R12kG12_pfac2[v] * src2_ptr[v];
156  }
157 #if LIBINT2_FLOP_COUNT
158  inteval->nflops[0] += 3 * veclen;
159 #endif
160  src2_ptr += veclen;
161 
162  END_FOR_CART
163  }
164 
165  if (K > 0) {
166  ++a[xyz];
167  const unsigned int ia = INT_CARTINDEX(La,a[0],a[1]);
168  const unsigned int a0c0_offset = ia * NcV;
169  --a[xyz];
170  const LIBINT2_REALTYPE* src3_ptr = src3 + a0c0_offset;
171  const LIBINT2_REALTYPE* src5_ptr = src5 + am10c0_offset;
172  unsigned int cv = 0;
173  for(unsigned int c = 0; c < Nc; ++c) {
174  for(unsigned int v=0; v<veclen; ++v, ++cv) {
175  target[cv] += K * inteval->R12kG12_pfac3_0[v] * (src3_ptr[cv] + pfac4_0[v] * src5_ptr[cv]);
176  }
177  }
178 #if LIBINT2_FLOP_COUNT
179  inteval->nflops[0] += 5 * NcV;
180 #endif
181 
182  // loop over c shell and include (a-1_xyz 0 | c+1_xyz 0) to (a 0 | c 0)
183  int cx, cy, cz;
184  cv = 0;
185  const unsigned int am10cp10_offset = iam1 * Ncp1V;
186  const LIBINT2_REALTYPE* src4_ptr0 = src4 + am10cp10_offset;
187  FOR_CART(cx, cy, cz, Lc)
188 
189  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
190  ++c[xyz];
191 
192  const unsigned int cc = INT_CARTINDEX(Lc+1,c[0],c[1]);
193  const unsigned int cc_offset = cc * veclen;
194  const LIBINT2_REALTYPE* src4_ptr = src4_ptr0 + cc_offset;
195  for(unsigned int v=0; v<veclen; ++v, ++cv) {
196  target[cv] -= K * inteval->R12kG12_pfac3_0[v] * src4_ptr[v];
197  }
198 #if LIBINT2_FLOP_COUNT
199  inteval->nflops[0] += 3 * veclen;
200 #endif
201 
202  END_FOR_CART
203 
204  }
205 
206  target += NcV;
207 
208  END_FOR_CART
209 
211  //inteval->nflops[0] = inteval->nflops[0] + 222 * 1 * 1 * veclen;
212 
213  }
214 
215  };
216 
217 };
218 
219 #endif // header guard
220 
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, const LIBINT2_REALTYPE *src5)
Definition: VRR_r12kg12_xs_xs.h:52
Definition: VRR_r12kg12_xs_xs.h:29
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
Definition: test_eri_rys.cc:46