Commit daf5fa5b authored by PIACIBELLO Cyrille's avatar PIACIBELLO Cyrille

Most efficient version is Kernel_Boucle3

parent e883708e
......@@ -469,7 +469,6 @@ public:
void P2M(CellClass* const pole,
const ContainerClass* const particles)
{
// printf("P2M %%%%%%%%%%%%%%%%%%%%%%%%%%\n");
//Variables computed for each power of Multipole
int a,b,c ;
FReal facto, coeff;
......@@ -815,7 +814,7 @@ public:
{
// printf("L2P %%%%%%%%%%%%%%%%%%%%%%%%%%\n");
FPoint locCenter = getLeafCenter(local->getCoordinate());
//Iterator over particles
int nbPart = particles->getNbParticles();
//
......@@ -842,6 +841,7 @@ public:
FReal dx = posX[i] - locCenter.getX();
FReal dy = posY[i] - locCenter.getY();
FReal dz = posZ[i] - locCenter.getZ();
//
// Precompute an arrays of Array[i] = dx^(i-1)
arrayDX[0] = 0.0 ;
......@@ -859,7 +859,10 @@ public:
}
FReal partPhyValue = phyValues[i];
//
FReal locPot = 0.0, locForceX = 0.0, locForceY = 0.0, locForceZ = 0.0 ;
FReal locPot = 0.0;
FReal locForceX = 0.0;
FReal locForceY = 0.0;
FReal locForceZ = 0.0;
int a=0,b=0,c=0;
for(int j=0 ; j<SizeVector ; ++j){
FReal locForce = iterLocal[j];
......@@ -869,7 +872,7 @@ public:
locForceX += FReal(a)*locForce*arrayDX[a]*arrayDY[b+1]*arrayDZ[c+1];
locForceY += FReal(b)*locForce*arrayDX[a+1]*arrayDY[b]*arrayDZ[c+1];
locForceZ += FReal(c)*locForce*arrayDX[a+1]*arrayDY[b+1]*arrayDZ[c];
//
incPowers(&a,&b,&c);
}
targetsPotentials[i] = partPhyValue*locPot ;
......
This diff is collapsed.
This diff is collapsed.
......@@ -113,8 +113,6 @@ private:
FPoint cCenter = FPoint(X,Y,Z);
//For debug purpose
//printf("%f,%f,%f\n",cCenter.getX(),cCenter.getY(),cCenter.getZ());
return cCenter;
}
......@@ -203,7 +201,6 @@ private:
void initDerivative(const FReal & dx ,const FReal & dy ,const FReal & dz , FReal * tab)
{
FReal R2 = dx*dx+dy*dy+dz*dz;
// printf("dx : %f dy : %f dz : %f\n",dx,dy,dz);
tab[0]=FReal(1)/FMath::Sqrt(R2);
FReal R3 = tab[0]/(R2);
tab[1]= -dx*R3; //Derivative in (1,0,0) il doit y avoir un -
......@@ -216,11 +213,6 @@ private:
tab[7] = FReal(3)*dy*dy*R5-R3; //Derivative in (0,2,0)
tab[8] = FReal(3)*dy*dz*R5; //Derivative in (0,1,1)
tab[9] = FReal(3)*dz*dz*R5-R3; //Derivative in (0,0,2)
// For debugging purpose : print the values computed
// for(int c=0 ; c<=9 ; ++c){
// printf("just computed %f, a=%d, b=%d, c=%d target: %d %f\n",dx,0,0,0,c,tab[c]);
// }
}
/** @brief Compute and store the derivative for a given tuple.
......@@ -487,17 +479,14 @@ public:
void P2M(CellClass* const pole,
const ContainerClass* const particles)
{
// printf("P2M %%%%%%%%%%%%%%%%%%%%%%%%%%\n");
//Copying cell center position once and for all
const FPoint& cellCenter = getLeafCenter(pole->getCoordinate());
// printf("P2M :: pole : X: %f, Y: %f, Z:%f \n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ());
FReal * FRestrict multipole = pole->getMultipole();
FMemUtils::memset(multipole,0,SizeVector*sizeof(FReal(0.0)));
//
FReal multipole2[SizeVector] ;
//
// Iterator over Particles
//
int nbPart = particles->getNbParticles(), i;
const FReal* const * positions = particles->getPositions();
const FReal* posX = positions[0];
......@@ -511,14 +500,13 @@ public:
FReal xc = cellCenter.getX(), yc = cellCenter.getY(), zc = cellCenter.getZ() ;
FReal dx[3] ;
for(int idPart=0 ; idPart<nbPart ; ++idPart){
// printf("P2M :: part : X: %f, Y: %f, Z:%f Q %f\n", posX[idPart] , posY[idPart] , posZ[idPart] ,phyValue[idPart]);
//
dx[0] = xc - posX[idPart] ;
dx[1] = yc - posY[idPart] ;
dx[2] = zc - posZ[idPart] ;
// printf(" : DX: %f, DY: %f, DZ:%f \n",dx[0] ,dx[01] ,dx[2] );
multipole2[0] = phyValue[idPart] ;
//
int leading[3] = {0,0,0 } ;
for (int k=1, t=1, tail=1; k <= P; ++k, tail=t)
{
......@@ -544,35 +532,7 @@ public:
for( i=0 ; i < SizeVector ; ++i)
{
multipole[i] *=_coeffPoly[i] ;
// printf(" M[%d] %f\n", i,multipole[i] );
}
// //
// // Nouveau parcours graded lexicographic order.
// // std::cout << std::endl;
// int leading[3] ;
// leading[0] = 0;
// leading[1] = 0;
// leading[2] = 0;
// std::cout << std::endl << " My ordering " <<std::endl;
// // Cas o = 0
// // o > 0
// for (int k=1, t=1, tail=1; k <= P; ++k, tail=t)
// {
// std::cout << "order: "<<k << " tail "<<tail <<std::endl;
// for (int i = 0; i < 3; ++i)
// {
// int head = leading[i];
// leading[i] = t;
// std::cout << " i " << i << " l "<<leading[i] << " tail " << tail<<std::endl;
// for ( int j = head; j < tail; ++j, ++t)
// { std::cout << " ( " << t<<" , " <<j<<" )" ; }
// std::cout <<std::endl;
// } // for i
// }// for k
}
}
......@@ -585,7 +545,7 @@ public:
const CellClass*const FRestrict *const FRestrict child,
const int inLevel)
{
// printf("M2M %%%%%%%%%%%%%%%%%%%%%%%%%%\n");
//Powers of expansions
int a=0,b=0,c=0;
......@@ -593,20 +553,14 @@ public:
int idx_a,idx_b,idx_c;
//Distance from current child to parent
// FReal boxSize = (boxWidth/FReal(1 << inLevel));
FReal dx = 0.0;
FReal dy = 0.0;
FReal dz = 0.0;
//Center point of parent cell
const FPoint& cellCenter = getCellCenter(pole->getCoordinate(),inLevel);
// printf("M2M :: fatherCell : X: %f, Y: %f, Z:%f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ());
FReal * FRestrict mult = pole->getMultipole();
// for(int idx=0 ; idx<SizeVector ; ++idx)
// {
// printf("Coeff Père avant mult[%d]: %f\n",idx,mult[idx]);
// }
// FMemUtils::memset(mult,FReal(0.0),SizeVector*FReal(0.0));
//Iteration over the eight children
int idxChild;
FReal coeff;
......@@ -615,15 +569,13 @@ public:
if(child[idxChild]){
//
const FPoint& childCenter = getCellCenter(child[idxChild]->getCoordinate(),inLevel+1);
// printf("M2M :: Father cell : (%f,%f,%f) child cells : X: %f, Y: %f, Z:%f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ(),childCenter.getX(),childCenter.getY(),childCenter.getZ());
const FReal * FRestrict multChild = child[idxChild]->getMultipole();
//Set the distance between centers of cells
dx = cellCenter.getX() - childCenter.getX();
dy = cellCenter.getY() - childCenter.getY();
dz = cellCenter.getZ() - childCenter.getZ();
// printf("M2M :: dx=%f, dy=%f, dz=%f\n",dx,dy,dz);
// Precompute the arrays of dx^i
arrayDX[0] = 1.0 ;
arrayDY[0] = 1.0 ;
......@@ -638,12 +590,11 @@ public:
b=0;
c=0;
//FReal value ;
FReal value;// Nk ;
FReal value;
//Iteration over parent multipole array
for(int idxMult = 0 ; idxMult<SizeVector ; idxMult++)
{
value = 0.0;
// Nk = factorials[a+b+c]/fact3int(a,b,c) ;
int idMultiChild;
//Iteration over the powers to find the cell multipole
//involved in the computation of the parent multipole
......@@ -654,29 +605,15 @@ public:
//Computation
//Child multipole involved
idMultiChild = powerToIdx(a-idx_a,b-idx_b,c-idx_c);
// coeff = fact3int(a-idx_a,b-idx_b,c-idx_c)/(factorials[a+b+c-idx_a-idx_b-idx_c]*fact3int(idx_a,idx_b,idx_c));
coeff = FReal(1.0)/fact3int(idx_a,idx_b,idx_c);
value += multChild[idMultiChild]*coeff*arrayDX[idx_a]*arrayDY[idx_b]*arrayDZ[idx_c];
}
}
}
// mult[idxMult] += Nk*value;
mult[idxMult] += value;
incPowers(&a,&b,&c);
}
//For debugging purposes
// int x=0,y=0,z=0;
// for(int idxChk=0 ; idxChk<SizeVector ; ++idxChk)
// {
// // printf("M2M in (%f,%f,%f) child (%f,%f,%f) MF %f and MP %f\n",cellCenter.getX(),cellCenter.getY(),cellCenter.getZ(),childCenter.getX(),childCenter.getY(),childCenter.getZ(), multChild[idxChk],mult[idxChk]);
// incPowers(&x,&y,&z);
// }
}
// Multiply by a coeff
for(int idxMult = 0 ; idxMult<SizeVector ; ++idxMult){
// mult[idxMult] *= this->_coeffNk[idxMult];
}
}
}
......@@ -700,34 +637,24 @@ public:
const CellClass* distantNeighbors[343], // Sources to be read
const int /*size*/, const int inLevel)
{
//printf("M2L %%%%%%%%%%%%%%%%%%%%%%%%%%\n");
//Iteration over distantNeighbors
//index for distant neighbors
int idxNeigh;
// int sizeDerivative = (2*P+1)*(P+1)*(2*P+3)/3;
const FPoint & locCenter = getCellCenter(local->getCoordinate(),inLevel);
FReal * FRestrict iterLocal = local->getLocal();
// FMemUtils::memset(iterLocal,0,SizeVector*sizeof(FReal(0.0)));
// FReal yetComputed[sizeDerivative];
for(idxNeigh=0 ; idxNeigh<343 ; ++idxNeigh){
//Need to test if current neighbor is one of the interaction list
if(distantNeighbors[idxNeigh]){
const FPoint curDistCenter = getCellCenter(distantNeighbors[idxNeigh]->getCoordinate(),inLevel);
//Derivatives are computed iteratively
//printf("M2L %%%%%%%%%%%%%%%% (%f,%f,%f) ----> (%f,%f,%f)\n",curDistCenter.getX(),curDistCenter.getY(),curDistCenter.getZ(),
// locCenter.getX(),locCenter.getY(),locCenter.getZ());
//Derivatives are computed iteratively
FMemUtils::memset(this->_PsiVector,0,sizeDerivative*sizeof(FReal(0.0)));
//
// Compute derivatives on locCenter - curDistCenter
// target source
FReal dx = locCenter.getX() - curDistCenter.getX();
FReal dy = locCenter.getY() - curDistCenter.getY();
FReal dz = locCenter.getZ() - curDistCenter.getZ();
......@@ -735,8 +662,8 @@ public:
//Computation of all the derivatives needed
this->computeFullDerivative(dx,dy,dz,this->_PsiVector);
const FReal * multipole = distantNeighbors[idxNeigh]->getMultipole();
//Iteration over Multipole / Local
int al=0,bl=0,cl=0; // For local array
int am,bm /*,cm*/; // For distant array
......@@ -760,23 +687,19 @@ public:
am=0; bm=0;
//Case (am,bm,cm) --> (0,0,0)
tmp += this->_PsiVector[i]*multipole[0];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",i,al,bl,cl,al,bl,cl,0,0,0,0);
for(int ord=1,id=1 ; ord<=P ; ++ord){
//Case (am,bm,cm) --> (ord,0,0)
int idxPsi = powerToIdx(ord+al,bl,cl);
tmp += this->_PsiVector[idxPsi]*multipole[id];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",idxPsi,ord+al,bl,cl,al,bl,cl,ord,0,0,id);
//Case (am,bm,cm) --> (ord-1,1,0)
idxPsi = powerToIdx((ord-1)+al,1+bl,cl);
tmp += this->_PsiVector[idxPsi]*multipole[++id];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",idxPsi,(ord-1)+al,1+bl,cl,al,bl,cl,ord-1,1,0,id);
//Case (am,bm,cm) --> (ord-1,0,1)
idxPsi = powerToIdx((ord-1)+al,bl,1+cl);
tmp += this->_PsiVector[idxPsi]*multipole[++id];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",idxPsi,(ord-1)+al,bl,1+cl,al,bl,cl,ord-1,0,1,id);
//End of specific case
//Start of general case :
......@@ -787,12 +710,10 @@ public:
{
idxPsi = powerToIdx(al+am,bl+bm,cl+(ord-(am+bm)));
tmp += this->_PsiVector[idxPsi]*multipole[id];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",idxPsi,al+am,bl+bm,cl+(ord-(am+bm)),al,bl,cl,am,bm,ord-am-bm,id);
}
}
}
iterLocal[i] += tmp*coeffL;
//printf("i==%d orn=%d, iterLocal[i]:: %f\n",i,orn,iterLocal[i]);
}
}
else
......@@ -808,22 +729,18 @@ public:
bm=0;
//Case (am,bm,cm) --> (0,0,0)
tmp += this->_PsiVector[i]*multipole[0];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",i,al,bl,cl,al,bl,cl,0,0,0,0);
for(int ord=1,id=1 ; ord<=P ; ++ord){
//Case (am,bm,cm) --> (ord,0,0)
int idxPsi = powerToIdx(ord+al,bl,cl);
tmp += this->_PsiVector[idxPsi]*multipole[id];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",idxPsi,ord+al,bl,cl,al,bl,cl,ord,0,0,id);
//Case (am,bm,cm) --> (ord-1,1,0)
idxPsi = powerToIdx((ord-1)+al,1+bl,cl);
tmp += this->_PsiVector[idxPsi]*multipole[++id];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",idxPsi,(ord-1)+al,1+bl,cl,al,bl,cl,ord-1,1,0,id);
//Case (am,bm,cm) --> (ord-1,0,1)
idxPsi = powerToIdx((ord-1)+al,bl,1+cl);
tmp += this->_PsiVector[idxPsi]*multipole[++id];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",idxPsi,(ord-1)+al,bl,1+cl,al,bl,cl,ord-1,0,1,id);
//End of specific case
//Start of general case :
......@@ -834,13 +751,10 @@ public:
{
idxPsi = powerToIdx(al+am,bl+bm,cl+(ord-(am+bm)));
tmp += this->_PsiVector[idxPsi]*multipole[id];
//printf("idxPsi: %d (%d,%d,%d) (%d,%d,%d) (%d,%d,%d) id==%d \n",idxPsi,al+am,bl+bm,cl+(ord-(am+bm)),al,bl,cl,am,bm,ord-am-bm,id);
}
}
}
iterLocal[i] += tmp*coeffL;
//printf("i==%d orn=%d, iterLocal[i]:: %f\n",i,orn,iterLocal[i]);
//printf("i==%d orn=%d\n",i,orn);
}
}
}
......@@ -868,7 +782,6 @@ public:
CellClass* FRestrict * const FRestrict childCell,
const int inLevel)
{
//printf("L2L %%%%%%%%%%%%%%%%%%%%%%%%%%\n");
FPoint locCenter = getCellCenter(fatherCell->getCoordinate(),inLevel);
// Get father local expansion
......@@ -885,8 +798,6 @@ public:
FReal* FRestrict childExpansion = childCell[idxChild]->getLocal() ;
const FPoint& childCenter = getCellCenter(childCell[idxChild]->getCoordinate(),inLevel+1);
//printf("L2L %%%%%%%%%%% (%f,%f,%f) ---> (%f,%f,%f) \n",locCenter.getX(),locCenter.getY(),locCenter.getZ(),childCenter.getX(),childCenter.getY(),childCenter.getZ());
//Set the distance between centers of cells
// Child - father
dx = childCenter.getX()-locCenter.getX();
......@@ -943,9 +854,7 @@ public:
void L2P(const CellClass* const local,
ContainerClass* const particles)
{
//printf("L2P %%%%%%%%%%%%%%%%%%%%%%%%%%\n");
FPoint locCenter = getLeafCenter(local->getCoordinate());
// printf("L2P :: locCenter (%f,%f,%f)\n",locCenter.getX(),locCenter.getY(),locCenter.getZ());
//Iterator over particles
int nbPart = particles->getNbParticles();
......@@ -953,10 +862,6 @@ public:
//Iteration over Local array
//
const FReal * iterLocal = local->getLocal();
for(int k=0 ; k<SizeVector ; ++k)
{
// printf("iterLocal[%d] :: %f\n",k,iterLocal[k]);
}
const FReal * const * positions = particles->getPositions();
const FReal * posX = positions[0];
const FReal * posY = positions[1];
......@@ -965,18 +870,18 @@ public:
FReal * forceX = particles->getForcesX();
FReal * forceY = particles->getForcesY();
FReal * forceZ = particles->getForcesZ();
//
FReal * targetsPotentials = particles->getPotentials();
FReal * phyValues = particles->getPhysicalValues();
//Iteration over particles
for(int i=0 ; i<nbPart ; ++i){
//
FReal dx = posX[i] - locCenter.getX();
FReal dy = posY[i] - locCenter.getY();
FReal dz = posZ[i] - locCenter.getZ();
// printf("Distance : (%f,%f,%f)\n",dx,dy,dz);
// Precompute an arrays of Array[i] = dx^(i-1)
arrayDX[0] = 0.0 ;
arrayDY[0] = 0.0 ;
......@@ -1010,9 +915,7 @@ public:
forceX[i] += partPhyValue*locForceX ;
forceY[i] += partPhyValue*locForceY ;
forceZ[i] += partPhyValue*locForceZ ;
// printf(" L2P part %d %f %f %f %f\n",i,targetsPotentials[i],forceX[i],forceY[i], forceZ[i]);
}
// printf(" L2P end\n");
}
/**
......@@ -1028,7 +931,6 @@ public:
ContainerClass* const FRestrict targets, const ContainerClass* const FRestrict /*sources*/,
ContainerClass* const directNeighborsParticles[27], const int /*size*/)
{
// printf("P2P %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n");
FP2P::FullMutual(targets,directNeighborsParticles,14);
}
......
......@@ -152,15 +152,15 @@ private:
int powerToIdx(const int a,const int b,const int c)
{
int t,res,p = a+b+c;
if(p==0) // si p=0 alors res=1/3=0 donc on peut l'enlever
{return 0;}
else
{
res = p*(p+1)*(p+2)/6;
t = p - a;
res += t*(t+1)/2+c;
return res;
}
// if(p==0) // si p=0 alors res=1/3=0 donc on peut l'enlever
// {return 0;}
// else
// {
res = p*(p+1)*(p+2)/6;
t = p - a;
res += t*(t+1)/2+c;
return res;
// }
}
/* Return the factorial of a number
......@@ -415,7 +415,8 @@ private:
}
for(c=2 ; c <= 2*P-b-a ; ++c){
//Computation of derivatives Psi_{a,b,c} with a >= 2
// |x-y|^2*Psi_{a,b,c} + (2*a-1)*dx*Psi_{a-1,b,c} + a*(a-2)*Psi_{a-2,b,c} + 2*b*dy*Psi_{a,b-1,c} + b*(b-1)*Psi_{a,b-2,c} + 2*c= 0
// |x-y|^2*Psi_{a,b,c} + (2*a-1)*dx*Psi_{a-1,b,c} + a*(a-2)*Psi_{a-2,b,c} + 2*b*dy*Psi_{a,b-1,c} + b*(b-1)*Psi_{a,b-2,c}
// + 2*c*dz*Psi_{a,b,c-1}} = 0
idxTarget = powerToIdx(a,b,c);
idxSrc1 = powerToIdx(a-1,b,c);
idxSrc2 = powerToIdx(a,b-1,c);
......@@ -749,21 +750,15 @@ public:
//Iterating over local array : n
for(int i=0 ; i< SizeVector ; ++i){
FReal fctl = fact3int(al,bl,cl);
//FReal coeffL = factorials[al+bl+cl]/(fctl*fctl);
FReal coeffL = FReal(1.0)/(fctl);
//Iterator over multipole array
FReal tmp = 0.0 ;
//Iterating over multipole array : k
// Loc(al,bl,cl) = N(al,bl,cl)/((al,bl,cl)!*(al,bl,cl)!) sum_(am,bm,cm) Psi[am+al,bm+bl,cm+cl] * M[am,bm,cm]
//
am=0; bm=0; cm=0;
//printf("al= %d, bl=%d, cl=%d ==> i =%d \n",al,bl,cl,i);
for(int j = 0 ; j < SizeVector ; ++j){ //corresponding powers am,bm,cm
int idxPsi = powerToIdx(al+am,bl+bm,cl+cm);
tmp += this->_PsiVector[idxPsi]*multipole[j];
// if((al+bl+cl)==0) {printf("coeffL:%f (%d,%d,%d)--> mult x psi ::> %f * %f tmp: %f\n",coeffL,am,bm,cm,multipole[j],_PsiVector[idxPsi],tmp);}
//updating a,b,c
incPowers(&am,&bm,&cm);
}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment