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