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