LIBINT  2.1.0-stable
VRR_GTG_1d_xx_xx.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_vrrgtg1dxxxx_h_
20 #define _libint2_src_lib_libint_vrrgtg1dxxxx_h_
21 
22 #include <cstdlib>
23 #include <cassert>
24 #include <libint2.h>
25 #include <util_types.h>
26 
27 namespace libint2 {
28 
36  template <unsigned int CartesianAxis, int La, int Lb, int Lc, int Ld, bool vectorize>
38 
39  static void compute(const Libint_t* inteval,
40  LIBINT2_REALTYPE* target,
41  const LIBINT2_REALTYPE* src0) {
42 
43  enum XYZ {x=0, y=1, z=2};
44  assert(CartesianAxis == x || CartesianAxis == y || CartesianAxis == z);
45  assert(vectorize == false);
46 
47  const unsigned int veclen = vectorize ? inteval->veclen : 1;
48 
49  // corner case: (00|00)
50  if (La == 0 && Lb == 0 && Lc == 0 && Ld == 0) {
51  for (unsigned int v=0; v!=veclen; ++v)
52  target[v] = src0[v];
53  return;
54  }
55 
56  //---------------------------------------------
57  // Part (1): build (a+b 0|c+d 0)
58  //---------------------------------------------
59 
60  LIBINT2_REALTYPE apb_0_GTG_cpd_0[La+Lb+1][Lc+Ld+1];
61  apb_0_GTG_cpd_0[0][0] = src0[0];
62 
63  const LIBINT2_REALTYPE *pfac0_0, *pfac0_1;
64  const LIBINT2_REALTYPE *pfac1_0 = inteval->R12kG12_pfac1_0;
65  const LIBINT2_REALTYPE *pfac1_1 = inteval->R12kG12_pfac1_1;
66  const LIBINT2_REALTYPE *pfac2 = inteval->R12kG12_pfac2;
67  switch (CartesianAxis) {
68  case x:
69  pfac0_0 = inteval->R12kG12_pfac0_0_x;
70  pfac0_1 = inteval->R12kG12_pfac0_1_x;
71  break;
72  case y:
73  pfac0_0 = inteval->R12kG12_pfac0_0_y;
74  pfac0_1 = inteval->R12kG12_pfac0_1_y;
75  break;
76  case z:
77  pfac0_0 = inteval->R12kG12_pfac0_0_z;
78  pfac0_1 = inteval->R12kG12_pfac0_1_z;
79  break;
80  default: assert(false);
81  }
82 
83  // build (0 0|1 0)
84  if (Lc+Ld > 0) {
85  apb_0_GTG_cpd_0[0][1] = pfac0_1[0] * apb_0_GTG_cpd_0[0][0];
86 #if LIBINT2_FLOP_COUNT
87  inteval->nflops[0] += 1;
88 #endif
89  }
90 
91  // build (0 0|c+d 0)
92  if (Lc+Ld > 1) {
93  for(int c_plus_d=1; c_plus_d!=Lc+Ld; ++c_plus_d) {
94  apb_0_GTG_cpd_0[0][c_plus_d+1] = pfac0_1[0] * apb_0_GTG_cpd_0[0][c_plus_d] +
95  c_plus_d * pfac1_1[0] * apb_0_GTG_cpd_0[0][c_plus_d-1];
96  }
97 #if LIBINT2_FLOP_COUNT
98  inteval->nflops[0] += 4*(Lc+Ld-1);
99 #endif
100  }
101 
102  // build (1 0|0 0)
103  if (La+Lb > 0) {
104  apb_0_GTG_cpd_0[1][0] = pfac0_0[0] * apb_0_GTG_cpd_0[0][0];
105 #if LIBINT2_FLOP_COUNT
106  inteval->nflops[0] += 1;
107 #endif
108  }
109 
110  // build (a+b 0|0 0)
111  if (La+Lb > 1) {
112  for(int a_plus_b=1; a_plus_b!=La+Lb; ++a_plus_b) {
113  apb_0_GTG_cpd_0[a_plus_b+1][0] = pfac0_0[0] * apb_0_GTG_cpd_0[a_plus_b][0] +
114  a_plus_b * pfac1_0[0] * apb_0_GTG_cpd_0[a_plus_b-1][0];
115  }
116 #if LIBINT2_FLOP_COUNT
117  inteval->nflops[0] += 4*(La+Lb-1);
118 #endif
119  }
120 
121  // build (1 0|c+d 0)
122  if (La+Lb > 0 && Lc+Ld > 0) {
123  for(int c_plus_d=1; c_plus_d<=Lc+Ld; ++c_plus_d) {
124  apb_0_GTG_cpd_0[1][c_plus_d] = pfac0_0[0] * apb_0_GTG_cpd_0[0][c_plus_d] +
125  c_plus_d * pfac2[0] * apb_0_GTG_cpd_0[0][c_plus_d-1];
126  }
127 #if LIBINT2_FLOP_COUNT
128  inteval->nflops[0] += 4*(Lc+Ld-1);
129 #endif
130  }
131 
132  // build (a+b 0|c+d 0)
133  if (La+Lb > 1 && Lc+Ld > 0) {
134  for(int a_plus_b=1; a_plus_b!=La+Lb; ++a_plus_b) {
135  for(int c_plus_d=1; c_plus_d<=Lc+Ld; ++c_plus_d) {
136  apb_0_GTG_cpd_0[a_plus_b+1][c_plus_d] = pfac0_0[0] * apb_0_GTG_cpd_0[a_plus_b][c_plus_d] +
137  a_plus_b * pfac1_0[0] * apb_0_GTG_cpd_0[a_plus_b-1][c_plus_d] +
138  c_plus_d * pfac2[0] * apb_0_GTG_cpd_0[a_plus_b][c_plus_d-1];
139  }
140  }
141 #if LIBINT2_FLOP_COUNT
142  inteval->nflops[0] += 7*(La+Lb-1)*(Lc+Ld-1);
143 #endif
144  }
145 
146  //---------------------------------------------
147  // Part (2): build (a b|c+d 0)
148  //---------------------------------------------
149 
150  const LIBINT2_REALTYPE* AB;
151  switch (CartesianAxis) {
152  case x:
153  AB = inteval->AB_x;
154  break;
155  case y:
156  AB = inteval->AB_y;
157  break;
158  case z:
159  AB = inteval->AB_z;
160  break;
161  default: assert(false);
162  }
163 
164  LIBINT2_REALTYPE a_b_GTG_cpd_0[La+1][Lb+1][Lc+Ld+1];
165  for(int c_plus_d=0; c_plus_d<=Lc+Ld; ++c_plus_d) {
166  // copy (a+b 0| to a local 0,a+b buffer
167  LIBINT2_REALTYPE b_a_GTG[La+Lb+1][La+Lb+1];
168  for(int a_plus_b=0; a_plus_b<=La+Lb; ++a_plus_b) {
169  b_a_GTG[0][a_plus_b] = apb_0_GTG_cpd_0[a_plus_b][c_plus_d];
170  }
171  // use HRR to compute b,a
172  for(int b=1; b<=Lb; ++b) {
173  for(int a=0; a<=La+Lb-b; ++a) {
174  b_a_GTG[b][a] = b_a_GTG[b-1][a+1] + AB[0] * b_a_GTG[b-1][a];
175  }
176 #if LIBINT2_FLOP_COUNT
177  inteval->nflops[0] += 2 * (La+Lb-b+1);
178 #endif
179  }
180  // copy b,a to (a b|
181  for(int b=0; b<=Lb; ++b) {
182  for(int a=0; a<=La; ++a) {
183  a_b_GTG_cpd_0[a][b][c_plus_d] = b_a_GTG[b][a];
184  }
185  }
186  }
187 
188  //---------------------------------------------
189  // Part (3): build (a b|c d)
190  //---------------------------------------------
191 
192  const LIBINT2_REALTYPE* CD;
193  switch (CartesianAxis) {
194  case x:
195  CD = inteval->CD_x;
196  break;
197  case y:
198  CD = inteval->CD_y;
199  break;
200  case z:
201  CD = inteval->CD_z;
202  break;
203  default: assert(false);
204  }
205 
206  LIBINT2_REALTYPE* target_a_b_blk_ptr = target;
207  const int Nd = (Ld+1);
208  const int Ncd = (Lc+1)*Nd;
209  for(int a=0; a<=La; ++a) {
210  for(int b=0; b<=Lb; ++b, target_a_b_blk_ptr+=Ncd) {
211  // copy |c+d 0) to a local 0,c+d buffer
212  LIBINT2_REALTYPE d_c_GTG[Lc+Ld+1][Lc+Ld+1];
213  for(int c_plus_d=0; c_plus_d<=Lc+Ld; ++c_plus_d) {
214  d_c_GTG[0][c_plus_d] = a_b_GTG_cpd_0[a][b][c_plus_d];
215  }
216  // use HRR to compute d,c
217  for(int d=1; d<=Ld; ++d) {
218  for(int c=0; c<=Lc+Ld-d; ++c) {
219  d_c_GTG[d][c] = d_c_GTG[d-1][c+1] + CD[0] * d_c_GTG[d-1][c];
220  }
221 #if LIBINT2_FLOP_COUNT
222  inteval->nflops[0] += 2 * (Lc+Ld-d+1);
223 #endif
224  }
225  // copy d,c to |c d)
226  for(int d=0; d<=Ld; ++d) {
227  for(int c=0, cd=d; c<=Lc; ++c, cd+=Nd) {
228  target_a_b_blk_ptr[cd] = d_c_GTG[d][c];
229  }
230  }
231  }
232  }
233 
234  // done
235  }
236 
237  };
238 
239 };
240 
241 #endif // header guard
242 
Defaults definitions for various parameters assumed by Libint.
Definition: algebra.cc:23
Definition: test_eri_rys.cc:46
builds (ab| GTG_1d |cd), the shell set of 2-dimensional integrals needed for Rys quadrature evaluatio...
Definition: VRR_GTG_1d_xx_xx.h:37