Attention une mise à jour du serveur va être effectuée le lundi 17 mai entre 13h et 13h30. Cette mise à jour va générer une interruption du service de quelques minutes.

Commit 0bb85290 authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Update in API in order to remove memory leaks (changes added to FChebInterpolator too.)

parent 0b6ed9f6
......@@ -465,7 +465,7 @@ void scalfmm_execute_fmm(scalfmm_handle Handle);
* The last param cellDestroyer is meaningless in case the user uses
* one of the provided kernel. (i.e. Chebyshev, Lagrange)
*/
void Scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell cellDestroyer);
void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell cellDestroyer);
#endif
......@@ -24,9 +24,10 @@
#include "FScalFMMEngine.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/P2P/FP2PLeafInterface.hpp"
//#include "Kernels/P2P/FP2PLeafInterface.hpp"
#include "Arranger/FOctreeArranger.hpp"
#include "Arranger/FArrangerPeriodic.hpp"
#include "Arranger/FBasicParticleContainerIndexedMover.hpp"
#include "Core/FFmmAlgorithmThread.hpp"
#include "Core/FFmmAlgorithm.hpp"
......@@ -45,12 +46,14 @@ class FInterEngine : public FScalFMMEngine{
private:
//Typedef on the octree class, in order to clarify following code
typedef FOctree<InterCell,ContainerClass,LeafClass> OctreeClass;
typedef FP2PLeafInterface<OctreeClass> LeafInterface;
//typedef FP2PLeafInterface<OctreeClass> LeafInterface;
//Typedef on Octree Arranger, in order to clarify following code
typedef FOctreeArranger<OctreeClass,ContainerClass,LeafInterface> ArrangerClass;
typedef FArrangerPeriodic<OctreeClass,ContainerClass,LeafInterface> ArrangerClassPeriodic;
typedef FBasicParticleContainerIndexedMover<OctreeClass, ContainerClass> MoverClass;
typedef FOctreeArranger<OctreeClass, ContainerClass, MoverClass> ArrangerClass;
typedef FArrangerPeriodic<OctreeClass, ContainerClass, MoverClass> ArrangerClassPeriodic;
//Pointer to the kernel to be executed
InterKernel * kernel;
MatrixKernelClass * matrix;
......@@ -78,8 +81,12 @@ public:
//TODO free kernel too
~FInterEngine(){
free(octree);
free(kernel);
delete matrix;
delete octree;
delete kernel;
if(arranger){
delete arranger;
}
}
//Inserting array of position
......@@ -528,7 +535,6 @@ public:
void execute_fmm(){
//typedef FOctree<InterCell,ContainerClass,LeafClass> OctreeClass;
switch(Algorithm){
case 0:
{
......@@ -557,7 +563,9 @@ public:
}
}
void intern_dealloc_handle(Callback_free_cell unUsed){
//this->~FInterEngine();
}
};
......
......@@ -51,11 +51,16 @@ protected:
scalfmm_kernel_type kernelType;
scalfmm_algorithm Algorithm;
FVector<bool> progress;
FVector<bool>* progress;
int nbPart;
public:
FScalFMMEngine() : Algorithm(multi_thread), progress(), nbPart(0){
FScalFMMEngine() : Algorithm(multi_thread), progress(nullptr), nbPart(0){
progress = new FVector<bool>();
}
virtual ~FScalFMMEngine() {
delete progress;
}
//First function displayed there are common function for every
......@@ -208,6 +213,10 @@ public:
FAssertLF(0,"No kernel set, cannot execute anything, exiting ...\n");
}
virtual void intern_dealloc_handle(Callback_free_cell userDeallocator){
FAssertLF(0,"No kernel set, cannot execute anything, exiting ...\n");
}
};
......@@ -378,6 +387,10 @@ extern "C" void scalfmm_execute_fmm(scalfmm_handle Handle){
((ScalFmmCoreHandle * ) Handle)->engine->execute_fmm();
}
extern "C" void scalfmm_user_kernel_config(scalfmm_handle Handle, Scalfmm_Kernel_Descriptor userKernel, void * userDatas){
((ScalFmmCoreHandle * ) Handle)->engine->user_kernel_config(userKernel,userDatas);
}
extern "C" void scalfmm_init_cell(scalfmm_handle Handle, Callback_init_cell user_cell_initializer){
((ScalFmmCoreHandle * ) Handle)->engine->init_cell(user_cell_initializer);
}
......@@ -386,5 +399,8 @@ extern "C" void scalfmm_free_cell(scalfmm_handle Handle, Callback_free_cell user
((ScalFmmCoreHandle * ) Handle)->engine->free_cell(user_cell_deallocator);
}
// extern "C" void scalfmm_intern_dealloc_handle(scalfmm_handle Handle,Callback_free_cell userDeallocator){
// ((ScalFmmCoreHandle * ) Handle)->engine->intern_dealloc_handle(userDeallocator);
// }
#endif
#ifndef CALL_HPP
#define CALL_HPP
/** It should be compiled with C export */
extern "C" {
#include "CScalfmmApi.h"
......@@ -15,7 +11,7 @@ extern "C" scalfmm_handle scalfmm_init(int TreeHeight,double BoxWidth,double* Bo
switch(KernelType){
case 0:
//handle->engine = new FUserKernelEngine(...);
handle->engine = new FUserKernelEngine(TreeHeight, BoxWidth, BoxCenter, KernelType);
std::cout<< "User Kernel type unsupported yet" << std::endl;
break;
......@@ -48,5 +44,8 @@ extern "C" scalfmm_handle scalfmm_init(int TreeHeight,double BoxWidth,double* Bo
return handle;
}
#endif
extern "C" void scalfmm_dealloc_handle(scalfmm_handle handle, Callback_free_cell userDeallocator){
((ScalFmmCoreHandle *) handle)->engine->intern_dealloc_handle(userDeallocator);
delete ((ScalFmmCoreHandle *) handle)->engine ;
delete (ScalFmmCoreHandle *) handle;
}
......@@ -153,10 +153,11 @@ private:
typedef FSimpleLeaf<ContainerClass> LeafClass;
typedef FOctree<CoreCell,ContainerClass,LeafClass> OctreeClass;
typedef CoreKernel<CoreCell,ContainerClass> CoreKernelClass;
typedef FP2PLeafInterface<OctreeClass> LeafInterface;
//For arranger classes
typedef FOctreeArranger<OctreeClass,ContainerClass,LeafInterface> ArrangerClass;
typedef FArrangerPeriodic<OctreeClass,ContainerClass,LeafInterface> ArrangerClassPeriodic;
typedef FBasicParticleContainerIndexedMover<OctreeClass, ContainerClass> MoverClass;
typedef FOctreeArranger<OctreeClass, ContainerClass, MoverClass> ArrangerClass;
typedef FArrangerPeriodic<OctreeClass, ContainerClass, MoverClass> ArrangerClassPeriodic;
//Attributes
......@@ -171,12 +172,23 @@ public:
kernelType = KernelType;
//Kernel is not set now because the user must provide a
//Scalfmm_Kernel_descriptor
//kernel = new CoreKernelClass();
}
~FUserKernelEngine(){
delete octree;
if(arranger){
delete arranger;
}
if(kernel){
delete kernel;
}
}
void user_kernel_config( Scalfmm_Kernel_Descriptor userKernel, void * userDatas){
kernel = new CoreKernelClass(userKernel,userDatas);
if(!kernel){
kernel = new CoreKernelClass(userKernel,userDatas);
}
}
......@@ -226,11 +238,11 @@ public:
}
void execute_fmm(){
//FAssertLF(kernel,"No kernel set, please use scalfmm_user_kernel_config before calling the execute routine ... Exiting \n");
FAssertLF(kernel,"No kernel set, please use scalfmm_user_kernel_config before calling the execute routine ... Exiting \n");
switch(Algorithm){
case 0:
{
typedef FFmmAlgorithm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafInterface> AlgoClassSeq;
typedef FFmmAlgorithm<OctreeClass,CoreCell,ContainerClass,CoreKernelClass,LeafClass> AlgoClassSeq;
AlgoClassSeq algoSeq(octree,kernel);
algoSeq.execute();
break;
......@@ -255,7 +267,10 @@ public:
}
}
void intern_dealloc_handle(Callback_free_cell userDeallocator){
free_cell(userDeallocator);
this->~FUserKernelEngine();
}
};
......
......@@ -23,7 +23,7 @@ int main(int argc, char ** av){
scalfmm_kernel_type myChoice = chebyshev;
//Octree configuration
int TreeHeight = 5;
int TreeHeight = 7;
double boxWidth = 2.0;
double boxCenter[3] = {0.0,0.0,0.0};
......@@ -56,7 +56,7 @@ int main(int argc, char ** av){
//Computation Part
int nb_iteration = 100;
int nb_iteration = 10;//atoi(av[1]);
int curr_iteration = 0;
//Array to store the output
......@@ -112,13 +112,16 @@ int main(int argc, char ** av){
curr_iteration++;
}
//End of Computation, useer get the position after a hundreds of iterations
//Free memory
free(positionsXYZ);
free(array_of_charge);
free(array_of_forces);
free(array_of_pot);
free(array_of_displacement);
//End of Computation, useer get the position after some iterations
scalfmm_dealloc_handle(handle,NULL);
return 0;
}
......@@ -236,10 +236,10 @@ int main(int argc, char ** argv){
scalfmm_handle handle = scalfmm_init(treeHeight, boxWidth, boxCenter,user_defined_kernel);
// Insert particles
printf("Inserting particles...\n");
Scalfmm_insert_array_of_particles(handle, nbParticles, particleIndexes, particleXYZ);
scalfmm_tree_insert_particles_xyz(handle, nbParticles, particleXYZ);
// Init cell
printf("Initizalizing the cells:\n");
Scalfmm_init_cell(handle, my_Callback_init_cell);
scalfmm_init_cell(handle, my_Callback_init_cell);
// Init our callback struct
struct User_Scalfmm_Kernel_Descriptor kernel;
......@@ -255,9 +255,11 @@ int main(int argc, char ** argv){
struct MyData my_data;
memset(&my_data, 0, sizeof(struct MyData));
my_data.insertedPositions = particleXYZ;
//Set my datas before calling fmm (this will set as well the kernel)
scalfmm_user_kernel_config(handle,kernel,&my_data);
// Execute the FMM
Scalfmm_execute_kernel(handle, kernel, &my_data);
scalfmm_execute_fmm(handle/*, kernel, &my_data*/);
// Print the results store in our callback
printf("There was %d P2M\n", my_data.countP2M);
......@@ -268,8 +270,10 @@ int main(int argc, char ** argv){
printf("There was %d P2PInner\n", my_data.countP2PInner);
printf("There was %d P2P\n", my_data.countP2P);
//Dealloc the cells
scalfmm_free_cell(handle,my_Callback_free_cell);
// Dealloc the handle
Scalfmm_dealloc_handle(handle,my_Callback_free_cell);
scalfmm_dealloc_handle(handle,my_Callback_free_cell);
return 0;
}
......@@ -4,13 +4,13 @@
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FCHEBINTERPOLATOR_HPP
......@@ -36,9 +36,9 @@
*
* The class @p FChebInterpolator defines the anterpolation (M2M) and
* interpolation (L2L) concerning operations.
*
* PB: MatrixKernelClass is passed as template in order to inform interpolator
* of the size of the vectorial interpolators. Default is the scalar
*
* PB: MatrixKernelClass is passed as template in order to inform interpolator
* of the size of the vectorial interpolators. Default is the scalar
* matrix kernel class of type ONE_OVER_R (NRHS=NLHS=1).
*/
template <int ORDER, class MatrixKernelClass = struct FInterpMatrixKernelR>
......@@ -58,7 +58,7 @@ protected: // PB for OptiDis
FReal T[ORDER * (ORDER-1)];
unsigned int node_ids[nnodes][3];
// 8 Non-leaf (i.e. M2M/L2L) interpolators
// 8 Non-leaf (i.e. M2M/L2L) interpolators
// x1 per level if box is extended
// only 1 is required for all levels if extension is 0
FReal*** ChildParentInterpolator;
......@@ -279,7 +279,7 @@ protected: // PB for OptiDis
FPoint ChildRoots[nnodes];
// Ratio of extended cell widths (definition: child ext / parent ext)
const FReal ExtendedCellRatio =
const FReal ExtendedCellRatio =
FReal(FReal(ParentWidth)/FReal(2.) + CellWidthExtension) / FReal(ParentWidth + CellWidthExtension);
// Child cell center and width
......@@ -301,7 +301,7 @@ protected: // PB for OptiDis
}
}
/**
* Initialize the child - parent - interpolator, it is basically the matrix
* S which is precomputed and reused for all M2M and L2L operations, ie for
......@@ -312,11 +312,11 @@ protected: // PB for OptiDis
*/
void initTensorM2MandL2L(const int TreeLevel, const FReal ParentWidth)
{
FReal ChildCoords[3][ORDER];
FReal ChildCoords[3][ORDER];
FPoint ChildCenter;
// Ratio of extended cell widths (definition: child ext / parent ext)
const FReal ExtendedCellRatio =
const FReal ExtendedCellRatio =
FReal(FReal(ParentWidth)/FReal(2.) + CellWidthExtension) / FReal(ParentWidth + CellWidthExtension);
// Child cell width
......@@ -325,7 +325,7 @@ protected: // PB for OptiDis
// loop: child cells
for (unsigned int child=0; child<8; ++child) {
// set child info
// set child info
FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter, ExtendedCellRatio);
FChebTensor<ORDER>::setPolynomialsRoots(ChildCenter, ChildWidth, ChildCoords);
// allocate memory
......@@ -359,13 +359,13 @@ public:
*
* PB: Input parameters ONLY affect the computation of the M2M/L2L ops.
* These parameters are ONLY required in the context of extended bbox.
* If no M2M/L2L is required then the interpolator can be built with
* If no M2M/L2L is required then the interpolator can be built with
* the default ctor.
*/
explicit FChebInterpolator(const int inTreeHeight=3,
const FReal inRootCellWidth=FReal(1.),
const FReal inCellWidthExtension=FReal(0.))
: TreeHeight(inTreeHeight),
const FReal inRootCellWidth=FReal(1.),
const FReal inCellWidthExtension=FReal(0.))
: TreeHeight(inTreeHeight),
RootCellWidth(inRootCellWidth),
CellWidthExtension(inCellWidthExtension)
{
......@@ -390,27 +390,27 @@ public:
for (unsigned int l=0; l<static_cast<unsigned int>(TreeHeight); ++l){
ChildParentInterpolator[l] = new FReal*[8];
for (unsigned int c=0; c<8; ++c)
ChildParentInterpolator[l][c]=nullptr;
ChildParentInterpolator[l][c]=nullptr;
}
// Set number of non-leaf ios that actually need to be computed
unsigned int reducedTreeHeight; // = 2 + nb of computed nl ios
if(CellWidthExtension==0.) // if no cell extension, then ...
reducedTreeHeight = 3; // cmp only 1 non-leaf io
else
else
reducedTreeHeight = TreeHeight; // cmp 1 non-leaf io per level
// Init non-leaf interpolators
FReal CellWidth = RootCellWidth / FReal(2.); // at level 1
CellWidth /= FReal(2.); // at level 2
for (unsigned int l=2; l<reducedTreeHeight; ++l) {
//this -> initM2MandL2L(l,CellWidth); // non tensor-product interpolation
this -> initTensorM2MandL2L(l,CellWidth); // tensor-product interpolation
// update cell width
CellWidth /= FReal(2.); // at level l+1
CellWidth /= FReal(2.); // at level l+1
}
}
......@@ -420,10 +420,15 @@ public:
*/
~FChebInterpolator()
{
for (unsigned int l=0; l<static_cast<unsigned int>(TreeHeight); ++l)
for (unsigned int child=0; child<8; ++child)
if(ChildParentInterpolator[l][child] != nullptr)
delete [] ChildParentInterpolator[l][child];
for (unsigned int l=0; l<static_cast<unsigned int>(TreeHeight); ++l){
for (unsigned int child=0; child<8; ++child){
if(ChildParentInterpolator[l][child] != nullptr){
delete [] ChildParentInterpolator[l][child];
}
}
delete[] ChildParentInterpolator[l];
}
delete[] ChildParentInterpolator;
}
......@@ -702,7 +707,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
} // flops: (ORDER-1) * 2
} // flops: (ORDER-1) * (6 + (ORDER-1) * 2)
} // flops: (ORDER-1) * (6 + (ORDER-1) * (6 + (ORDER-1) * 2))
} // flops: ... * NRHS
} // flops: ... * NRHS
} // flops: N * (18 + (ORDER-2) * 6 + (ORDER-1) * (6 + (ORDER-1) * (6 + (ORDER-1) * 2)))
////////////////////////////////////////////////////////////////////
......@@ -768,9 +773,9 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
//template <int ORDER>
//template <class ContainerClass>
//inline void FChebInterpolator<ORDER>::applyP2M(const FPoint& center,
// const FReal width,
// FReal *const multipoleExpansion,
// const ContainerClass *const sourceParticles) const
// const FReal width,
// FReal *const multipoleExpansion,
// const ContainerClass *const sourceParticles) const
//{
// // set all multipole expansions to zero
// FBlas::setzero(nnodes, multipoleExpansion);
......@@ -803,12 +808,12 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyP2M(const FPoint& c
// T_of_x[o][1] = ypy * T_of_x[o-1][1] - T_of_x[o-2][1]; // 2 flops
// T_of_x[o][2] = zpz * T_of_x[o-1][2] - T_of_x[o-2][2]; // 2 flops
// } // flops: (ORDER-1) * 6
//
//
// // anterpolate
// const FReal sourceValue = iter.data().getPhysicalValue();
// for (unsigned int n=0; n<nnodes; ++n) {
// const unsigned int j[3] = {node_ids[n][0], node_ids[n][1], node_ids[n][2]};
// S[0] = FReal(0.5) + T_of_x[1][0] * T_of_roots[1][j[0]]; // 2 flops
// S[0] = FReal(0.5) + T_of_x[1][0] * T_of_roots[1][j[0]]; // 2 flops
// S[1] = FReal(0.5) + T_of_x[1][1] * T_of_roots[1][j[1]]; // 2 flops
// S[2] = FReal(0.5) + T_of_x[1][2] * T_of_roots[1][j[2]]; // 2 flops
// for (unsigned int o=2; o<ORDER; ++o) {
......@@ -842,12 +847,12 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
FReal W2[nLhs][3][ ORDER-1];
FReal W4[nLhs][3][(ORDER-1)*(ORDER-1)];
FReal W8[nLhs][ (ORDER-1)*(ORDER-1)*(ORDER-1)];
{
{
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// sum over interpolation points
f1[idxLhs] = FReal(0.);
for(unsigned int i=0; i<ORDER-1; ++i) W2[idxLhs][0][i] = W2[idxLhs][1][i] = W2[idxLhs][2][i] = FReal(0.);
for(unsigned int i=0; i<ORDER-1; ++i) W2[idxLhs][0][i] = W2[idxLhs][1][i] = W2[idxLhs][2][i] = FReal(0.);
for(unsigned int i=0; i<(ORDER-1)*(ORDER-1); ++i) W4[idxLhs][0][i] = W4[idxLhs][1][i] = W4[idxLhs][2][i] = FReal(0.);
for(unsigned int i=0; i<(ORDER-1)*(ORDER-1)*(ORDER-1); ++i) W8[idxLhs][i] = FReal(0.);
......@@ -916,13 +921,13 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
}
}
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
// distribution over potential components:
// We sum the multidim contribution of PhysValue
// This was originally done at M2L step but moved here
// This was originally done at M2L step but moved here
// because their storage is required by the force computation.
// In fact : f_{ik}(x)=w_j(x) \nabla_{x_i} K_{ij}(x,y)w_j(y))
const unsigned int idxPot = idxLhs / nPV;
const unsigned int idxPot = idxLhs / nPV;
FReal*const potentials = inParticles->getPotentials(idxPot);
......@@ -970,7 +975,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
// FReal F2[(ORDER-1) * ORDER*ORDER];
// FBlas::gemtm(ORDER, ORDER-1, ORDER*ORDER, FReal(1.), const_cast<FReal*>(T), ORDER,
// const_cast<FReal*>(localExpansion), ORDER, F2, ORDER-1);
// const_cast<FReal*>(localExpansion), ORDER, F2, ORDER-1);
// FReal F[ORDER-1]; FBlas::setzero(ORDER-1, F);
// for (unsigned int i=0; i<ORDER-1; ++i)
// for (unsigned int j=0; j<ORDER*ORDER; ++j) F[i] += F2[j*(ORDER-1) + i];
......@@ -984,9 +989,9 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
//template <int ORDER>
//template <class ContainerClass>
//inline void FChebInterpolator<ORDER>::applyL2P(const FPoint& center,
// const FReal width,
// const FReal *const localExpansion,
// ContainerClass *const localParticles) const
// const FReal width,
// const FReal *const localExpansion,
// ContainerClass *const localParticles) const
//{
// // allocate stuff
// const map_glob_loc map(center, width);
......@@ -998,7 +1003,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
// c1 = FReal(8.) / nnodes ;
// typename ContainerClass::BasicIterator iter(*localParticles);
// while(iter.hasNotFinished()){
//
//
// // map global position to [-1,1]
// map(iter.data().getPosition(), localPosition); // 15 flops
//
......@@ -1026,7 +1031,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2P(const FPoint& c
// S[0] += T_of_x[o][0] * T_of_roots[o][j[0]]; // 2 flops
// S[1] += T_of_x[o][1] * T_of_roots[o][j[1]]; // 2 flops
// S[2] += T_of_x[o][2] * T_of_roots[o][j[2]]; // 2 flops
// } // (ORDER-2) * 6 flops
// } // (ORDER-2) * 6 flops
// // gather contributions
// S[0] += FReal(0.5); // 1 flops
// S[1] += FReal(0.5); // 1 flops
......@@ -1115,7 +1120,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
}
// apply P and increment forces
FReal forces[nLhs][3];
FReal forces[nLhs][3];
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs)
for (unsigned int i=0; i<3; ++i)
forces[idxLhs][i] = FReal(0.);
......@@ -1172,8 +1177,8 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PGradient(const F
}
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
const unsigned int idxPot = idxLhs / nPV;
const unsigned int idxPV = idxLhs % nPV;
const unsigned int idxPot = idxLhs / nPV;
const unsigned int idxPV = idxLhs % nPV;
// scale forces
forces[idxLhs][0] *= jacobian[0] / nnodes;
......@@ -1213,7 +1218,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi
//{ // sum over interpolation points
// f1 = FReal(0.);
// for(unsigned int i=0; i<ORDER-1; ++i) W2[0][i] = W2[1][i] = W2[2][i] = FReal(0.);
// for(unsigned int i=0; i<ORDER-1; ++i) W2[0][i] = W2[1][i] = W2[2][i] = FReal(0.);
// for(unsigned int i=0; i<(ORDER-1)*(ORDER-1); ++i) W4[0][i] = W4[1][i] = W4[2][i] = FReal(0.);
// for(unsigned int i=0; i<(ORDER-1)*(ORDER-1)*(ORDER-1); ++i) W8[i] = FReal(0.);
//
......@@ -1263,7 +1268,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
f1[idxLhs] = FReal(0.);
for (unsigned int idx=0; idx<nnodes; ++idx) f1[idxLhs] += localExpansion[2*idxLhs*nnodes + idx]; // 1 flop
//////////////////////////////////////////////////////////////////
// IMPORTANT: NOT CHANGE ORDER OF COMPUTATIONS!!! ////////////////
//////////////////////////////////////////////////////////////////
......@@ -1361,7 +1366,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi
// apply P and increment forces
FReal potential[nLhs];
FReal forces[nLhs][3];
FReal forces[nLhs][3];
for(int idxLhs = 0 ; idxLhs < nLhs ; ++idxLhs){
potential[idxLhs]= FReal(0.);
for (unsigned int i=0; i<3; ++i)
......@@ -1409,7 +1414,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi
const int idxPot = idxLhs / nPV;
const int idxPV = idxLhs % nPV;
// get potentials, physValues and forces components
// get potentials, physValues and forces components
const FReal*const physicalValues = inParticles->getPhysicalValues(idxPV);
FReal*const forcesX = inParticles->getForcesX(idxPot);
FReal*const forcesY = inParticles->getForcesY(idxPot);
......@@ -1445,7 +1450,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi
// const map_glob_loc map(center, width);
// FPoint Jacobian;
// map.computeJacobian(Jacobian); // 6 flops
// const FReal jacobian[3] = {Jacobian.getX(), Jacobian.getY(), Jacobian.getZ()};
// const FReal jacobian[3] = {Jacobian.getX(), Jacobian.getY(), Jacobian.getZ()};
// FPoint localPosition;
// FReal T_of_x[ORDER][3];
// FReal U_of_x[ORDER][3];
......@@ -1456,10 +1461,10 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi
// //
// typename ContainerClass::BasicIterator iter(*localParticles);
// while(iter.hasNotFinished()){
//
//
// // map global position to [-1,1]
// map(iter.data().getPosition(), localPosition); // 15 flops
//
//
// // evaluate chebyshev polynomials of source particle
// // T_0(x_i) and T_1(x_i)
// xpx = FReal(2.) * localPosition.getX(); // 1 flop
......@@ -1473,7 +1478,7 @@ inline void FChebInterpolator<ORDER,MatrixKernelClass>::applyL2PTotal(const FPoi
// U_of_x[0][1] = FReal(1.); U_of_x[1][1] = ypy;
// U_of_x[0][2] = FReal(1.); U_of_x[1][2] = zpz;