LIBINT  2.6.0
VRR_r12kg12_xs_xs.h
1 /*
2  * Copyright (C) 2004-2019 Edward F. Valeev
3  *
4  * This file is part of Libint.
5  *
6  * Libint is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU Lesser General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * Libint is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public License
17  * along with Libint. If not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 #ifndef _libint2_src_lib_libint_vrrr12kg12xsxs_h_
22 #define _libint2_src_lib_libint_vrrr12kg12xsxs_h_
23 
24 #include <cstdlib>
25 #include <libint2.h>
26 #include <util_types.h>
27 #include <libint2/cgshell_ordering.h>
28 
29 namespace libint2 {
30 
31  template <int part, int La, int Lc, int K, bool vectorize> struct VRR_r12kg12_xs_xs {
32  static void compute(const Libint_t* inteval,
33  LIBINT2_REALTYPE* target,
34  const LIBINT2_REALTYPE* src0,
35  const LIBINT2_REALTYPE* src1,
36  const LIBINT2_REALTYPE* src2,
37  const LIBINT2_REALTYPE* src3,
38  const LIBINT2_REALTYPE* src4,
39  const LIBINT2_REALTYPE* src5);
40  };
41 
52  template <int La, int Lc, int K, bool vectorize> struct VRR_r12kg12_xs_xs<0,La,Lc,K,vectorize> {
53 
54  static void compute(const Libint_t* inteval,
55  LIBINT2_REALTYPE* target,
56  const LIBINT2_REALTYPE* src0,
57  const LIBINT2_REALTYPE* src1,
58  const LIBINT2_REALTYPE* src2,
59  const LIBINT2_REALTYPE* src3,
60  const LIBINT2_REALTYPE* src4,
61  const LIBINT2_REALTYPE* src5) {
62 
63  // works for (ds|ps) and higher, K >= 0
64  if ((La < 2 && Lc < 1) || K < 0)
65  abort();
66 
67  const unsigned int veclen = vectorize ? inteval->veclen : 1;
68 
69  const unsigned int Nc = INT_NCART(Lc);
70  const unsigned int NcV = Nc * veclen;
71  const unsigned int Ncm1 = INT_NCART(Lc-1);
72  const unsigned int Ncm1V = Ncm1 * veclen;
73  const unsigned int Ncp1 = INT_NCART(Lc+1);
74  const unsigned int Ncp1V = Ncp1 * veclen;
75 
76  int ax, ay, az;
77  FOR_CART(ax, ay, az, La)
78 
79  int a[3]; a[0] = ax; a[1] = ay; a[2] = az;
80 
81  enum XYZ {x=0, y=1, z=2};
82  // Build along x, if possible
83  XYZ xyz = z;
84  if (ay != 0) xyz = y;
85  if (ax != 0) xyz = x;
86  --a[xyz];
87 
88  // redirect
89  const LIBINT2_REALTYPE *pfac0_0, *pfac4_0;
90  switch(xyz) {
91  case x:
92  pfac0_0 = inteval->R12kG12_pfac0_0_x;
93  pfac4_0 = inteval->R12kG12_pfac4_0_x;
94  break;
95  case y:
96  pfac0_0 = inteval->R12kG12_pfac0_0_y;
97  pfac4_0 = inteval->R12kG12_pfac4_0_y;
98  break;
99  case z:
100  pfac0_0 = inteval->R12kG12_pfac0_0_z;
101  pfac4_0 = inteval->R12kG12_pfac4_0_z;
102  break;
103  }
104 
105  const unsigned int iam1 = INT_CARTINDEX(La-1,a[0],a[1]);
106  const unsigned int am10c0_offset = iam1 * NcV;
107  const LIBINT2_REALTYPE* src0_ptr = src0 + am10c0_offset;
108 
109  // if a-2_xyz exists, include (a-2_xyz 0 | c 0)
110  if (a[xyz] > 0) {
111  --a[xyz];
112  const unsigned int iam2 = INT_CARTINDEX(La-2,a[0],a[1]);
113  const unsigned int am20c0_offset = iam2 * NcV;
114  ++a[xyz];
115  const LIBINT2_REALTYPE* src1_ptr = src1 + am20c0_offset;
116  const LIBINT2_REALTYPE axyz = (LIBINT2_REALTYPE)a[xyz];
117 
118  unsigned int cv = 0;
119  for(unsigned int c = 0; c < Nc; ++c) {
120  for(unsigned int v=0; v<veclen; ++v, ++cv) {
121  target[cv] = pfac0_0[v] * src0_ptr[cv] + axyz * inteval->R12kG12_pfac1_0[v] * src1_ptr[cv];
122  }
123  }
124 #if LIBINT2_FLOP_COUNT
125  inteval->nflops[0] += 4 * 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  target[cv] = pfac0_0[v] * src0_ptr[cv];
134  }
135  }
136 #if LIBINT2_FLOP_COUNT
137  inteval->nflops[0] += NcV;
138 #endif
139  }
140 
141  {
142  const unsigned int am10cm10_offset = iam1 * Ncm1V;
143  const LIBINT2_REALTYPE* src2_ptr = src2 + am10cm10_offset;
144 
145  // loop over c-1 shell and include (a-1_xyz 0 | c-1_xyz 0) to (a 0 | c 0)
146  int cx, cy, cz;
147  FOR_CART(cx, cy, cz, Lc-1)
148 
149  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
150  ++c[xyz];
151 
152  const unsigned int cc = INT_CARTINDEX(Lc,c[0],c[1]);
153  const unsigned int cc_offset = cc * veclen;
154  LIBINT2_REALTYPE* tptr = target + cc_offset;
155  const LIBINT2_REALTYPE cxyz = (LIBINT2_REALTYPE)c[xyz];
156  for(unsigned int v=0; v<veclen; ++v) {
157  tptr[v] += cxyz * inteval->R12kG12_pfac2[v] * src2_ptr[v];
158  }
159 #if LIBINT2_FLOP_COUNT
160  inteval->nflops[0] += 3 * veclen;
161 #endif
162  src2_ptr += veclen;
163 
164  END_FOR_CART
165  }
166 
167  if (K > 0) {
168  ++a[xyz];
169  const unsigned int ia = INT_CARTINDEX(La,a[0],a[1]);
170  const unsigned int a0c0_offset = ia * NcV;
171  --a[xyz];
172  const LIBINT2_REALTYPE* src3_ptr = src3 + a0c0_offset;
173  const LIBINT2_REALTYPE* src5_ptr = src5 + am10c0_offset;
174  unsigned int cv = 0;
175  for(unsigned int c = 0; c < Nc; ++c) {
176  for(unsigned int v=0; v<veclen; ++v, ++cv) {
177  target[cv] += K * inteval->R12kG12_pfac3_0[v] * (src3_ptr[cv] + pfac4_0[v] * src5_ptr[cv]);
178  }
179  }
180 #if LIBINT2_FLOP_COUNT
181  inteval->nflops[0] += 5 * NcV;
182 #endif
183 
184  // loop over c shell and include (a-1_xyz 0 | c+1_xyz 0) to (a 0 | c 0)
185  int cx, cy, cz;
186  cv = 0;
187  const unsigned int am10cp10_offset = iam1 * Ncp1V;
188  const LIBINT2_REALTYPE* src4_ptr0 = src4 + am10cp10_offset;
189  FOR_CART(cx, cy, cz, Lc)
190 
191  int c[3]; c[0] = cx; c[1] = cy; c[2] = cz;
192  ++c[xyz];
193 
194  const unsigned int cc = INT_CARTINDEX(Lc+1,c[0],c[1]);
195  const unsigned int cc_offset = cc * veclen;
196  const LIBINT2_REALTYPE* src4_ptr = src4_ptr0 + cc_offset;
197  for(unsigned int v=0; v<veclen; ++v, ++cv) {
198  target[cv] -= K * inteval->R12kG12_pfac3_0[v] * src4_ptr[v];
199  }
200 #if LIBINT2_FLOP_COUNT
201  inteval->nflops[0] += 3 * veclen;
202 #endif
203 
204  END_FOR_CART
205 
206  }
207 
208  target += NcV;
209 
210  END_FOR_CART
211 
213  //inteval->nflops[0] = inteval->nflops[0] + 222 * 1 * 1 * veclen;
214 
215  }
216 
217  };
218 
219 };
220 
221 #endif // header guard
222 
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:54
Definition: VRR_r12kg12_xs_xs.h:31
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:24