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 a036c5dd authored by COULAUD Olivier's avatar COULAUD Olivier

Merge commit '87a868c8'

* commit '87a868c8':
  Making ffmper.sh a bit more maintainable. No need to list all the numbers of the interval you want to use.
  Add test for adaptive Uniform kernel and update FUnifCell.
  Provided copy ctor for all Uniform M2L handlers.
parents f06cfcd1 87a868c8
#ifndef FADAPTUNIFKERNEL_HPP
#define FADAPTUNIFKERNEL_HPP
// ===================================================================================
// Copyright ScalFmm 2011 INRIA,
// 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.
//
// 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.gnu.org/licenses".
// ===================================================================================
#include "Utils/FGlobal.hpp"
#include "Utils/FTrace.hpp"
#include "Utils/FPoint.hpp"
#include "Adaptative/FAdaptiveCell.hpp"
#include "Adaptative/FAdaptiveKernelWrapper.hpp"
#include "Adaptative/FAbstractAdaptiveKernel.hpp"
#include "Kernels/Uniform/FUnifKernel.hpp"
#include "Kernels/Uniform/FUnifM2LHandler.hpp"
class FTreeCoordinate;
// ==== CMAKE =====
// @FUSE_BLAS
// ================
// for verbosity only!!!
//#define COUNT_BLOCKED_INTERACTIONS
// if timings should be logged
//#define LOG_TIMINGS
/**
* @author O. Coulaud
* @class FAdaptUnifKernel
* @brief
* Please read the license
*
* This kernels implement the Lagrange interpolation based FMM operators.
* It implements all interfaces (P2P, P2M, M2M, M2L, L2L, L2P) which are
* required by the FFmmAlgorithm and FFmmAlgorithmThread.
*
* @tparam CellClass Type of cell
* @tparam ContainerClass Type of container to store particles
* @tparam MatrixKernelClass Type of matrix kernel function
* @tparam ORDER Chebyshev interpolation order
*/
template< class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FAdaptiveUnifKernel : FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
, public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
//
typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
/// Needed for M2L operator
typedef FUnifM2LHandler<ORDER,MatrixKernelClass::Type> M2LHandlerClass;
const M2LHandlerClass M2LHandler;
public:
using KernelBaseClass::P2M;
using KernelBaseClass::M2M;
using KernelBaseClass::M2L;
using KernelBaseClass::finishedLevelM2L;
using KernelBaseClass::L2L;
using KernelBaseClass::L2P;
using KernelBaseClass::P2P;
using KernelBaseClass::P2PRemote;
// /**
// * The constructor initializes all constant attributes and it reads the
// * precomputed and compressed M2L operators from a binary file (an
// * runtime_error is thrown if the required file is not valid).
// */
FAdaptiveUnifKernel(const int inTreeHeight, const FReal inBoxWidth,
const FPoint& inBoxCenter) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter), M2LHandler(KernelBaseClass::MatrixKernel.getPtr(), inTreeHeight, inBoxWidth)
{}
// /** Copy constructor */
FAdaptiveUnifKernel(const FAdaptiveUnifKernel& other)
: KernelBaseClass(other), M2LHandler(other.M2LHandler)
{ }
//
// /** Destructor */
~FAdaptiveUnifKernel()
{
//this->~KernelBaseClass() ;
}
void P2M(CellClass* const pole, const int cellLevel, const ContainerClass* const particles) override {
//pole->setDataUp(pole->getDataUp() + particles->getNbParticles());
//
const FPoint LeafCellCenter(KernelBaseClass::getLeafCellCenter(pole->getCoordinate()));
const FReal BoxWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-cellLevel);
//
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
KernelBaseClass::Interpolator->applyP2M(LeafCellCenter, BoxWidth,
pole->getMultipole(idxRhs), particles);
// 2) apply Discrete Fourier Transform
M2LHandler.applyZeroPaddingAndDFT(pole->getMultipole(idxRhs),
pole->getTransformedMultipole(idxRhs));
}
}
void M2M(CellClass* const pole, const int /*poleLevel*/, const CellClass* const subCell, const int /*subCellLevel*/) override {
// pole->setDataUp(pole->getDataUp() + subCell->getDataUp());
}
void P2L(CellClass* const local, const int /*localLevel*/, const ContainerClass* const particles) override {
// local->setDataDown(local->getDataDown() + particles->getNbParticles());
}
void M2L(CellClass* const local, const int /*localLevel*/, const CellClass* const aNeighbor, const int /*neighborLevel*/) override {
// local->setDataDown(local->getDataDown() + aNeighbor->getDataUp());
}
void M2P(const CellClass* const pole, const int /*poleLevel*/, ContainerClass* const particles) override {
// long long int*const particlesAttributes = particles->getDataDown();
// for(int idxPart = 0 ; idxPart < particles->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += pole->getDataUp();
// }
}
void L2L(const CellClass* const local, const int /*localLevel*/, CellClass* const subCell, const int /*subCellLevel*/) override {
// subCell->setDataDown(local->getDataDown() + subCell->getDataDown());
}
void L2P(const CellClass* const local, const int /*cellLevel*/, ContainerClass* const particles) override {
// long long int*const particlesAttributes = particles->getDataDown();
// for(int idxPart = 0 ; idxPart < particles->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += local->getDataDown();
// }
}
void P2P(ContainerClass* target, const ContainerClass* sources) override {
// long long int*const particlesAttributes = target->getDataDown();
// for(int idxPart = 0 ; idxPart < target->getNbParticles() ; ++idxPart){
// particlesAttributes[idxPart] += sources->getNbParticles();
// }
}
bool preferP2M(const ContainerClass* const particles) override {
return particles->getNbParticles() < 10;
}
bool preferP2M(const int /*atLevel*/, const ContainerClass*const particles[], const int nbContainers) override {
int counterParticles = 0;
for(int idxContainer = 0 ; idxContainer < nbContainers ; ++idxContainer){
counterParticles += particles[idxContainer]->getNbParticles();
}
return counterParticles < 10;
}
};
//
//template < class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
//class FAdaptUnifKernel
// : public FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
//{
// typedef FUnifKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
//
//#ifdef LOG_TIMINGS
// FTic time;
// FReal t_m2l_1, t_m2l_2, t_m2l_3;
//#endif
//
//public:
// /**
// * The constructor initializes all constant attributes and it reads the
// * precomputed and compressed M2L operators from a binary file (an
// * runtime_error is thrown if the required file is not valid).
// */
// FAdaptUnifKernel(const int inTreeHeight,
// const FReal inBoxWidth,
// const FPoint& inBoxCenter)
//: KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter)
//{
//
//#ifdef LOG_TIMINGS
// t_m2l_1 = FReal(0.);
// t_m2l_2 = FReal(0.);
// t_m2l_3 = FReal(0.);
//#endif
//}
//
//
// /** Copy constructor */
// FAdaptUnifKernel(const FAdaptUnifKernel& other)
// : KernelBaseClass(other)
// { }
//
//
//
// /** Destructor */
// ~FAdaptUnifKernel()
// {
// this->~KernelBaseClass() ;
//#ifdef LOG_TIMINGS
// std::cout << "- Permutation took " << t_m2l_1 << "s"
// << "\n- GEMMT and GEMM took " << t_m2l_2 << "s"
// << "\n- Unpermutation took " << t_m2l_3 << "s"
// << std::endl;
//#endif
// }
//
//
// void P2MAdapt(CellClass* const ParentCell, const int &level)
// {
// const FPoint LeafCellCenter(KernelBaseClass::getLeafCellCenter(ParentCell->getCoordinate()));
// const FReal BoxWidth = KernelBaseClass::BoxWidthLeaf*FMath::pow(2.0,KernelBaseClass::TreeHeight-level);
// //
// for(int i = 0 ; i <ParentCell->getLeavesSize(); ++i ){
// //
// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// KernelBaseClass::Interpolator->applyP2M(LeafCellCenter, BoxWidth,
// ParentCell->getMultipole(idxRhs), ParentCell->getLeaf(i)->getSrc());
// }
// }
// }
// void M2MAdapt(CellClass* const FRestrict ParentCell, const int &TreeLevel, const int &numberOfM2M,
// const int * FRestrict ChildLevel , const CellClass*const FRestrict *const FRestrict ChildCells)
// {
// for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // // apply Sy
// for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
// if (ChildCells[ChildIndex]){
// // KernelBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs));
// }
// }
// }
// }
//
//
//
// void M2L(CellClass* const FRestrict TargetCell,
// const CellClass* SourceCells[343],
// const int /*NumSourceCells*/,
// const int TreeLevel)
// {
//
// }
//
//
// void L2L(const CellClass* const FRestrict ParentCell,
// CellClass* FRestrict *const FRestrict ChildCells,
// const int /*TreeLevel*/)
// {
// // for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// // // apply Sx
// // for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
// // if (ChildCells[ChildIndex]){
// // AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxRhs), ChildCells[ChildIndex]->getLocal(idxRhs));
// // }
// // }
// // }
// }
//
// void L2P(const CellClass* const LeafCell,
// ContainerClass* const TargetParticles)
// {
// KernelBaseClass::L2P(LeafCell,TargetParticles) ;
// }
//
// // void P2P(const FTreeCoordinate& /* LeafCellCoordinate */, // needed for periodic boundary conditions
// // ContainerClass* const FRestrict TargetParticles,
// // const ContainerClass* const FRestrict /*SourceParticles*/,
// // ContainerClass* const NeighborSourceParticles[27],
// // const int /* size */)
// // {
// // DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
// // }
// //
// //
// // void P2PRemote(const FTreeCoordinate& /*inPosition*/,
// // ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
// // ContainerClass* const inNeighbors[27], const int /*inSize*/){
// // DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
// // }
//
//};
//
//
#endif //FADAPTUNIFKERNELS_HPP
// [--END--]
......@@ -170,6 +170,21 @@ public:
return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplexe);
}
template <class StreamClass>
friend StreamClass& operator<<(StreamClass& output, const FUnifCell<ORDER, NRHS, NLHS, NVALS>& cell){
output <<" Multipole exp NRHS " << NRHS <<" NVALS " <<NVALS << " VectorSize " << cell.getVectorSize() << std::endl;
for (int rhs= 0 ; rhs < NRHS ; ++rhs) {
const FReal* pole = cell.getMultipole(rhs);
for (int val= 0 ; val < NVALS ; ++val) {
output<< " val : " << val << " exp: " ;
for (int i= 0 ; i < cell.getVectorSize() ; ++i) {
output<< pole[i] << " ";
}
output << std::endl;
}
}
return output;
}
};
......
......@@ -211,6 +211,20 @@ public:
}
/*
* Copy constructor
*/
FUnifM2LHandler(const FUnifM2LHandler& other)
: FC(other.FC), Dft(), opt_rc(other.opt_rc)
{
// init DFT
const int steps[dimfft] = {rc};
Dft.buildDFT(steps);
// copy node_diff
memcpy(node_diff,other.node_diff,sizeof(unsigned int)*nnodes*nnodes);
}
~FUnifM2LHandler()
{ }
......@@ -350,8 +364,8 @@ public:
FUnifM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth)
: TreeHeight(inTreeHeight),
RootCellWidth(inRootCellWidth),
opt_rc(rc/2+1),
Dft()
Dft(), opt_rc(rc/2+1)
{
// init DFT
const int steps[dimfft] = {rc};
......@@ -370,6 +384,25 @@ public:
}
/*
* Copy constructor
*/
FUnifM2LHandler(const FUnifM2LHandler& other)
: FC(other.FC),
TreeHeight(other.TreeHeight),
RootCellWidth(other.RootCellWidth),
Dft(), opt_rc(other.opt_rc)
{
// init DFT
const int steps[dimfft] = {rc};
Dft.buildDFT(steps);
// copy node_diff
memcpy(node_diff,other.node_diff,sizeof(unsigned int)*nnodes*nnodes);
}
~FUnifM2LHandler()
{
for (unsigned int l=0; l<TreeHeight; ++l)
......
......@@ -184,7 +184,7 @@ static void Compute(const MatrixKernelClass *const MatrixKernel,
template <int ORDER, class MatrixKernelClass, KERNEL_FUNCTION_TYPE TYPE> class FUnifTensorialM2LHandler;
template <int ORDER, class MatrixKernelClass>
class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyable
class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,HOMOGENEOUS>
{
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes,
......@@ -239,6 +239,22 @@ public:
ComputeAndSet(MatrixKernel);
}
/*
* Copy constructor
*/
FUnifTensorialM2LHandler(const FUnifTensorialM2LHandler& other)
: FC(other.FC),
CellWidthExtension(other.CellWidthExtension),
Dft(), opt_rc(other.opt_rc)
{
// init DFT
const int steps[dimfft] = {rc};
Dft.buildDFT(steps);
// copy node_diff
memcpy(node_diff,other.node_diff,sizeof(unsigned int)*nnodes*nnodes);
}
~FUnifTensorialM2LHandler()
{
for (unsigned int d=0; d<ncmp; ++d)
......@@ -358,7 +374,7 @@ public:
template <int ORDER, class MatrixKernelClass>
class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCopyable
class FUnifTensorialM2LHandler<ORDER,MatrixKernelClass,NON_HOMOGENEOUS>
{
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes,
......@@ -422,6 +438,25 @@ public:
ComputeAndSet(MatrixKernel);
}
/*
* Copy constructor
*/
FUnifTensorialM2LHandler(const FUnifTensorialM2LHandler& other)
: FC(other.FC),
TreeHeight(other.TreeHeight),
RootCellWidth(other.RootCellWidth),
CellWidthExtension(other.CellWidthExtension),
Dft(), opt_rc(other.opt_rc)
{
// init DFT
const int steps[dimfft] = {rc};
Dft.buildDFT(steps);
// copy node_diff
memcpy(node_diff,other.node_diff,sizeof(unsigned int)*nnodes*nnodes);
}
~FUnifTensorialM2LHandler()
{
for (unsigned int l=0; l<TreeHeight; ++l)
......
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// 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.
//
// 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.gnu.org/licenses".
// ===================================================================================
// ==== CMAKE =====
// @FUSE_BLAS
// ================
#include <iostream>
#include <cstdio>
#include "Utils/FParameters.hpp"
#include "Utils/FTic.hpp"
#include "Containers/FOctree.hpp"
//#include "Containers/FVector.hpp"
//#include "Components/FSimpleLeaf.hpp"
#include "Utils/FPoint.hpp"
#include "Files/FFmaGenericLoader.hpp"
#include "Files/FRandomLoader.hpp"
#include "Components/FBasicKernels.hpp"
#include "Components/FSimpleIndexedLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Adaptative/FAdaptiveCell.hpp"
#include "Adaptative/FAdaptiveKernelWrapper.hpp"
#include "Adaptative/FAbstractAdaptiveKernel.hpp"
//
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Uniform/FUnifCell.hpp"
#include "AdaptiveTree/FAdaptUnifKernel.hpp"
#include "AdaptiveTree/FAdaptTools.hpp"
//
//
#include "Core/FFmmAlgorithm.hpp"
//#include "Core/FFmmAlgorithmThread.hpp"
//#include "Core/FFmmAlgorithmTask.hpp"
/** This program show an example of use of the fmm basic algo
* it also check that each particles is impacted each other particles
*/
void usage() {
std::cout << "Driver to obtain statistics on the octree" << std::endl;
std::cout << "Options "<< std::endl
<< " -help to see the parameters " << std::endl
<< " -depth the depth of the octree "<< std::endl
<< " -subdepth specifies the size of the sub octree " << std::endl
<< " -fin name specifies the name of the particle distribution" << std::endl
<< " -sM s_min^M threshold for Multipole (l+1)^2 for Spherical harmonics"<<std::endl
<< " -sL s_min^L threshold for Local (l+1)^2 for Spherical harmonics"<<std::endl;
}
// Simply create particles and try the kernels
int main(int argc, char ** argv){
//
// accuracy
const unsigned int P = 3 ;
// typedef FTestCell CellClass;
// typedef FAdaptiveTestKernel< CellClass, ContainerClass > KernelClass;
typedef FUnifCell<P> CellClass;
typedef FP2PParticleContainerIndexed<> ContainerClass;
typedef FSimpleIndexedLeaf<ContainerClass> LeafClass;
typedef FInterpMatrixKernelR MatrixKernelClass;
//
typedef FAdaptiveUnifKernel<CellClass,ContainerClass,MatrixKernelClass,P> KernelClass;
//
//
typedef FAdaptiveCell< CellClass, ContainerClass > CellWrapperClass;
typedef FAdaptiveKernelWrapper< KernelClass, CellClass, ContainerClass > KernelWrapperClass;
typedef FOctree< CellWrapperClass, ContainerClass , LeafClass > OctreeClass;
// FFmmAlgorithmTask FFmmAlgorithmThread
typedef FFmmAlgorithm<OctreeClass, CellWrapperClass, ContainerClass, KernelWrapperClass, LeafClass > FmmClass;
///////////////////////What we do/////////////////////////////
std::cout << ">> This executable has to be used to test the FMM algorithm.\n";
//////////////////////////////////////////////////////////////
//
const int NbLevels = FParameters::getValue(argc,argv,"-depth", 7);
const int SizeSubLevels = FParameters::getValue(argc,argv,"subdepth", 3);
const int sminM = FParameters::getValue(argc,argv,"-sM", P*P*P);
const int sminL = FParameters::getValue(argc,argv,"-sL", P*P*P);
//
FTic counter;
//////////////////////////////////////////////////////////////////////////////////
// Not Random Loader
//////////////////////////////////////////////////////////////////////////////////
const std::string fileName(FParameters::getStr(argc,argv,"-fin", "../Data/prolate50.fma"));
FFmaGenericLoader loader(fileName);
const long int NbPart = loader.getNumberOfParticles() ;
// Random Loader
//const int NbPart = FParameters::getValue(argc,argv,"-nb", 2000000);
// FRandomLoader loader(NbPart, 1, FPoint(0.5,0.5,0.5), 1);
//////////////////////////////////////////////////////////////////////////////////
OctreeClass tree(NbLevels, SizeSubLevels, loader.getBoxWidth(), loader.getCenterOfBox());
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Creating & Inserting " << NbPart << " particles ..." << std::endl;
std::cout << "\tHeight : " << NbLevels << " \t sub-height : " << SizeSubLevels << std::endl;
std::cout << " criteria SM: "<< sminM <<std::endl
<< " criteria SL: "<< sminL <<std::endl <<std::endl;
//
{
counter.tic();
FReal L= loader.getBoxWidth();
FmaRParticle* particles= new FmaRParticle[NbPart];
FPoint minPos(L,L,L), maxPos(-L,-L,-L);
//
loader.fillParticle(particles,NbPart);
for(int idxPart = 0 ; idxPart < NbPart; ++idxPart){
const FPoint PP(particles[idxPart].getPosition() ) ;
//
minPos.setX(FMath::Min(minPos.getX(),PP.getX())) ;
minPos.setY(FMath::Min(minPos.getY(),PP.getY())) ;
minPos.setZ(FMath::Min(minPos.getZ(),PP.getZ())) ;
maxPos.setX(FMath::Max(maxPos.getX(),PP.getX())) ;
maxPos.setY(FMath::Max(maxPos.getY(),PP.getY())) ;
maxPos.setZ(FMath::Max(maxPos.getZ(),PP.getZ())) ;
//
tree.insert(PP, idxPart, particles[idxPart].getPhysicalValue());
}
counter.tac();
std::cout << "Data are inside the box delimited by "<<std::endl
<< " Min corner: "<< minPos<<std::endl
<< " Max corner: "<< maxPos<<std::endl <<std::endl;
std::cout << "Done " << "(@Creating and Inserting Particles = " << counter.elapsed() << " s)." << std::endl;
}
//////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////
std::cout << "Working on particles ..." << std::endl;
counter.tic();
KernelWrapperClass kernels(NbLevels, loader.getBoxWidth(), loader.getCenterOfBox());; // FTestKernels FBasicKernels
FmmClass algo(&tree,&kernels); //FFmmAlgorithm FFmmAlgorithmThread
algo.execute();
counter.tac();