1 /* Ergo, version 3.8, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2019 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU 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 * This program 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 General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /** @file integrals_2el_util_funcs.cc
31
32 \brief Code for utility functions used by 2-electron integral
33 computation (i.e. computation of J and K matrices).
34
35 @author: Elias Rudberg <em>responsible</em>.
36 */
37
38 #include "integrals_2el_util_funcs.h"
39
40 /* ELIAS NOTE 2014-05-29: The do_summedIntegralList_contribs_std
41 routine defined in this file is responsible for a large part of the
42 computational effort for both J and K matrix
43 construction. Therefore, this routine would be a good candidate for
44 further optimization attempts. */
45
46
47 /* This is the simple implementation, without unrolling. It turned
48 out that this could be optimized significantly by unrolling the
49 outer loop. */
50 /*
51 ELIAS NOTE 2014-07-13: Commented out this unused routine to silence compiler warning.
52 static void do_summedIntegralList_contribs_std_simple(const i_j_val_struct* conv_mat_1_sp, int conv_mat_1_sp_nnz,
53 const i_j_val_struct* conv_mat_2_sp, int conv_mat_2_sp_nnz,
54 int noOfMonomials_1, int noOfMonomials_2,
55 const ergo_real* primitiveIntegralList,
56 int noOfBasisFuncPairs_1, int noOfBasisFuncPairs_2,
57 ergo_real* summedIntegralList) {
58 for(int idx_i = 0; idx_i < conv_mat_1_sp_nnz; idx_i++) {
59 int idx_1 = conv_mat_1_sp[idx_i].i;
60 int ii = conv_mat_1_sp[idx_i].j;
61 ergo_real value_i = conv_mat_1_sp[idx_i].value;
62 const ergo_real* primitiveIntegralListPtr = &primitiveIntegralList[ii*noOfMonomials_2];
63 ergo_real* summedIntegralListPtr = &summedIntegralList[idx_1*noOfBasisFuncPairs_2];
64 int idx_j = 0;
65 while(idx_j < conv_mat_2_sp_nnz) {
66 int nn = conv_mat_2_sp[idx_j].same_i_count;
67 ergo_real sum = 0;
68 for(int kk = 0; kk < nn; kk++) {
69 int jj = conv_mat_2_sp[idx_j+kk].j;
70 sum += primitiveIntegralListPtr[jj] * conv_mat_2_sp[idx_j+kk].value;
71 }
72 int idx_2 = conv_mat_2_sp[idx_j].i;
73 summedIntegralListPtr[idx_2] += sum * value_i;
74 idx_j += nn;
75 }
76 }
77 }
78 */
79
80
do_summedIntegralList_contribs_std(const i_j_val_struct * conv_mat_1_sp,int conv_mat_1_sp_nnz,const i_j_val_struct * conv_mat_2_sp,int conv_mat_2_sp_nnz,int noOfMonomials_1,int noOfMonomials_2,const ergo_real * primitiveIntegralList,int noOfBasisFuncPairs_1,int noOfBasisFuncPairs_2,ergo_real * summedIntegralList)81 void do_summedIntegralList_contribs_std(const i_j_val_struct* conv_mat_1_sp, int conv_mat_1_sp_nnz,
82 const i_j_val_struct* conv_mat_2_sp, int conv_mat_2_sp_nnz,
83 int noOfMonomials_1, int noOfMonomials_2,
84 const ergo_real* primitiveIntegralList,
85 int noOfBasisFuncPairs_1, int noOfBasisFuncPairs_2,
86 ergo_real* summedIntegralList) {
87 int idx_i = 0;
88 while(idx_i < conv_mat_1_sp_nnz) {
89 if(idx_i == conv_mat_1_sp_nnz-1) {
90 // Only one left; treat in the old way
91 int idx_1 = conv_mat_1_sp[idx_i].i;
92 int ii = conv_mat_1_sp[idx_i].j;
93 ergo_real value_i = conv_mat_1_sp[idx_i].value;
94 const ergo_real* primitiveIntegralListPtr = &primitiveIntegralList[ii*noOfMonomials_2];
95 ergo_real* summedIntegralListPtr = &summedIntegralList[idx_1*noOfBasisFuncPairs_2];
96 int idx_j = 0;
97 while(idx_j < conv_mat_2_sp_nnz) {
98 int nn = conv_mat_2_sp[idx_j].same_i_count;
99 ergo_real sum = 0;
100 for(int kk = 0; kk < nn; kk++) {
101 int jj = conv_mat_2_sp[idx_j+kk].j;
102 sum += primitiveIntegralListPtr[jj] * conv_mat_2_sp[idx_j+kk].value;
103 }
104 int idx_2 = conv_mat_2_sp[idx_j].i;
105 summedIntegralListPtr[idx_2] += sum * value_i;
106 idx_j += nn;
107 }
108 idx_i++;
109 }
110 else if(idx_i == conv_mat_1_sp_nnz-2) {
111 // Unroll by 2
112 int idx_1_A = conv_mat_1_sp[idx_i].i;
113 int ii_A = conv_mat_1_sp[idx_i].j;
114 ergo_real value_i_A = conv_mat_1_sp[idx_i].value;
115 const ergo_real* primitiveIntegralListPtr_A = &primitiveIntegralList[ii_A*noOfMonomials_2];
116 ergo_real* summedIntegralListPtr_A = &summedIntegralList[idx_1_A*noOfBasisFuncPairs_2];
117 int idx_1_B = conv_mat_1_sp[idx_i+1].i;
118 int ii_B = conv_mat_1_sp[idx_i+1].j;
119 ergo_real value_i_B = conv_mat_1_sp[idx_i+1].value;
120 const ergo_real* primitiveIntegralListPtr_B = &primitiveIntegralList[ii_B*noOfMonomials_2];
121 ergo_real* summedIntegralListPtr_B = &summedIntegralList[idx_1_B*noOfBasisFuncPairs_2];
122 int idx_j = 0;
123 while(idx_j < conv_mat_2_sp_nnz) {
124 int nn = conv_mat_2_sp[idx_j].same_i_count;
125 ergo_real sum_A = 0;
126 ergo_real sum_B = 0;
127 for(int kk = 0; kk < nn; kk++) {
128 int jj = conv_mat_2_sp[idx_j+kk].j;
129 sum_A += primitiveIntegralListPtr_A[jj] * conv_mat_2_sp[idx_j+kk].value;
130 sum_B += primitiveIntegralListPtr_B[jj] * conv_mat_2_sp[idx_j+kk].value;
131 }
132 int idx_2 = conv_mat_2_sp[idx_j].i;
133 summedIntegralListPtr_A[idx_2] += sum_A * value_i_A;
134 summedIntegralListPtr_B[idx_2] += sum_B * value_i_B;
135 idx_j += nn;
136 }
137 idx_i += 2;
138 }
139 else if(idx_i == conv_mat_1_sp_nnz-3) {
140 // Unroll by 3
141 int idx_1_A = conv_mat_1_sp[idx_i].i;
142 int ii_A = conv_mat_1_sp[idx_i].j;
143 ergo_real value_i_A = conv_mat_1_sp[idx_i].value;
144 const ergo_real* primitiveIntegralListPtr_A = &primitiveIntegralList[ii_A*noOfMonomials_2];
145 ergo_real* summedIntegralListPtr_A = &summedIntegralList[idx_1_A*noOfBasisFuncPairs_2];
146 int idx_1_B = conv_mat_1_sp[idx_i+1].i;
147 int ii_B = conv_mat_1_sp[idx_i+1].j;
148 ergo_real value_i_B = conv_mat_1_sp[idx_i+1].value;
149 const ergo_real* primitiveIntegralListPtr_B = &primitiveIntegralList[ii_B*noOfMonomials_2];
150 ergo_real* summedIntegralListPtr_B = &summedIntegralList[idx_1_B*noOfBasisFuncPairs_2];
151 int idx_1_C = conv_mat_1_sp[idx_i+2].i;
152 int ii_C = conv_mat_1_sp[idx_i+2].j;
153 ergo_real value_i_C = conv_mat_1_sp[idx_i+2].value;
154 const ergo_real* primitiveIntegralListPtr_C = &primitiveIntegralList[ii_C*noOfMonomials_2];
155 ergo_real* summedIntegralListPtr_C = &summedIntegralList[idx_1_C*noOfBasisFuncPairs_2];
156 int idx_j = 0;
157 while(idx_j < conv_mat_2_sp_nnz) {
158 int nn = conv_mat_2_sp[idx_j].same_i_count;
159 ergo_real sum_A = 0;
160 ergo_real sum_B = 0;
161 ergo_real sum_C = 0;
162 for(int kk = 0; kk < nn; kk++) {
163 int jj = conv_mat_2_sp[idx_j+kk].j;
164 sum_A += primitiveIntegralListPtr_A[jj] * conv_mat_2_sp[idx_j+kk].value;
165 sum_B += primitiveIntegralListPtr_B[jj] * conv_mat_2_sp[idx_j+kk].value;
166 sum_C += primitiveIntegralListPtr_C[jj] * conv_mat_2_sp[idx_j+kk].value;
167 }
168 int idx_2 = conv_mat_2_sp[idx_j].i;
169 summedIntegralListPtr_A[idx_2] += sum_A * value_i_A;
170 summedIntegralListPtr_B[idx_2] += sum_B * value_i_B;
171 summedIntegralListPtr_C[idx_2] += sum_C * value_i_C;
172 idx_j += nn;
173 }
174 idx_i += 3;
175 }
176 else if(idx_i == conv_mat_1_sp_nnz-4) {
177 // Unroll by 4
178 int idx_1_A = conv_mat_1_sp[idx_i].i;
179 int ii_A = conv_mat_1_sp[idx_i].j;
180 ergo_real value_i_A = conv_mat_1_sp[idx_i].value;
181 const ergo_real* primitiveIntegralListPtr_A = &primitiveIntegralList[ii_A*noOfMonomials_2];
182 ergo_real* summedIntegralListPtr_A = &summedIntegralList[idx_1_A*noOfBasisFuncPairs_2];
183 int idx_1_B = conv_mat_1_sp[idx_i+1].i;
184 int ii_B = conv_mat_1_sp[idx_i+1].j;
185 ergo_real value_i_B = conv_mat_1_sp[idx_i+1].value;
186 const ergo_real* primitiveIntegralListPtr_B = &primitiveIntegralList[ii_B*noOfMonomials_2];
187 ergo_real* summedIntegralListPtr_B = &summedIntegralList[idx_1_B*noOfBasisFuncPairs_2];
188 int idx_1_C = conv_mat_1_sp[idx_i+2].i;
189 int ii_C = conv_mat_1_sp[idx_i+2].j;
190 ergo_real value_i_C = conv_mat_1_sp[idx_i+2].value;
191 const ergo_real* primitiveIntegralListPtr_C = &primitiveIntegralList[ii_C*noOfMonomials_2];
192 ergo_real* summedIntegralListPtr_C = &summedIntegralList[idx_1_C*noOfBasisFuncPairs_2];
193 int idx_1_D = conv_mat_1_sp[idx_i+3].i;
194 int ii_D = conv_mat_1_sp[idx_i+3].j;
195 ergo_real value_i_D = conv_mat_1_sp[idx_i+3].value;
196 const ergo_real* primitiveIntegralListPtr_D = &primitiveIntegralList[ii_D*noOfMonomials_2];
197 ergo_real* summedIntegralListPtr_D = &summedIntegralList[idx_1_D*noOfBasisFuncPairs_2];
198 int idx_j = 0;
199 while(idx_j < conv_mat_2_sp_nnz) {
200 int nn = conv_mat_2_sp[idx_j].same_i_count;
201 ergo_real sum_A = 0;
202 ergo_real sum_B = 0;
203 ergo_real sum_C = 0;
204 ergo_real sum_D = 0;
205 for(int kk = 0; kk < nn; kk++) {
206 int jj = conv_mat_2_sp[idx_j+kk].j;
207 sum_A += primitiveIntegralListPtr_A[jj] * conv_mat_2_sp[idx_j+kk].value;
208 sum_B += primitiveIntegralListPtr_B[jj] * conv_mat_2_sp[idx_j+kk].value;
209 sum_C += primitiveIntegralListPtr_C[jj] * conv_mat_2_sp[idx_j+kk].value;
210 sum_D += primitiveIntegralListPtr_D[jj] * conv_mat_2_sp[idx_j+kk].value;
211 }
212 int idx_2 = conv_mat_2_sp[idx_j].i;
213 summedIntegralListPtr_A[idx_2] += sum_A * value_i_A;
214 summedIntegralListPtr_B[idx_2] += sum_B * value_i_B;
215 summedIntegralListPtr_C[idx_2] += sum_C * value_i_C;
216 summedIntegralListPtr_D[idx_2] += sum_D * value_i_D;
217 idx_j += nn;
218 }
219 idx_i += 4;
220 }
221 else if(idx_i == conv_mat_1_sp_nnz-5) {
222 // Unroll by 5
223 int idx_1_A = conv_mat_1_sp[idx_i].i;
224 int ii_A = conv_mat_1_sp[idx_i].j;
225 ergo_real value_i_A = conv_mat_1_sp[idx_i].value;
226 const ergo_real* primitiveIntegralListPtr_A = &primitiveIntegralList[ii_A*noOfMonomials_2];
227 ergo_real* summedIntegralListPtr_A = &summedIntegralList[idx_1_A*noOfBasisFuncPairs_2];
228 int idx_1_B = conv_mat_1_sp[idx_i+1].i;
229 int ii_B = conv_mat_1_sp[idx_i+1].j;
230 ergo_real value_i_B = conv_mat_1_sp[idx_i+1].value;
231 const ergo_real* primitiveIntegralListPtr_B = &primitiveIntegralList[ii_B*noOfMonomials_2];
232 ergo_real* summedIntegralListPtr_B = &summedIntegralList[idx_1_B*noOfBasisFuncPairs_2];
233 int idx_1_C = conv_mat_1_sp[idx_i+2].i;
234 int ii_C = conv_mat_1_sp[idx_i+2].j;
235 ergo_real value_i_C = conv_mat_1_sp[idx_i+2].value;
236 const ergo_real* primitiveIntegralListPtr_C = &primitiveIntegralList[ii_C*noOfMonomials_2];
237 ergo_real* summedIntegralListPtr_C = &summedIntegralList[idx_1_C*noOfBasisFuncPairs_2];
238 int idx_1_D = conv_mat_1_sp[idx_i+3].i;
239 int ii_D = conv_mat_1_sp[idx_i+3].j;
240 ergo_real value_i_D = conv_mat_1_sp[idx_i+3].value;
241 const ergo_real* primitiveIntegralListPtr_D = &primitiveIntegralList[ii_D*noOfMonomials_2];
242 ergo_real* summedIntegralListPtr_D = &summedIntegralList[idx_1_D*noOfBasisFuncPairs_2];
243 int idx_1_E = conv_mat_1_sp[idx_i+4].i;
244 int ii_E = conv_mat_1_sp[idx_i+4].j;
245 ergo_real value_i_E = conv_mat_1_sp[idx_i+4].value;
246 const ergo_real* primitiveIntegralListPtr_E = &primitiveIntegralList[ii_E*noOfMonomials_2];
247 ergo_real* summedIntegralListPtr_E = &summedIntegralList[idx_1_E*noOfBasisFuncPairs_2];
248 int idx_j = 0;
249 while(idx_j < conv_mat_2_sp_nnz) {
250 int nn = conv_mat_2_sp[idx_j].same_i_count;
251 ergo_real sum_A = 0;
252 ergo_real sum_B = 0;
253 ergo_real sum_C = 0;
254 ergo_real sum_D = 0;
255 ergo_real sum_E = 0;
256 for(int kk = 0; kk < nn; kk++) {
257 int jj = conv_mat_2_sp[idx_j+kk].j;
258 sum_A += primitiveIntegralListPtr_A[jj] * conv_mat_2_sp[idx_j+kk].value;
259 sum_B += primitiveIntegralListPtr_B[jj] * conv_mat_2_sp[idx_j+kk].value;
260 sum_C += primitiveIntegralListPtr_C[jj] * conv_mat_2_sp[idx_j+kk].value;
261 sum_D += primitiveIntegralListPtr_D[jj] * conv_mat_2_sp[idx_j+kk].value;
262 sum_E += primitiveIntegralListPtr_E[jj] * conv_mat_2_sp[idx_j+kk].value;
263 }
264 int idx_2 = conv_mat_2_sp[idx_j].i;
265 summedIntegralListPtr_A[idx_2] += sum_A * value_i_A;
266 summedIntegralListPtr_B[idx_2] += sum_B * value_i_B;
267 summedIntegralListPtr_C[idx_2] += sum_C * value_i_C;
268 summedIntegralListPtr_D[idx_2] += sum_D * value_i_D;
269 summedIntegralListPtr_E[idx_2] += sum_E * value_i_E;
270 idx_j += nn;
271 }
272 idx_i += 5;
273 }
274 else {
275 // Unroll by 6
276 int idx_1_A = conv_mat_1_sp[idx_i].i;
277 int ii_A = conv_mat_1_sp[idx_i].j;
278 ergo_real value_i_A = conv_mat_1_sp[idx_i].value;
279 const ergo_real* primitiveIntegralListPtr_A = &primitiveIntegralList[ii_A*noOfMonomials_2];
280 ergo_real* summedIntegralListPtr_A = &summedIntegralList[idx_1_A*noOfBasisFuncPairs_2];
281 int idx_1_B = conv_mat_1_sp[idx_i+1].i;
282 int ii_B = conv_mat_1_sp[idx_i+1].j;
283 ergo_real value_i_B = conv_mat_1_sp[idx_i+1].value;
284 const ergo_real* primitiveIntegralListPtr_B = &primitiveIntegralList[ii_B*noOfMonomials_2];
285 ergo_real* summedIntegralListPtr_B = &summedIntegralList[idx_1_B*noOfBasisFuncPairs_2];
286 int idx_1_C = conv_mat_1_sp[idx_i+2].i;
287 int ii_C = conv_mat_1_sp[idx_i+2].j;
288 ergo_real value_i_C = conv_mat_1_sp[idx_i+2].value;
289 const ergo_real* primitiveIntegralListPtr_C = &primitiveIntegralList[ii_C*noOfMonomials_2];
290 ergo_real* summedIntegralListPtr_C = &summedIntegralList[idx_1_C*noOfBasisFuncPairs_2];
291 int idx_1_D = conv_mat_1_sp[idx_i+3].i;
292 int ii_D = conv_mat_1_sp[idx_i+3].j;
293 ergo_real value_i_D = conv_mat_1_sp[idx_i+3].value;
294 const ergo_real* primitiveIntegralListPtr_D = &primitiveIntegralList[ii_D*noOfMonomials_2];
295 ergo_real* summedIntegralListPtr_D = &summedIntegralList[idx_1_D*noOfBasisFuncPairs_2];
296 int idx_1_E = conv_mat_1_sp[idx_i+4].i;
297 int ii_E = conv_mat_1_sp[idx_i+4].j;
298 ergo_real value_i_E = conv_mat_1_sp[idx_i+4].value;
299 const ergo_real* primitiveIntegralListPtr_E = &primitiveIntegralList[ii_E*noOfMonomials_2];
300 ergo_real* summedIntegralListPtr_E = &summedIntegralList[idx_1_E*noOfBasisFuncPairs_2];
301 int idx_1_F = conv_mat_1_sp[idx_i+5].i;
302 int ii_F = conv_mat_1_sp[idx_i+5].j;
303 ergo_real value_i_F = conv_mat_1_sp[idx_i+5].value;
304 const ergo_real* primitiveIntegralListPtr_F = &primitiveIntegralList[ii_F*noOfMonomials_2];
305 ergo_real* summedIntegralListPtr_F = &summedIntegralList[idx_1_F*noOfBasisFuncPairs_2];
306 int idx_j = 0;
307 while(idx_j < conv_mat_2_sp_nnz) {
308 int nn = conv_mat_2_sp[idx_j].same_i_count;
309 ergo_real sum_A = 0;
310 ergo_real sum_B = 0;
311 ergo_real sum_C = 0;
312 ergo_real sum_D = 0;
313 ergo_real sum_E = 0;
314 ergo_real sum_F = 0;
315 for(int kk = 0; kk < nn; kk++) {
316 int jj = conv_mat_2_sp[idx_j+kk].j;
317 sum_A += primitiveIntegralListPtr_A[jj] * conv_mat_2_sp[idx_j+kk].value;
318 sum_B += primitiveIntegralListPtr_B[jj] * conv_mat_2_sp[idx_j+kk].value;
319 sum_C += primitiveIntegralListPtr_C[jj] * conv_mat_2_sp[idx_j+kk].value;
320 sum_D += primitiveIntegralListPtr_D[jj] * conv_mat_2_sp[idx_j+kk].value;
321 sum_E += primitiveIntegralListPtr_E[jj] * conv_mat_2_sp[idx_j+kk].value;
322 sum_F += primitiveIntegralListPtr_F[jj] * conv_mat_2_sp[idx_j+kk].value;
323 }
324 int idx_2 = conv_mat_2_sp[idx_j].i;
325 summedIntegralListPtr_A[idx_2] += sum_A * value_i_A;
326 summedIntegralListPtr_B[idx_2] += sum_B * value_i_B;
327 summedIntegralListPtr_C[idx_2] += sum_C * value_i_C;
328 summedIntegralListPtr_D[idx_2] += sum_D * value_i_D;
329 summedIntegralListPtr_E[idx_2] += sum_E * value_i_E;
330 summedIntegralListPtr_F[idx_2] += sum_F * value_i_F;
331 idx_j += nn;
332 }
333 idx_i += 6;
334 }
335 }
336 }
337
338
do_summedIntegralList_contribs_self(const i_j_val_struct * conv_mat_1_sp,int conv_mat_1_sp_nnz,const i_j_val_struct * conv_mat_2_sp,int conv_mat_2_sp_nnz,int noOfMonomials_1,int noOfMonomials_2,const ergo_real * primitiveIntegralList,int noOfBasisFuncPairs_1,int noOfBasisFuncPairs_2,ergo_real * summedIntegralList)339 void do_summedIntegralList_contribs_self(const i_j_val_struct* conv_mat_1_sp, int conv_mat_1_sp_nnz,
340 const i_j_val_struct* conv_mat_2_sp, int conv_mat_2_sp_nnz,
341 int noOfMonomials_1, int noOfMonomials_2,
342 const ergo_real* primitiveIntegralList,
343 int noOfBasisFuncPairs_1, int noOfBasisFuncPairs_2,
344 ergo_real* summedIntegralList) {
345 // Special interactionWithSelf case
346 for(int idx_i = 0; idx_i < conv_mat_1_sp_nnz; idx_i++) {
347 int idx_1 = conv_mat_1_sp[idx_i].i;
348 int ii = conv_mat_1_sp[idx_i].j;
349 ergo_real value_i = conv_mat_1_sp[idx_i].value;
350 const ergo_real* primitiveIntegralListPtr = &primitiveIntegralList[ii*noOfMonomials_2];
351 ergo_real* summedIntegralListPtr = &summedIntegralList[idx_1*noOfBasisFuncPairs_2];
352 int idx_j = 0;
353 while(idx_j < conv_mat_2_sp_nnz) {
354 int nn = conv_mat_2_sp[idx_j].same_i_count;
355 ergo_real sum = 0;
356 for(int kk = 0; kk < nn; kk++) {
357 int jj = conv_mat_2_sp[idx_j+kk].j;
358 sum += primitiveIntegralListPtr[jj] * conv_mat_2_sp[idx_j+kk].value;
359 }
360 int idx_2 = conv_mat_2_sp[idx_j].i;
361 if(idx_1 == idx_2)
362 summedIntegralListPtr[idx_2] += sum * value_i * 0.5;
363 else if(idx_1 > idx_2)
364 summedIntegralListPtr[idx_2] += sum * value_i;
365 idx_j += nn;
366 }
367 }
368 }
369