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 organize_distrs.cc
31
32 @brief Code for organizing a given set of primitive Gaussian
33 distributions (typically coming from basis function products); the
34 distributions are grouped according to their location in space,
35 their exponents, etc.
36
37 @author: Elias Rudberg <em>responsible</em>
38 */
39
40 #include <stdlib.h>
41 #include <memory.h>
42 #include <algorithm>
43
44 #include "organize_distrs.h"
45 #include "pi.h"
46 #include "serialization_tools.h"
47
48 #include <cstdio>
49
50 /* distr_org_struct functions */
51
Data()52 distr_org_struct::Data::Data():
53 maxExtent(0),
54 maxDistanceOutsideBox(0),
55 maxNoOfMonomials(0)
56 {}
57
writeToBuffer(char * dataBuffer,size_t const bufferSize) const58 void distr_org_struct::writeToBuffer(char* dataBuffer, size_t const bufferSize) const {
59 assert(bufferSize >= getSize());
60 char* p = dataBuffer;
61 memcpy(p, &data, sizeof(data));
62 p += sizeof(data);
63 std_vector_writeToBuffer_and_move_ptr(minimalDistrList, p);
64 std_vector_writeToBuffer_and_move_ptr(groupList, p);
65 std_vector_writeToBuffer_and_move_ptr(clusterList, p);
66 std_vector_writeToBuffer_and_move_ptr(batchList, p);
67 std_vector_writeToBuffer_and_move_ptr(basisFuncPairList, p);
68 std_vector_writeToBuffer_and_move_ptr(basisFuncListForBatchs, p);
69 std_vector_writeToBuffer_and_move_ptr(basisFuncListForBatchs_map, p);
70 std_vector_writeToBuffer_and_move_ptr(basisFuncList, p);
71 std_vector_writeToBuffer_and_move_ptr(spMatElementList, p);
72 std_vector_writeToBuffer_and_move_ptr(spMatCountList, p);
73 std_vector_writeToBuffer_and_move_ptr(spMatIdxList, p);
74 std_vector_writeToBuffer_and_move_ptr(basisFuncGroupInfoListForK, p);
75 }
76
getSize() const77 size_t distr_org_struct::getSize() const {
78 size_t size = sizeof(distr_org_struct::Data);
79 size += std_vector_getSize(minimalDistrList);
80 size += std_vector_getSize(groupList);
81 size += std_vector_getSize(clusterList);
82 size += std_vector_getSize(batchList);
83 size += std_vector_getSize(basisFuncPairList);
84 size += std_vector_getSize(basisFuncListForBatchs);
85 size += std_vector_getSize(basisFuncListForBatchs_map);
86 size += std_vector_getSize(basisFuncList);
87 size += std_vector_getSize(spMatElementList);
88 size += std_vector_getSize(spMatCountList);
89 size += std_vector_getSize(spMatIdxList);
90 size += std_vector_getSize(basisFuncGroupInfoListForK);
91 return size;
92 }
93
assignFromBuffer(char const * dataBuffer,size_t const bufferSize)94 void distr_org_struct::assignFromBuffer(char const * dataBuffer, size_t const bufferSize) {
95 const char* p = dataBuffer;
96 size_t remainingBytes = bufferSize;
97 assert(remainingBytes >= sizeof(data));
98 memcpy(&data, p, sizeof(data));
99 p += sizeof(data);
100 const char* bufEndPtr = &dataBuffer[bufferSize];
101 std_vector_assignFromBuffer_and_move_ptr(minimalDistrList, p, bufEndPtr);
102 std_vector_assignFromBuffer_and_move_ptr(groupList, p, bufEndPtr);
103 std_vector_assignFromBuffer_and_move_ptr(clusterList, p, bufEndPtr);
104 std_vector_assignFromBuffer_and_move_ptr(batchList, p, bufEndPtr);
105 std_vector_assignFromBuffer_and_move_ptr(basisFuncPairList, p, bufEndPtr);
106 std_vector_assignFromBuffer_and_move_ptr(basisFuncListForBatchs, p, bufEndPtr);
107 std_vector_assignFromBuffer_and_move_ptr(basisFuncListForBatchs_map, p, bufEndPtr);
108 std_vector_assignFromBuffer_and_move_ptr(basisFuncList, p, bufEndPtr);
109 std_vector_assignFromBuffer_and_move_ptr(spMatElementList, p, bufEndPtr);
110 std_vector_assignFromBuffer_and_move_ptr(spMatCountList, p, bufEndPtr);
111 std_vector_assignFromBuffer_and_move_ptr(spMatIdxList, p, bufEndPtr);
112 std_vector_assignFromBuffer_and_move_ptr(basisFuncGroupInfoListForK, p, bufEndPtr);
113 }
114
115 /* ****************************** */
116
117 static void
do_sort_int_list(int * list,int n)118 do_sort_int_list(int* list, int n)
119 {
120 for(int i = 0; i < n; i++)
121 for(int j = 0; j < n-i-1; j++)
122 {
123 if(list[j+1] < list[j])
124 {
125 int temp = list[j];
126 list[j] = list[j+1];
127 list[j+1] = temp;
128 }
129 } // END FOR i j
130 }
131
132
133
get_conversion_matrix_for_group(const IntegralInfo & integralInfo,const distr_group_struct & group,int n1max,const minimal_distr_struct * minimalDistrList_1,int noOfBasisFuncPairs_1,const i_j_val_struct * convMat1_sp,int convMat1_nnz,i_j_val_struct * BB1_x_Ai1_x_convMat1_sp_result,int & BB1_x_Ai1_x_convMat1_nnz_result)134 static void get_conversion_matrix_for_group(
135 const IntegralInfo & integralInfo,
136 const distr_group_struct & group,
137 int n1max,
138 const minimal_distr_struct* minimalDistrList_1,
139 int noOfBasisFuncPairs_1,
140 const i_j_val_struct* convMat1_sp,
141 int convMat1_nnz,
142 i_j_val_struct* BB1_x_Ai1_x_convMat1_sp_result, // result
143 int & BB1_x_Ai1_x_convMat1_nnz_result) // result
144 {
145 int noOfMonomials_1 = integralInfo.monomial_info.no_of_monomials_list[n1max];
146 int distrCount_i = group.distrCount;
147 i_j_val_struct Ai1_sp[distrCount_i];
148 for(int i = 0; i < distrCount_i; i++) {
149 int monomialIndex = minimalDistrList_1[group.startIndex+i].monomialIndex;
150 ergo_real value = minimalDistrList_1[group.startIndex+i].coeff;
151 Ai1_sp[i].i = i;
152 Ai1_sp[i].j = monomialIndex;
153 Ai1_sp[i].value = value;
154 }
155 i_j_val_struct BB1_sp[distrCount_i];
156 for(int kk = 0; kk < distrCount_i; kk++) {
157 int idx = minimalDistrList_1[kk+group.startIndex].basisFuncPairIndex;
158 BB1_sp[kk].i = idx;
159 BB1_sp[kk].j = kk;
160 BB1_sp[kk].value = 1;
161 }
162 // Multiply Ai1 by convMat1. Dimensions: (distrCount_i*noOfMonomials_1) x (noOfMonomials_1*noOfMonomials_1)
163 i_j_val_struct* Ai1_x_convMat1_sp = new i_j_val_struct[distrCount_i*noOfMonomials_1];
164 int Ai1_x_convMat1_nnz = spmat_multiply_matrices(Ai1_sp, distrCount_i, convMat1_sp, convMat1_nnz, Ai1_x_convMat1_sp, distrCount_i, noOfMonomials_1);
165 // Multiply BB1 by Ai1_x_convMat1. Dimensions: (noOfBasisFuncPairs_1*distrCount_i) x (distrCount_i*noOfMonomials_1)
166 BB1_x_Ai1_x_convMat1_nnz_result = spmat_multiply_matrices(BB1_sp, distrCount_i, Ai1_x_convMat1_sp, Ai1_x_convMat1_nnz, BB1_x_Ai1_x_convMat1_sp_result, noOfBasisFuncPairs_1, noOfMonomials_1);
167 delete [] Ai1_x_convMat1_sp;
168 }
169
170
171 template <typename T>
copy_vector(std::vector<T> & dest,std::vector<T> & src,int count)172 void copy_vector(std::vector<T> & dest, std::vector<T> & src, int count) {
173 // Do memcpy only if count > 0 to avoid problem if -D_GLIBCXX_ASSERTIONS is used
174 if(count > 0)
175 memcpy(&dest[0], &src[0], count*sizeof(T));
176 }
177
178
179 int
organize_distributions(const IntegralInfo & integralInfo,DistributionSpecStructLabeled * distrList_in,int distrCount,distr_org_struct * result,const ergo_real * boxCenterCoords,ergo_real boxWidth)180 organize_distributions(const IntegralInfo & integralInfo,
181 DistributionSpecStructLabeled* distrList_in,
182 int distrCount,
183 distr_org_struct* result,
184 const ergo_real* boxCenterCoords,
185 ergo_real boxWidth)
186 {
187 std::vector<DistributionSpecStructLabeled> distrList(distrCount);
188
189 // sort list of distributions by center, type and exponent
190 // first group the ones that have same center and same exponent.
191 std::vector<int> groupCountList(distrCount);
192 std::vector<int> groupIndexList(distrCount);
193
194
195 // start by bucket sort based on "best" coordinate.
196 const ergo_real HUGE_NUMBER = 888888888;
197 ergo_real xminList[3];
198 ergo_real xmaxList[3];
199 ergo_real xdiffList[3];
200 for(int kk = 0; kk < 3; kk++)
201 {
202 xminList[kk] = HUGE_NUMBER;
203 xmaxList[kk] = -HUGE_NUMBER;
204 }
205 for(int i = 0; i < distrCount; i++)
206 {
207 for(int kk = 0; kk < 3; kk++)
208 {
209 ergo_real x = distrList_in[i].distr.centerCoords[kk];
210 if(x < xminList[kk])
211 xminList[kk] = x;
212 if(x > xmaxList[kk])
213 xmaxList[kk] = x;
214 }
215 } // END FOR i
216 int bestCoordIndex = 0;
217 for(int kk = 0; kk < 3; kk++)
218 {
219 xdiffList[kk] = xmaxList[kk] - xminList[kk];
220 if(xdiffList[kk] > xdiffList[bestCoordIndex])
221 bestCoordIndex = kk;
222 }
223 #define NO_OF_SORT_BUCKETS 30
224 ergo_real splitterList[NO_OF_SORT_BUCKETS-1];
225 for(int i = 0; i < NO_OF_SORT_BUCKETS-1; i++)
226 splitterList[i] = xminList[bestCoordIndex] + ((ergo_real)i + 1) * xdiffList[bestCoordIndex] / NO_OF_SORT_BUCKETS;
227 int* bucketList[NO_OF_SORT_BUCKETS];
228 int bucketCounterList[NO_OF_SORT_BUCKETS];
229 for(int i = 0; i < NO_OF_SORT_BUCKETS; i++)
230 {
231 bucketList[i] = new int[distrCount];
232 bucketCounterList[i] = 0;
233 }
234 for(int i = 0; i < distrCount; i++)
235 {
236 int bucketIndex = -1;
237 for(int j = 0; j < NO_OF_SORT_BUCKETS-1; j++)
238 {
239 if(distrList_in[i].distr.centerCoords[bestCoordIndex] < splitterList[j])
240 {
241 bucketIndex = j;
242 break;
243 }
244 }
245 if(bucketIndex == -1)
246 bucketIndex = NO_OF_SORT_BUCKETS-1;
247 bucketList[bucketIndex][bucketCounterList[bucketIndex]] = i;
248 bucketCounterList[bucketIndex]++;
249 } // END FOR i
250
251 int destCount = 0;
252 int groupCount = 0;
253
254 // create groups for one bucket at a time
255 for(int bucketIndex = 0; bucketIndex < NO_OF_SORT_BUCKETS; bucketIndex++)
256 {
257 int nLeft = bucketCounterList[bucketIndex];
258 while(nLeft > 0)
259 {
260 int i = 0;
261 int remainingIndex = 0;
262 int destCountSaved = destCount;
263 distrList[destCount] = distrList_in[bucketList[bucketIndex][i]];
264 destCount++;
265 // now find all that belong to same group
266 for(int k = i+1; k < nLeft; k++)
267 {
268 ergo_real dx, dy, dz;
269 dx = distrList_in[bucketList[bucketIndex][k]].distr.centerCoords[0] - distrList[destCountSaved].distr.centerCoords[0];
270 dy = distrList_in[bucketList[bucketIndex][k]].distr.centerCoords[1] - distrList[destCountSaved].distr.centerCoords[1];
271 dz = distrList_in[bucketList[bucketIndex][k]].distr.centerCoords[2] - distrList[destCountSaved].distr.centerCoords[2];
272 ergo_real r2 = dx*dx + dy*dy + dz*dz;
273 ergo_real absExponentDiff = distrList_in[bucketList[bucketIndex][k]].distr.exponent - distrList[destCountSaved].distr.exponent;
274 if(absExponentDiff < 0)
275 absExponentDiff *= -1;
276 if(absExponentDiff < 1e-11 && r2 < 1e-10)
277 {
278 // OK, close enough, we regard this as being same center and same exponent.
279 // add to distrList, and remove from distrList_in.
280 distrList[destCount] = distrList_in[bucketList[bucketIndex][k]];
281 destCount++;
282 }
283 else
284 {
285 // no, different center or exponent
286 if(remainingIndex != k)
287 bucketList[bucketIndex][remainingIndex] = bucketList[bucketIndex][k];
288 remainingIndex++;
289 }
290 } // END FOR k find all that belong to same group
291 int noOfDistrsInGroup = destCount - destCountSaved;
292 nLeft -= noOfDistrsInGroup;
293 groupCountList[groupCount] = noOfDistrsInGroup;
294 groupCount++;
295 if(remainingIndex == 0)
296 break;
297 } // END WHILE group the ones that have same center and same exponent.
298 } // END FOR bucketIndex
299
300 for(int i = 0; i < NO_OF_SORT_BUCKETS; i++)
301 {
302 delete [] bucketList[i];
303 bucketList[i] = NULL;
304 }
305
306 // set groupIndexList
307 int currGroupIndex = 0;
308 for(int i = 0; i < groupCount; i++)
309 {
310 groupIndexList[i] = currGroupIndex;
311 currGroupIndex += groupCountList[i];
312 }
313
314 // Set groupID
315 for(int i = 0; i < groupCount; i++)
316 {
317 DistributionSpecStructLabeled* groupPtr = &distrList[groupIndexList[i]];
318 int currCount = groupCountList[i];
319 for(int j = 0; j < currCount; j++)
320 groupPtr[j].groupID = i + 1;
321 } // END FOR i
322
323 // Within each group, sort by monomialInts and basisFuncIndeces
324 for(int i = 0; i < groupCount; i++)
325 {
326 DistributionSpecStructLabeled* groupPtr = &distrList[groupIndexList[i]];
327 int currCount = groupCountList[i];
328 for(int k = 0; k < currCount; k++)
329 for(int m = 0; m < currCount - 1 - k; m++)
330 {
331 int doSwitch = 0;
332 if(doSwitch == 0 && groupPtr[m].distr.monomialInts[0] > groupPtr[m+1].distr.monomialInts[0])
333 doSwitch = 1;
334 else
335 doSwitch = -1;
336 if(doSwitch == 0 && groupPtr[m].distr.monomialInts[1] > groupPtr[m+1].distr.monomialInts[1])
337 doSwitch = 1;
338 else
339 doSwitch = -1;
340 if(doSwitch == 0 && groupPtr[m].distr.monomialInts[2] > groupPtr[m+1].distr.monomialInts[2])
341 doSwitch = 1;
342 else
343 doSwitch = -1;
344 if(doSwitch == 0 && groupPtr[m].basisFuncIndex_1 > groupPtr[m+1].basisFuncIndex_1)
345 doSwitch = 1;
346 else
347 doSwitch = -1;
348 if(doSwitch == 0 && groupPtr[m].basisFuncIndex_2 > groupPtr[m+1].basisFuncIndex_2)
349 doSwitch = 1;
350 else
351 doSwitch = -1;
352 if(doSwitch == 1)
353 {
354 // switch
355 DistributionSpecStructLabeled temp;
356 temp = groupPtr[m];
357 groupPtr[m] = groupPtr[m+1];
358 groupPtr[m+1] = temp;
359 }
360 } // END FOR k m
361 } // END FOR i
362
363
364 result->groupList.resize(groupCount);
365 std::vector<distr_group_struct> & groupList = result->groupList;
366
367 for(int i = 0; i < groupCount; i++)
368 {
369 groupList[i].distrCount = groupCountList[i];
370 groupList[i].startIndex = groupIndexList[i];
371 // get nmax
372 int nmax = 0;
373 for(int ii = groupIndexList[i]; ii < groupIndexList[i] + groupCountList[i]; ii++)
374 {
375 int sum = 0;
376 for(int kk = 0; kk < 3; kk++)
377 sum += distrList[ii].distr.monomialInts[kk];
378 if(sum > nmax)
379 nmax = sum;
380 }
381 groupList[i].nmax = nmax;
382 // get centerCoords and exponent
383 for(int ii = 0; ii < 3; ii++)
384 groupList[i].centerCoords[ii] = distrList[groupIndexList[i]].distr.centerCoords[ii];
385 groupList[i].exponent = distrList[groupIndexList[i]].distr.exponent;
386 // get maxSize, maxLimitingFactor, maxExtent for this group.
387 ergo_real maxSize = 0;
388 ergo_real maxLimitingFactor = 0;
389 ergo_real maxExtent = 0;
390 for(int ii = groupIndexList[i]; ii < groupIndexList[i] + groupCountList[i]; ii++)
391 {
392 ergo_real size = template_blas_fabs(template_blas_pow((ergo_real)pi/distrList[ii].distr.exponent, (ergo_real)1.5) * distrList[ii].distr.coeff);
393 if(size > maxSize)
394 maxSize = size;
395 ergo_real limitingFactor = distrList[ii].limitingFactor;
396 if(limitingFactor > maxLimitingFactor)
397 maxLimitingFactor = limitingFactor;
398 ergo_real extent = distrList[ii].distr.extent;
399 if(extent > maxExtent)
400 maxExtent = extent;
401 }
402 groupList[i].maxSizeGroup = maxSize;
403 groupList[i].maxLimitingFactorGroup = maxLimitingFactor;
404 groupList[i].maxExtentGroup = maxExtent;
405
406 // Get maxAbsDmatElementGroup
407 ergo_real maxabs = 0;
408 for(int ii = groupIndexList[i]; ii < groupIndexList[i] + groupCountList[i]; ii++)
409 {
410 ergo_real absval = template_blas_fabs(distrList[ii].dmatElement);
411 if(absval > maxabs)
412 maxabs = absval;
413 }
414 groupList[i].maxAbsDmatElementGroup = maxabs;
415 } // END FOR i
416
417 #define MAX_NO_OF_GROUPS_PER_CLUSTER 10
418
419
420 // create clusters and batchs.
421 // move groups into new list, one cluster at a time.
422 int batchCount = 0;
423 int clusterCount = 0;
424 int basisFuncPairCount = 0;
425
426 std::vector<distr_group_struct> groupList2(groupCount);
427
428 std::vector<cluster_struct> clusterList(groupCount);
429 std::vector<batch_struct> batchList(groupCount);
430
431 std::vector<basis_func_pair_struct> basisFuncPairList(distrCount);
432
433 int noOfGroupsInNewList = 0;
434 int noOfGroupsLeftInOldList = groupCount;
435 while(noOfGroupsInNewList < groupCount)
436 {
437 // the group that is first now will define the beginning of a new cluster, and a new batch.
438 batch_struct newBatch;
439 memset(&newBatch, 0, sizeof(batch_struct));
440 clusterList[clusterCount].groupStartIndex = noOfGroupsInNewList;
441 newBatch.clusterStartIndex = clusterCount;
442 newBatch.basisFuncPairListIndex = basisFuncPairCount;
443
444 // add basisFuncPairs for first group to newBatch
445 for(int i = groupList[0].startIndex; i < groupList[0].startIndex + groupList[0].distrCount; i++)
446 {
447 int alreadyInList = 0;
448 for(int kk = 0; kk < newBatch.noOfBasisFuncPairs; kk++)
449 {
450 if(distrList[i].basisFuncIndex_1 == basisFuncPairList[newBatch.basisFuncPairListIndex+kk].index_1 &&
451 distrList[i].basisFuncIndex_2 == basisFuncPairList[newBatch.basisFuncPairListIndex+kk].index_2)
452 {
453 alreadyInList = 1;
454 break;
455 }
456 } // END FOR kk
457 if(alreadyInList == 0)
458 {
459 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].index_1 = distrList[i].basisFuncIndex_1;
460 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].index_2 = distrList[i].basisFuncIndex_2;
461 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].pairIndex = distrList[i].pairIndex;
462 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].dmatElement = distrList[i].dmatElement;
463 newBatch.noOfBasisFuncPairs++;
464 basisFuncPairCount++;
465 if(newBatch.noOfBasisFuncPairs >= MAX_NO_OF_BASIS_FUNC_PAIRS_PER_BATCH)
466 {
467 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: (newBatch.noOfBasisFuncPairs >= MAX_NO_OF_BASIS_FUNC_PAIRS_PER_BATCH)");
468 return -1;
469 }
470 }
471 } // END FOR i add basisFuncPairs for first group to newBatch
472
473 int noOfClustersInCurrBatch = 1;
474 int oldListIndex = 0;
475 memcpy(&groupList2[noOfGroupsInNewList], &groupList[0], sizeof(distr_group_struct));
476 noOfGroupsInNewList++;
477 int noOfGroupsInCurrCluster = 1;
478 // now find other groups with same exponent and same nmax
479 ergo_real exponent = groupList[0].exponent;
480 int nmax = groupList[0].nmax;
481 for(int i = 1; i < noOfGroupsLeftInOldList; i++)
482 {
483 ergo_real absexponentDiff = template_blas_fabs(exponent - groupList[i].exponent);
484 if(absexponentDiff < 1e-11 && groupList[i].nmax == nmax && noOfGroupsInCurrCluster < MAX_NO_OF_GROUPS_PER_CLUSTER)
485 {
486 // same exponent and nmax found, add this group to cluster
487 memcpy(&groupList2[noOfGroupsInNewList], &groupList[i], sizeof(distr_group_struct));
488 noOfGroupsInNewList++;
489 noOfGroupsInCurrCluster++;
490 // add basisFuncPairs for group to newBatch
491 for(int ii = groupList[i].startIndex; ii < groupList[i].startIndex + groupList[i].distrCount; ii++)
492 {
493 int alreadyInList = 0;
494 for(int kk = 0; kk < newBatch.noOfBasisFuncPairs; kk++)
495 {
496 if(distrList[ii].basisFuncIndex_1 == basisFuncPairList[newBatch.basisFuncPairListIndex+kk].index_1 &&
497 distrList[ii].basisFuncIndex_2 == basisFuncPairList[newBatch.basisFuncPairListIndex+kk].index_2)
498 {
499 alreadyInList = 1;
500 break;
501 }
502 }
503 if(alreadyInList == 0)
504 {
505 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].index_1 = distrList[ii].basisFuncIndex_1;
506 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].index_2 = distrList[ii].basisFuncIndex_2;
507 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].pairIndex = distrList[ii].pairIndex;
508 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].dmatElement = distrList[ii].dmatElement;
509 newBatch.noOfBasisFuncPairs++;
510 basisFuncPairCount++;
511 if(newBatch.noOfBasisFuncPairs >= MAX_NO_OF_BASIS_FUNC_PAIRS_PER_BATCH)
512 {
513 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: (newBatch.noOfBasisFuncPairs >= MAX_NO_OF_BASIS_FUNC_PAIRS_PER_BATCH)");
514 return -1;
515 }
516 }
517 }
518 }
519 else
520 {
521 memcpy(&groupList[oldListIndex], &groupList[i], sizeof(distr_group_struct));
522 oldListIndex++;
523 }
524 } // END FOR i
525 noOfGroupsLeftInOldList -= noOfGroupsInCurrCluster;
526 clusterList[clusterCount].noOfGroups = noOfGroupsInCurrCluster;
527 clusterCount++;
528 // the cluster just created is the first one in a new batch.
529 // if possible, we want to add more clusters for that batch.
530 int definingClusterStartGrIndex = clusterList[clusterCount-1].groupStartIndex;
531 int definingClusterGrCount = clusterList[clusterCount-1].noOfGroups;
532 // look for other clusters to put in the same batch.
533 noOfGroupsInCurrCluster = 0;
534 while(noOfGroupsInNewList < groupCount)// && noOfGroupsInCurrCluster < MAX_NO_OF_GROUPS_PER_CLUSTER)
535 {
536 // look for a group that has the right basis funcs.
537 int foundIndex = -1;
538 for(int i = 0; i < noOfGroupsLeftInOldList; i++)
539 {
540 // we demand that all basisfuncpairs must be present in the batch (defined by first cluster)
541 int allPresentSoFar = 1;
542 for(int ii = 0; ii < groupList[i].distrCount; ii++)
543 {
544 // check if this distr is present in the batch
545 int bfidx1 = distrList[groupList[i].startIndex+ii].basisFuncIndex_1;
546 int bfidx2 = distrList[groupList[i].startIndex+ii].basisFuncIndex_2;
547 int found = 0;
548 for(int gr = definingClusterStartGrIndex; gr < definingClusterStartGrIndex + definingClusterGrCount; gr++)
549 {
550 int idistr;
551 for(idistr = 0; idistr < groupList2[gr].distrCount; idistr++)
552 {
553 if(distrList[groupList2[gr].startIndex+idistr].basisFuncIndex_1 == bfidx1 && distrList[groupList2[gr].startIndex+idistr].basisFuncIndex_2 == bfidx2)
554 {
555 found = 1;
556 break;
557 }
558 }
559 if(found == 1)
560 break;
561 }
562 if(found == 0)
563 {
564 allPresentSoFar = 0;
565 break;
566 }
567 } // END FOR ii
568 if(allPresentSoFar == 1)
569 {
570 // OK, use this group
571 foundIndex = i;
572 break;
573 }
574 } // END FOR i look for a group that has the right basis funcs.
575 if(foundIndex == -1)
576 break;
577 // OK, we have a group with accepted basis funcs.
578 // This group will be the first in a new cluster.
579
580 clusterList[clusterCount].groupStartIndex = noOfGroupsInNewList;
581 int oldListIndex = 0;
582 memcpy(&groupList2[noOfGroupsInNewList], &groupList[foundIndex], sizeof(distr_group_struct));
583 noOfGroupsInNewList++;
584 noOfGroupsInCurrCluster = 1;
585
586 // add basisFuncPairs for group to newBatch
587 for(int ii = groupList[foundIndex].startIndex; ii < groupList[foundIndex].startIndex + groupList[foundIndex].distrCount; ii++)
588 {
589 int alreadyInList = 0;
590 for(int kk = 0; kk < newBatch.noOfBasisFuncPairs; kk++)
591 {
592 if(distrList[ii].basisFuncIndex_1 == basisFuncPairList[newBatch.basisFuncPairListIndex+kk].index_1 &&
593 distrList[ii].basisFuncIndex_2 == basisFuncPairList[newBatch.basisFuncPairListIndex+kk].index_2)
594 {
595 alreadyInList = 1;
596 break;
597 }
598 }
599 if(alreadyInList == 0)
600 {
601 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].index_1 = distrList[ii].basisFuncIndex_1;
602 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].index_2 = distrList[ii].basisFuncIndex_2;
603 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].pairIndex = distrList[ii].pairIndex;
604 basisFuncPairList[newBatch.basisFuncPairListIndex+newBatch.noOfBasisFuncPairs].dmatElement = distrList[ii].dmatElement;
605 newBatch.noOfBasisFuncPairs++;
606 basisFuncPairCount++;
607 if(newBatch.noOfBasisFuncPairs >= MAX_NO_OF_BASIS_FUNC_PAIRS_PER_BATCH)
608 {
609 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: (newBatch.noOfBasisFuncPairs >= MAX_NO_OF_BASIS_FUNC_PAIRS_PER_BATCH)");
610 return -1;
611 }
612 }
613 }
614
615 ergo_real exponent = groupList[foundIndex].exponent;
616 int nmax = groupList[foundIndex].nmax;
617
618 // we have copied the entry at foundIndex to new list, all after that must be moved one step.
619 for(int i = foundIndex+1; i < noOfGroupsLeftInOldList; i++)
620 memcpy(&groupList[i-1], &groupList[i], sizeof(distr_group_struct));
621 noOfGroupsLeftInOldList--;
622
623 int noOfGroupsInCurrCluster = 1;
624 // now find other groups with same exponent and same nmax and accepted basis funcs
625 oldListIndex = 0;
626 for(int i = 0; i < noOfGroupsLeftInOldList; i++)
627 {
628 int addToCluster = 0;
629 ergo_real absexponentDiff = template_blas_fabs(exponent - groupList[i].exponent);
630 if(absexponentDiff < 1e-11 && groupList[i].nmax == nmax && noOfGroupsInCurrCluster < MAX_NO_OF_GROUPS_PER_CLUSTER)
631 {
632 // same exponent and nmax found, now check basis funcs
633 int allPresentSoFar = 1;
634 for(int ii = 0; ii < groupList[i].distrCount; ii++)
635 {
636 // check if this distr is present in the batch
637 int bfidx1 = distrList[groupList[i].startIndex+ii].basisFuncIndex_1;
638 int bfidx2 = distrList[groupList[i].startIndex+ii].basisFuncIndex_2;
639 int found = 0;
640 for(int gr = definingClusterStartGrIndex; gr < definingClusterStartGrIndex + definingClusterGrCount; gr++)
641 {
642 for(int idistr = 0; idistr < groupList2[gr].distrCount; idistr++)
643 {
644 if(distrList[groupList2[gr].startIndex+idistr].basisFuncIndex_1 == bfidx1 &&
645 distrList[groupList2[gr].startIndex+idistr].basisFuncIndex_2 == bfidx2)
646 {
647 found = 1;
648 break;
649 }
650 }
651 if(found == 1)
652 break;
653 }
654 if(found == 0)
655 {
656 allPresentSoFar = 0;
657 break;
658 }
659 } // END FOR ii
660 if(allPresentSoFar == 1)
661 addToCluster = 1;
662 }
663 if(addToCluster == 1)
664 {
665 // same exponent and nmax found and accepted funcs, add this group to cluster
666 memcpy(&groupList2[noOfGroupsInNewList], &groupList[i], sizeof(distr_group_struct));
667 noOfGroupsInNewList++;
668 noOfGroupsInCurrCluster++;
669 }
670 else
671 {
672 if(i != oldListIndex)
673 memcpy(&groupList[oldListIndex], &groupList[i], sizeof(distr_group_struct));
674 oldListIndex++;
675 }
676 } // END FOR i
677 noOfGroupsLeftInOldList -= noOfGroupsInCurrCluster-1;
678 clusterList[clusterCount].noOfGroups = noOfGroupsInCurrCluster;
679 clusterCount++;
680 noOfClustersInCurrBatch++;
681 } // END WHILE look for other clusters to put in the same batch
682
683 newBatch.noOfClusters = noOfClustersInCurrBatch;
684 batchList[batchCount] = newBatch;
685 batchCount++;
686
687 } // END WHILE create clusters
688
689 // check all batchs
690 for(int i = 0; i < batchCount; i++)
691 {
692 for(int j = 0; j < batchList[i].noOfBasisFuncPairs; j++)
693 {
694 for(int k = 0; k < batchList[i].noOfBasisFuncPairs; k++)
695 {
696 if(j != k &&
697 basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_1 == basisFuncPairList[batchList[i].basisFuncPairListIndex+k].index_1 &&
698 basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_2 == basisFuncPairList[batchList[i].basisFuncPairListIndex+k].index_2)
699 {
700 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: basisFuncPairs not unique in batch");
701 return -1;
702 }
703 }
704 }
705 }
706
707
708
709 copy_vector<distr_group_struct>(groupList, groupList2, groupCount);
710
711
712 // OK, clusters and batchs done.
713
714
715
716 // set nmax and exponent for all clusters
717 for(int i = 0; i < clusterCount; i++)
718 {
719 int groupStartIndex = clusterList[i].groupStartIndex;
720 int nGroups = clusterList[i].noOfGroups;
721 int nmax = 0;
722 ergo_real exponent = groupList[groupStartIndex].exponent;
723 for(int j = groupStartIndex; j < groupStartIndex + nGroups; j++)
724 {
725 if(groupList[j].nmax > nmax)
726 nmax = groupList[j].nmax;
727 ergo_real exponentdiff = template_blas_fabs(groupList[j].exponent - exponent);
728 if(exponentdiff > 1e-11)
729 {
730 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: different exponents found in same cluster");
731 return -1;
732 }
733 } // END FOR j
734 clusterList[i].nmax = nmax;
735 clusterList[i].exponent = exponent;
736 } // END FOR i set nmax for all clusters
737
738
739
740
741 // Sort clusters according to batchs
742 std::vector<cluster_struct> tempClusterList(clusterCount);
743 int count = 0;
744 for(int i = 0; i < batchCount; i++)
745 {
746 int savedCount = count;
747 for(int j = batchList[i].clusterStartIndex; j < batchList[i].clusterStartIndex + batchList[i].noOfClusters; j++)
748 {
749 tempClusterList[count] = clusterList[j];
750 count++;
751 } // END FOR j
752 batchList[i].clusterStartIndex = savedCount;
753 } // END FOR i
754 copy_vector<cluster_struct>(clusterList, tempClusterList, clusterCount);
755 tempClusterList.clear();
756
757 // Sort groups according to clusters, and set maxLimitingFactorForCluster
758 std::vector<distr_group_struct> tempGroupList(groupCount);
759 count = 0;
760 for(int i = 0; i < clusterCount; i++)
761 {
762 ergo_real maxLimitingFactorForCluster = 0;
763 int savedCount = count;
764 for(int j = clusterList[i].groupStartIndex; j < clusterList[i].groupStartIndex + clusterList[i].noOfGroups; j++)
765 {
766 ergo_real maxLimitingFactor = groupList[j].maxLimitingFactorGroup;
767 if(maxLimitingFactor > maxLimitingFactorForCluster)
768 maxLimitingFactorForCluster = maxLimitingFactor;
769 tempGroupList[count] = groupList[j];
770 count++;
771 } // END FOR j
772 clusterList[i].groupStartIndex = savedCount;
773 clusterList[i].maxLimitingFactorForCluster = maxLimitingFactorForCluster;
774 } // END FOR i
775 copy_vector<distr_group_struct>(groupList, tempGroupList, groupCount);
776 tempGroupList.clear();
777
778 // Sort distrs according to groups
779 std::vector<DistributionSpecStructLabeled> tempDistrList(distrCount);
780 //output_current_memory_usage("organize_distributions after allocating tempDistrList");
781 count = 0;
782 for(int i = 0; i < groupCount; i++)
783 {
784 int savedCount = count;
785 for(int j = groupList[i].startIndex; j < groupList[i].startIndex + groupList[i].distrCount; j++)
786 {
787 tempDistrList[count] = distrList[j];
788 count++;
789 } // END FOR j
790 groupList[i].startIndex = savedCount;
791 } // END FOR i
792 copy_vector<DistributionSpecStructLabeled>(distrList, tempDistrList, distrCount);
793 tempDistrList.clear();
794
795
796 result->minimalDistrList.resize(distrCount);
797 std::vector<minimal_distr_struct> & minimalDistrList = result->minimalDistrList;
798 for(int i = 0; i < distrCount; i++)
799 {
800 minimalDistrList[i].coeff = distrList[i].distr.coeff;
801 minimalDistrList[i].monomialIndex = integralInfo.monomial_info.monomial_index_list
802 [(int)distrList[i].distr.monomialInts[0]]
803 [(int)distrList[i].distr.monomialInts[1]]
804 [(int)distrList[i].distr.monomialInts[2]];
805 }
806
807
808 // get maxExtent
809 ergo_real maxExtent = 0;
810 for(int i = 0; i < distrCount; i++)
811 {
812 if(distrList[i].distr.extent > maxExtent)
813 maxExtent = distrList[i].distr.extent;
814 }
815 result->data.maxExtent = maxExtent;
816
817 // get maxDistanceOutsideBox
818 ergo_real maxDistanceOutsideBox = 0;
819 for(int i = 0; i < distrCount; i++)
820 {
821 // get minWallDist : minimum wall distance
822 ergo_real minWallDist = boxWidth;
823 int coordIndex;
824 for(coordIndex = 0; coordIndex< 3; coordIndex++)
825 {
826 // get wall distance for this coordinate
827 ergo_real dx = distrList[i].distr.centerCoords[coordIndex] - boxCenterCoords[coordIndex];
828 ergo_real wallDist = boxWidth / 2 - template_blas_fabs(dx);
829 if(wallDist < minWallDist)
830 minWallDist = wallDist;
831 } // END FOR coordIndex
832 if(minWallDist < -0.00001)
833 {
834 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: (minWallDist < -0.00001)");
835 return -1;
836 }
837 // Check that extent is nonzero, it should have been set before calling this routine, otherwise the maxDistanceOutsideBox info cannot be determined.
838 assert(distrList[i].distr.extent > 0);
839 ergo_real distanceOutsideBox = distrList[i].distr.extent - minWallDist;
840 if(distanceOutsideBox > maxDistanceOutsideBox)
841 maxDistanceOutsideBox = distanceOutsideBox;
842 }
843 result->data.maxDistanceOutsideBox = maxDistanceOutsideBox;
844
845 // Get maxNoOfMonomials
846 int maxNoOfMonomials = 0;
847 for(int i = 0; i < distrCount; i++) {
848 int degree = 0;
849 for(int j = 0; j < 3; j++)
850 degree += distrList[i].distr.monomialInts[j];
851 int noOfMonomials = integralInfo.monomial_info.no_of_monomials_list[degree];
852 if(noOfMonomials > maxNoOfMonomials)
853 maxNoOfMonomials = noOfMonomials;
854 }
855 result->data.maxNoOfMonomials = maxNoOfMonomials;
856
857 for(int i = 0; i < batchCount; i++)
858 for(int j = batchList[i].clusterStartIndex; j < batchList[i].clusterStartIndex + batchList[i].noOfClusters; j++)
859 {
860 int k_start = clusterList[j].groupStartIndex;
861 int k_end = k_start + clusterList[j].noOfGroups;
862 for(int k = k_start; k < k_end; k++)
863 {
864 int m_start = groupList[k].startIndex;
865 int m_end = m_start + groupList[k].distrCount;
866 for(int m = m_start; m < m_end; m++)
867 {
868 int foundIndex = -1;
869 for(int kk = 0; kk < batchList[i].noOfBasisFuncPairs; kk++)
870 {
871 if(basisFuncPairList[batchList[i].basisFuncPairListIndex+kk].index_1 == distrList[m].basisFuncIndex_1 &&
872 basisFuncPairList[batchList[i].basisFuncPairListIndex+kk].index_2 == distrList[m].basisFuncIndex_2)
873 {
874 foundIndex = kk;
875 break;
876 }
877 } // END FOR kk
878 if(foundIndex < 0)
879 {
880 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error setting basisFuncPairIndex");
881 return -1;
882 }
883 minimalDistrList[m].basisFuncPairIndex = foundIndex;
884 }
885 }
886 }
887
888
889 // within each group, sort minimalDistrList by monomialIndex and basisFuncPairIndex,
890 // join distrs that differ only in coefficient.
891 for(int i = 0; i < groupCount; i++)
892 {
893 minimal_distr_struct* p = & minimalDistrList[groupList[i].startIndex];
894 int count = groupList[i].distrCount;
895 for(int j = 0; j < count; j++)
896 for(int k = 0; k < count - 1 - j; k++)
897 {
898 int doSwitch = 0;
899 if(p[k].monomialIndex > p[k+1].monomialIndex)
900 doSwitch = 1;
901 if(p[k].monomialIndex == p[k+1].monomialIndex)
902 {
903 if(p[k].basisFuncPairIndex > p[k+1].basisFuncPairIndex)
904 doSwitch = 1;
905 }
906 if(doSwitch == 1)
907 {
908 minimal_distr_struct temp;
909 temp = p[k];
910 p[k] = p[k+1];
911 p[k+1] = temp;
912 }
913 } // END FOR j k
914 // OK, list sorted.
915 // We want to join together any entries that differ only in coefficient.
916 int j = 0;
917 int ii = 0;
918 while(ii < count)
919 {
920 ergo_real coeffSum = p[ii].coeff;
921 int k = ii + 1;
922 while(k < count)
923 {
924 if(p[k].monomialIndex != p[ii].monomialIndex || p[k].basisFuncPairIndex != p[ii].basisFuncPairIndex)
925 break;
926 coeffSum += p[k].coeff;
927 k++;
928 }
929 p[j] = p[ii];
930 p[j].coeff = coeffSum;
931 j++;
932 int nResult = k - ii;
933 ii += nResult;
934 }
935 groupList[i].distrCount = j;
936 } // END FOR i
937 // Now go through groups again to move the distrs together now that the groups are smaller.
938 count = 0;
939 for(int i = 0; i < groupCount; i++)
940 {
941 int oldStartIndex = groupList[i].startIndex;
942 groupList[i].startIndex = count;
943 int distrCount = groupList[i].distrCount;
944 for(int j = 0; j < distrCount; j++)
945 {
946 minimalDistrList[count] = minimalDistrList[oldStartIndex+j];
947 count++;
948 }
949 } // END FOR i
950 // check that no group contains repeating distrs
951 for(int i = 0; i < groupCount; i++)
952 {
953 minimal_distr_struct* p = & minimalDistrList[groupList[i].startIndex];
954 int distrCount = groupList[i].distrCount;
955 for(int j = 0; j < distrCount; j++)
956 for(int k = j+1; k < distrCount; k++)
957 {
958 if(p[k].monomialIndex == p[j].monomialIndex && p[k].basisFuncPairIndex == p[j].basisFuncPairIndex)
959 {
960 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: identical distrs found in same group.");
961 return -1;
962 }
963 }
964 }
965
966
967 int basisFuncForBatchsCount = 0;
968 // Now get list of basis func indeces occurring in each batch, store in basisFuncListForBatchs.
969 std::vector<int> basisFuncListForBatchs(2*distrCount);
970
971 for(int i = 0; i < batchCount; i++)
972 {
973
974 int count = 0;
975 for(int j = 0; j < batchList[i].noOfBasisFuncPairs; j++)
976 {
977 int i1 = basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_1;
978 int i2 = basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_2;
979 // Check if i1 and i2 are already present.
980 int i1_found = 0;
981 int i2_found = 0;
982 for(int k = 0; k < count; k++)
983 {
984 if(basisFuncListForBatchs[basisFuncForBatchsCount+k] == i1)
985 i1_found = 1;
986 if(basisFuncListForBatchs[basisFuncForBatchsCount+k] == i2)
987 i2_found = 1;
988 } // END FOR k
989 if(i1_found == 0)
990 {
991 basisFuncListForBatchs[basisFuncForBatchsCount+count] = i1;
992 count++;
993 }
994 if(i2_found == 0 && i1 != i2)
995 {
996 basisFuncListForBatchs[basisFuncForBatchsCount+count] = i2;
997 count++;
998 }
999 } // END FOR j
1000 // sort list for this batch
1001 do_sort_int_list(&basisFuncListForBatchs[basisFuncForBatchsCount], count);
1002 // now "rename" index_1 and index_2 using basisFuncListForBatchs.
1003 for(int j = 0; j < batchList[i].noOfBasisFuncPairs; j++)
1004 {
1005 int i1 = basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_1;
1006 int i2 = basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_2;
1007 // find positions of i1 and i2..
1008 int i1_index = -1;
1009 int i2_index = -1;
1010 for(int k = 0; k < count; k++)
1011 {
1012 if(basisFuncListForBatchs[basisFuncForBatchsCount+k] == i1)
1013 i1_index = k;
1014 if(basisFuncListForBatchs[basisFuncForBatchsCount+k] == i2)
1015 i2_index = k;
1016 } // END FOR k
1017 if(i1_index < 0 || i2_index < 0)
1018 {
1019 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: error 1!!!");
1020 return -1;
1021 }
1022 basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_1_mod = i1_index;
1023 basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_2_mod = i2_index;
1024 } // END FOR j
1025
1026 batchList[i].basisFuncForBatchsIndex = basisFuncForBatchsCount;
1027 batchList[i].basisFuncForBatchCount = count;
1028 basisFuncForBatchsCount += count;
1029 } // END FOR i
1030
1031 // Check result
1032 for(int i = 0; i < batchCount; i++)
1033 {
1034 for(int j = 0; j < batchList[i].noOfBasisFuncPairs; j++)
1035 {
1036 int i1 = basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_1;
1037 int i2 = basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_2;
1038 // Check if i1 and i2 are present.
1039 int i1_found = 0;
1040 int i2_found = 0;
1041 for(int k = 0; k < batchList[i].basisFuncForBatchCount; k++)
1042 {
1043 if(basisFuncListForBatchs[batchList[i].basisFuncForBatchsIndex+k] == i1)
1044 i1_found = 1;
1045 if(basisFuncListForBatchs[batchList[i].basisFuncForBatchsIndex+k] == i2)
1046 i2_found = 1;
1047 } // END FOR k
1048 if(i1_found == 0 || i2_found == 0)
1049 {
1050 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error: error !!");
1051 return -1;
1052 }
1053 } // END FOR j
1054 } // END FOR i
1055
1056
1057
1058 int basisFuncListCount = 0;
1059 // Now get list of basis func indices occurring, store in basisFuncList.
1060 // Use basisFuncListForBatchs to do this.
1061 std::vector<int> basisFuncList(basisFuncForBatchsCount);
1062 copy_vector<int>(basisFuncList, basisFuncListForBatchs, basisFuncForBatchsCount);
1063 std::sort(basisFuncList.begin(), basisFuncList.begin()+basisFuncForBatchsCount);
1064
1065 int prevIndex = -1;
1066 int i = 0;
1067 while(i < basisFuncForBatchsCount)
1068 {
1069 // now i points to a new basis func index.
1070 // check that sort order is OK.
1071 if(basisFuncList[i] < prevIndex)
1072 {
1073 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error 11! i = %i, basisFuncList[i] = %i, prevIndex = %i", i, basisFuncList[i], prevIndex);
1074 return -1;
1075 }
1076 basisFuncList[basisFuncListCount] = basisFuncList[i];
1077 basisFuncListCount++;
1078 prevIndex = basisFuncList[i];
1079 do i++; while(i < basisFuncForBatchsCount &&
1080 basisFuncList[i] == prevIndex);
1081 }
1082
1083 // Now go through batchs again to "rename" indices according to basisFuncList.
1084 for(int i = 0; i < batchCount; i++)
1085 {
1086 for(int j = 0; j < batchList[i].noOfBasisFuncPairs; j++)
1087 {
1088 int i1 = basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_1;
1089 int i2 = basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_2;
1090 // find positions of i1 and i2..
1091 int i1_index = -1;
1092 int i2_index = -1;
1093 for(int k = 0; k < basisFuncListCount; k++)
1094 {
1095 if(basisFuncList[k] == i1)
1096 i1_index = k;
1097 if(basisFuncList[k] == i2)
1098 i2_index = k;
1099 } // END FOR k
1100 if(i1_index < 0 || i2_index < 0)
1101 {
1102 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error 3!!!");
1103 return -1;
1104 }
1105 basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_inbox_1 = i1_index;
1106 basisFuncPairList[batchList[i].basisFuncPairListIndex+j].index_inbox_2 = i2_index;
1107 } // END FOR j
1108
1109 } // END FOR i
1110
1111 // take care of basisFuncListForBatchs_map
1112 result->basisFuncListForBatchs_map.resize(basisFuncForBatchsCount);
1113 for(int i = 0; i < basisFuncForBatchsCount; i++)
1114 {
1115 for(int k = 0; k < basisFuncListCount; k++)
1116 {
1117 if(basisFuncListForBatchs[i] == basisFuncList[k])
1118 result->basisFuncListForBatchs_map[i] = k;
1119 }
1120 }
1121
1122
1123 // take care of integral conversion matrices
1124 int noOfBasisFuncPairs_max = 0;
1125 for(int i = 0; i < batchCount; i++) {
1126 int noOfBasisFuncPairs = batchList[i].noOfBasisFuncPairs;
1127 if(noOfBasisFuncPairs > noOfBasisFuncPairs_max)
1128 noOfBasisFuncPairs_max = noOfBasisFuncPairs;
1129 }
1130 int noOfMonomials_max = 0;
1131 for(int i = 0; i < batchCount; i++)
1132 for(int j = batchList[i].clusterStartIndex; j < batchList[i].clusterStartIndex + batchList[i].noOfClusters; j++) {
1133 int noOfMonomials = integralInfo.monomial_info.no_of_monomials_list[clusterList[j].nmax];
1134 if(noOfMonomials > noOfMonomials_max)
1135 noOfMonomials_max = noOfMonomials;
1136 }
1137 // We do not know the final size of the spMatElementList yet. We give it some size to start with, and then increase the size when needed.
1138 std::vector<i_j_val_struct> spMatElementList(5*noOfBasisFuncPairs_max*noOfMonomials_max);
1139 int spMatElementListCount = 0;
1140 result->spMatCountList.resize(groupCount);
1141 result->spMatIdxList.resize(groupCount);
1142 for(int i = 0; i < batchCount; i++)
1143 for(int j = batchList[i].clusterStartIndex; j < batchList[i].clusterStartIndex + batchList[i].noOfClusters; j++) {
1144 // Now we are dealing with cluster j
1145 int noOfBasisFuncPairs = batchList[i].noOfBasisFuncPairs;
1146 int nmax = clusterList[j].nmax;
1147 int noOfMonomials = integralInfo.monomial_info.no_of_monomials_list[nmax];
1148 int group_k_start = clusterList[j].groupStartIndex;
1149 int group_k_end = group_k_start + clusterList[j].noOfGroups;
1150 // Get conversion matrices
1151 ergo_real alpha = groupList[group_k_start].exponent;
1152 i_j_val_struct convMat_sp[noOfMonomials*noOfMonomials];
1153 int convMat_nnz = integralInfo.get_hermite_conversion_matrix_right_sparse(nmax, 1.0/alpha, convMat_sp);
1154 for(int group_k = group_k_start; group_k < group_k_end; group_k++) {
1155 i_j_val_struct BB1_x_Ai1_x_convMat1_sp[noOfBasisFuncPairs*noOfMonomials];
1156 int BB1_x_Ai1_x_convMat1_nnz = 0;
1157 get_conversion_matrix_for_group(integralInfo,
1158 groupList[group_k],
1159 nmax,
1160 &minimalDistrList[0],
1161 noOfBasisFuncPairs,
1162 convMat_sp,
1163 convMat_nnz,
1164 BB1_x_Ai1_x_convMat1_sp, // result
1165 BB1_x_Ai1_x_convMat1_nnz); // result
1166 spmat_sort_elements(BB1_x_Ai1_x_convMat1_sp, BB1_x_Ai1_x_convMat1_nnz);
1167 // Check if the size of spMatElementList needs to be extended.
1168 while((int)(spMatElementList.size()) < spMatElementListCount + BB1_x_Ai1_x_convMat1_nnz)
1169 spMatElementList.resize(2*spMatElementList.size());
1170 memcpy(&spMatElementList[spMatElementListCount], BB1_x_Ai1_x_convMat1_sp, BB1_x_Ai1_x_convMat1_nnz*sizeof(i_j_val_struct));
1171 result->spMatCountList[group_k] = BB1_x_Ai1_x_convMat1_nnz;
1172 result->spMatIdxList[group_k] = spMatElementListCount;
1173 spMatElementListCount += BB1_x_Ai1_x_convMat1_nnz;
1174 if(spMatElementListCount > groupCount*noOfBasisFuncPairs_max*noOfMonomials_max)
1175 return -1;
1176 }
1177 }
1178
1179
1180 // Generate multipole limits for each group.
1181 for(int i = 0; i < batchCount; i++) {
1182 for(int j = batchList[i].clusterStartIndex; j < batchList[i].clusterStartIndex + batchList[i].noOfClusters; j++) {
1183 // Now we are dealing with cluster j
1184 int group_start = clusterList[j].groupStartIndex;
1185 int group_end = group_start + clusterList[j].noOfGroups;
1186 for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++)
1187 clusterList[j].multipoleEuclNormListForK[l] = 0;
1188 for(int groupIndex = group_start; groupIndex < group_end; groupIndex++) {
1189 distr_group_struct* currGroup = &groupList[groupIndex];
1190 ergo_real maxMomentVectorNormForDistrsList[MAX_MULTIPOLE_DEGREE_BASIC+1];
1191 for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++)
1192 maxMomentVectorNormForDistrsList[l] = 0;
1193 int distr_start = currGroup->startIndex;
1194 int distr_end = distr_start + currGroup->distrCount;
1195 for(int distrIndex = distr_start; distrIndex < distr_end; distrIndex++) {
1196 int monomialIndex = minimalDistrList[distrIndex].monomialIndex;
1197 ergo_real coeff = minimalDistrList[distrIndex].coeff;
1198 // get monomialInts from monomialIndex
1199 DistributionSpecStruct distr;
1200 distr.monomialInts[0] = integralInfo.monomial_info.monomial_list[monomialIndex].ix;
1201 distr.monomialInts[1] = integralInfo.monomial_info.monomial_list[monomialIndex].iy;
1202 distr.monomialInts[2] = integralInfo.monomial_info.monomial_list[monomialIndex].iz;
1203 distr.coeff = coeff;
1204 distr.exponent = currGroup->exponent;
1205 distr.centerCoords[0] = currGroup->centerCoords[0];
1206 distr.centerCoords[1] = currGroup->centerCoords[1];
1207 distr.centerCoords[2] = currGroup->centerCoords[2];
1208 multipole_struct_small multipole;
1209 if(compute_multipole_moments(integralInfo, &distr, &multipole) != 0) {
1210 do_output(LOG_CAT_ERROR, LOG_AREA_INTEGRALS, "error in compute_multipole_moments");
1211 return -1;
1212 }
1213 // modfy maxMomentVectorNormForDistrsList if needed.
1214 for(int l = 0; l <= multipole.degree; l++) {
1215 int startIndex = l*l;
1216 int endIndex = (l+1)*(l+1);
1217 ergo_real sum = 0;
1218 for(int A = startIndex; A < endIndex; A++)
1219 sum += multipole.momentList[A]*multipole.momentList[A];
1220 ergo_real subNorm = template_blas_sqrt(sum);
1221 if(subNorm > maxMomentVectorNormForDistrsList[l])
1222 maxMomentVectorNormForDistrsList[l] = subNorm;
1223 }
1224 } // END FOR distrIndex
1225 for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++)
1226 currGroup->multipoleEuclNormListForK[l] = maxMomentVectorNormForDistrsList[l];
1227 for(int l = 0; l <= MAX_MULTIPOLE_DEGREE_BASIC; l++) {
1228 if(currGroup->multipoleEuclNormListForK[l] > clusterList[j].multipoleEuclNormListForK[l])
1229 clusterList[j].multipoleEuclNormListForK[l] = currGroup->multipoleEuclNormListForK[l];
1230 }
1231 } // end for groupIndex -- loop over groups
1232 } // end for j -- loop over clusters
1233 } // end for i -- loop over batches
1234
1235
1236 result->spMatElementList.resize(spMatElementListCount);
1237 copy_vector<i_j_val_struct>(result->spMatElementList, spMatElementList, spMatElementListCount);
1238
1239 result->clusterList.resize(clusterCount);
1240 copy_vector<cluster_struct>(result->clusterList, clusterList, clusterCount);
1241
1242 result->batchList.resize(batchCount);
1243 copy_vector<batch_struct>(result->batchList, batchList, batchCount);
1244
1245 result->basisFuncPairList.resize(basisFuncPairCount);
1246 copy_vector<basis_func_pair_struct>(result->basisFuncPairList, basisFuncPairList, basisFuncPairCount);
1247
1248 result->basisFuncListForBatchs.resize(basisFuncForBatchsCount);
1249 copy_vector<int>(result->basisFuncListForBatchs, basisFuncListForBatchs, basisFuncForBatchsCount);
1250
1251 result->basisFuncList.resize(basisFuncListCount);
1252 copy_vector<int>(result->basisFuncList, basisFuncList, basisFuncListCount);
1253
1254 // Do memcpy only if distrCount > 0 to avoid problem if -D_GLIBCXX_ASSERTIONS is used
1255 if(distrCount > 0)
1256 memcpy(&distrList_in[0], &distrList[0], distrCount*sizeof(DistributionSpecStructLabeled));
1257
1258 return 0;
1259 } // END organize_distributions
1260
1261