Skip to content

GitLab

Commit cbdbdcdf authored by BLANCHARD Pierre's avatar BLANCHARD Pierre
Browse files

Implemented UniformGrid/Lagrange interpolation and factorized interpolation... 

Implemented UniformGrid/Lagrange interpolation and factorized interpolation specific elements such as the MatrixKernel, Tensor and Mapping classes + Added required tests.
parent 96d44b04
Hide whitespace changes
Inline Side-by-side
CMakeLists.txt
View file @ cbdbdcdf
......@@ -42,6 +42,7 @@ MESSAGE(STATUS " CXX ${CMAKE_CXX_COMPILER_ID}" )
#
# Options
OPTION( ScalFMM_USE_BLAS "Set to ON to build ScaFMM with BLAS" OFF )
OPTION( ScalFMM_USE_FFTW "Set to ON to build ScaFMM with FFTW" OFF )
OPTION( ScalFMM_USE_TRACE "Set to ON to print trace or use itac trace" OFF )
OPTION( ScalFMM_BUILD_TESTS "Set to ON to build fonctionnalities Tests" OFF )
OPTION( ScalFMM_BUILD_UTESTS "Set to ON to build UTests" OFF )
......@@ -138,7 +139,8 @@ endif()
if( ScalFMM_USE_BLAS )
OPTION( ScalFMM_USE_MKL_AS_BLAS "Set to ON to use MKL CBLAS" OFF )
if( ScalFMM_USE_MKL_AS_BLAS )
SET(BLAS_LIBRARIES "-L$ENV{MKLROOT}/lib;-lmkl_intel_lp64;-lmkl_sequential;-lmkl_core" CACHE STRING "Set your MKL flags")
SET(BLAS_LIBRARIES
"-L$ENV{MKLROOT}/lib;-lmkl_intel_lp64;-lmkl_sequential;-lmkl_core" CACHE STRING "Set your MKL flags")
SET(LAPACK_LIBRARIES "")
elseif(ScalFMM_USE_EXTERNAL_BLAS)
MESSAGE(STATUS "BLAS SET BY EXTERNAL PROGRAM = ${BLAS_LIBRARIES}")
......@@ -150,6 +152,18 @@ if( ScalFMM_USE_BLAS )
MESSAGE(STATUS "SCALFMM_LIBRARIES = ${SCALFMM_LIBRARIES}")
endif(ScalFMM_USE_BLAS)
# FFTW option
if( ScalFMM_USE_FFTW )
OPTION( ScalFMM_USE_MKL_AS_FFTW "Set to ON to use MKL FFTW" ON )
if( ScalFMM_USE_MKL_AS_FFTW )
SET(FFTW_LIBRARIES "-I$ENV{MKLROOT}/include/fftw; -L$ENV{MKLROOT}/lib; -lmkl_intel_lp64; -lmkl_sequential; -lmkl_core; -lpthread; -lm" CACHE STRING "Set your MKL flags")
else()
SET(FFTW_LIBRARIES "-lfftw3" CACHE STRING "Use LIBFFTW")
endif()
SET(SCALFMM_LIBRARIES "${SCALFMM_LIBRARIES}; ${FFTW_LIBRARIES}")
MESSAGE(STATUS "SCALFMM_LIBRARIES = ${SCALFMM_LIBRARIES}")
endif(ScalFMM_USE_FFTW)
# Compile option
#ADD_DEFINITIONS(-Wall -Wshadow -Wpointer-arith -Wcast-qual -Wconversion -fpic )
#
......
Src/Kernels/Chebyshev/FAbstractChebKernel.hpp
View file @ cbdbdcdf
......@@ -22,7 +22,7 @@
#include "../../Components/FAbstractKernels.hpp"
#include "FChebP2PKernels.hpp"
#include "../Interpolation/FInterpP2PKernels.hpp"
#include "./FChebInterpolator.hpp"
#include "../../Containers/FTreeCoordinate.hpp"
......
Src/Kernels/Chebyshev/FChebInterpolator.hpp
View file @ cbdbdcdf
......@@ -17,7 +17,7 @@
#define FCHEBINTERPOLATOR_HPP
#include "./FChebMapping.hpp"
#include "./../Interpolation/FInterpMapping.hpp"
#include "./FChebTensor.hpp"
#include "./FChebRoots.hpp"
......@@ -296,7 +296,7 @@ class FChebInterpolator : FNoCopyable
// set child info
FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter);
FChebTensor<ORDER>::setChebyshevRoots(ChildCenter, ChildWidth, ChildCoords);
FChebTensor<ORDER>::setPolynomialsRoots(ChildCenter, ChildWidth, ChildCoords);
// allocate memory
ChildParentInterpolator[child] = new FReal [3 * ORDER*ORDER];
......
Src/Kernels/Chebyshev/FChebKernel.hpp
View file @ cbdbdcdf
......@@ -70,8 +70,8 @@ public:
M2LHandler(new M2LHandlerClass(Epsilon))
{
// read precomputed compressed m2l operators from binary file
M2LHandler->ReadFromBinaryFileAndSet();
//M2LHandler->ComputeAndCompressAndSet();
//M2LHandler->ReadFromBinaryFileAndSet();
M2LHandler->ComputeAndCompressAndSet();
}
......@@ -209,14 +209,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
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*/){
DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
}
};
......
Src/Kernels/Chebyshev/FChebM2LHandler.hpp
View file @ cbdbdcdf
......@@ -99,8 +99,12 @@ public:
// check if aready set
if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
rank = ComputeAndCompress(epsilon, U, C, B);
unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);
// write info
std::cout << "Compressed and set M2L operators (" << rank << ") in "
std::cout << "Compressed and set M2L operators (" << long(sizeM2L) << " B) in "
<< time.tacAndElapsed() << "sec." << std::endl;
}
......
Src/Kernels/Chebyshev/FChebMatrixKernel.hpp deleted 100755 → 0
View file @ 96d44b04
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger 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".
// ===================================================================================
#ifndef FCHEBMATRIXKERNEL_HPP
#define FCHEBMATRIXKERNEL_HPP
#include "../../Utils/FPoint.hpp"
#include "../../Utils/FNoCopyable.hpp"
#include "../../Utils/FMath.hpp"
#include "../../Utils/FBlas.hpp"
// extendable
enum KERNEL_FUNCCTION_IDENTIFIER {ONE_OVER_R,
ONE_OVER_R_SQUARED,
LEONARD_JONES_POTENTIAL};
// probably not extedable :)
enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
/**
* @author Matthias Messner (matthias.messner@inria.fr)
* @class FChebMatrixKernels
* Please read the license
*/
struct FChebAbstractMatrixKernel : FNoCopyable
{
virtual ~FChebAbstractMatrixKernel(){} // to remove warning
virtual FReal evaluate(const FPoint&, const FPoint&) const = 0;
// I need both functions because required arguments are not always given
virtual FReal getScaleFactor(const FReal, const int) const = 0;
virtual FReal getScaleFactor(const FReal) const = 0;
};
/// One over r
struct FChebMatrixKernelR : FChebAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCCTION_IDENTIFIER Identifier = ONE_OVER_R;
FChebMatrixKernelR() {}
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
return FReal(1.) / FMath::Sqrt(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ());
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
return getScaleFactor(CellWidth);
}
FReal getScaleFactor(const FReal CellWidth) const
{
return FReal(2.) / CellWidth;
}
};
/// One over r^2
struct FChebMatrixKernelRR : FChebAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCCTION_IDENTIFIER Identifier = ONE_OVER_R_SQUARED;
FChebMatrixKernelRR() {}
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
return FReal(1.) / FReal(xy.getX()*xy.getX() +
xy.getY()*xy.getY() +
xy.getZ()*xy.getZ());
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
return getScaleFactor(CellWidth);
}
FReal getScaleFactor(const FReal CellWidth) const
{
return FReal(4.) / CellWidth;
}
};
/// One over r^12 - One over r^6
struct FChebMatrixKernelLJ : FChebAbstractMatrixKernel
{
static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
static const KERNEL_FUNCCTION_IDENTIFIER Identifier = LEONARD_JONES_POTENTIAL;
FChebMatrixKernelLJ() {}
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
const FReal r = xy.norm();
const FReal r3 = r*r*r;
const FReal one_over_r6 = FReal(1.) / (r3*r3);
//return one_over_r6 * one_over_r6;
//return one_over_r6;
return one_over_r6 * one_over_r6 - one_over_r6;
}
FReal getScaleFactor(const FReal, const int) const
{
// return 1 because non homogeneous kernel functions cannot be scaled!!!
return FReal(1.);
}
FReal getScaleFactor(const FReal) const
{
// return 1 because non homogeneous kernel functions cannot be scaled!!!
return FReal(1.);
}
};
/*! Functor which provides the interface to assemble a matrix based on the
number of rows and cols and on the coordinates x and y and the type of the
generating matrix-kernel function.
*/
template <typename MatrixKernel>
class EntryComputer
{
const MatrixKernel Kernel;
const unsigned int nx, ny;
const FPoint *const px, *const py;
const FReal *const weights;
public:
explicit EntryComputer(const unsigned int _nx, const FPoint *const _px,
const unsigned int _ny, const FPoint *const _py,
const FReal *const _weights = NULL)
: Kernel(), nx(_nx), ny(_ny), px(_px), py(_py), weights(_weights) {}
// template <typename Point>
// void operator()(const unsigned int nx, const Point *const px,
// const unsigned int ny, const Point *const py,
// FReal *const data) const
// {
// for (unsigned int j=0; j<ny; ++j)
// for (unsigned int i=0; i<nx; ++i)
// data[j*nx + i] = Kernel.evaluate(px[i], py[j]);
// }
void operator()(const unsigned int xbeg, const unsigned int xend,
const unsigned int ybeg, const unsigned int yend,
FReal *const data) const
{
unsigned int idx = 0;
if (weights) {
for (unsigned int j=ybeg; j<yend; ++j)
for (unsigned int i=xbeg; i<xend; ++i)
data[idx++] = weights[i] * weights[j] * Kernel.evaluate(px[i], py[j]);
} else {
for (unsigned int j=ybeg; j<yend; ++j)
for (unsigned int i=xbeg; i<xend; ++i)
data[idx++] = Kernel.evaluate(px[i], py[j]);
}
/*
// apply weighting matrices
if (weights) {
if ((xend-xbeg) == (yend-ybeg) && (xend-xbeg) == nx)
for (unsigned int n=0; n<nx; ++n) {
FBlas::scal(nx, weights[n], data + n, nx); // scale rows
FBlas::scal(nx, weights[n], data + n * nx); // scale cols
}
else if ((xend-xbeg) == 1 && (yend-ybeg) == ny)
for (unsigned int j=0; j<ny; ++j) data[j] *= weights[j];
else if ((yend-ybeg) == 1 && (xend-xbeg) == nx)
for (unsigned int i=0; i<nx; ++i) data[i] *= weights[i];
}
*/
}
};
#endif // FCHEBMATRIXKERNEL_HPP
// [--END--]
Src/Kernels/Chebyshev/FChebP2PKernels.hpp deleted 100644 → 0
View file @ 96d44b04
#ifndef FCHEBP2PKERNELS_HPP
#define FCHEBP2PKERNELS_HPP
#include "../P2P/FP2P.hpp"
template <KERNEL_FUNCCTION_IDENTIFIER Identifier, int NVALS>
struct DirectInteactionComputer;
///////////////////////////////////////////////////////
// P2P Wrappers
///////////////////////////////////////////////////////
/*! Specialization for Laplace potential */
template <>
struct DirectInteactionComputer<ONE_OVER_R, 1>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
}
};
/*! Specialization for Leonard-Jones potential */
template <>
struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL, 1>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutualLJ(TargetParticles,NeighborSourceParticles,14);
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemoteLJ(inTargets,inNeighbors,inSize);
}
};
///////////////////////////////////////////////////////
// In case of multi right hand side
///////////////////////////////////////////////////////
template <int NVALS>
struct DirectInteactionComputer<ONE_OVER_R, NVALS>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
}
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
}
}
};
/*! Specialization for Leonard-Jones potential */
template <int NVALS>
struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL, NVALS>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutualLJ(TargetParticles,NeighborSourceParticles,14);
}
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullRemoteLJ(inTargets,inNeighbors,inSize);
}
}
};
#endif // FCHEBP2PKERNELS_HPP
Src/Kernels/Chebyshev/FChebSymKernel.hpp
View file @ cbdbdcdf
......@@ -454,14 +454,14 @@ public:
ContainerClass* const NeighborSourceParticles[27],
const int /* size */)
{
DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
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*/){
DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
}
};
......
Src/Kernels/Chebyshev/FChebSymM2LHandler.hpp
View file @ cbdbdcdf
......@@ -479,9 +479,12 @@ static void precompute(const MatrixKernelClass *const MatrixKernel, const FReal
}
}
}
// std::cout << "The approximation of the " << counter
// << " far-field interactions (overall rank " << overall_rank << ") took "
// << overall_time << "s\n" << std::endl;
std::cout << "The approximation of the " << counter
<< " far-field interactions (overall rank " << overall_rank
<< " / " << 16*nnodes
<< " , sizeM2L= " << 2*overall_rank*nnodes*sizeof(FReal) << ""
<< " / " << 16*nnodes*nnodes*sizeof(FReal) << " B"
<< ") took " << overall_time << "s\n" << std::endl;
delete [] U;
delete [] WORK;
delete [] VT;
......
Src/Kernels/Chebyshev/FChebTensor.hpp
View file @ cbdbdcdf
......@@ -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 FCHEBTENSOR_HPP
......@@ -19,7 +19,7 @@
#include "../../Utils/FMath.hpp"
#include "./FChebRoots.hpp"
#include "./FChebMapping.hpp"
#include "./../Interpolation/FInterpTensor.hpp"
/**
......@@ -27,19 +27,6 @@
* Please read the license
*/
/**
* @class TensorTraits
*
* The class @p TensorTraits gives the number of interpolation nodes per
* cluster in 3D, depending on the interpolation order.
*
* @tparam ORDER interpolation order
*/
template <int ORDER> struct TensorTraits
{
enum {nnodes = ORDER*ORDER*ORDER};
};
/**
......@@ -51,123 +38,37 @@ template <int ORDER> struct TensorTraits
* @tparam ORDER interpolation order \f$\ell\f$
*/
template <int ORDER>
class FChebTensor : FNoCopyable
class FChebTensor : public FInterpTensor<ORDER,FChebRoots<ORDER>>
{
enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FChebRoots<ORDER> BasisType;
public:
/**
* Sets the ids of the coordinates of all \f$\ell^3\f$ interpolation
* nodes
*
* @param[out] NodeIds ids of coordinates of interpolation nodes
*/
static
void setNodeIds(unsigned int NodeIds[nnodes][3])
{
for (unsigned int n=0; n<nnodes; ++n) {
NodeIds[n][0] = n % ORDER;
NodeIds[n][1] = (n/ ORDER) % ORDER;
NodeIds[n][2] = n/(ORDER * ORDER);
}
}
/**
* Sets the roots of the Chebyshev quadrature weights defined as \f$w_i =
* \frac{\pi}{\ell}\sqrt{1-\bar x_i^2}\f$ with the Chebyshev roots \f$\bar
* x\f$.
*
* @param weights[out] the root of the weights \f$\sqrt{w_i}\f$
*/
static
void setRootOfWeights(FReal weights[nnodes])
{
// weights in 1d
FReal weights_1d[ORDER];
for (unsigned int o=0; o<ORDER; ++o)
weights_1d[o] = FMath::FPi/ORDER * FMath::Sqrt(FReal(1.)-FReal(BasisType::roots[o])*FReal(BasisType::roots[o]));
// weights in 3d (tensor structure)
unsigned int node_ids[nnodes][3];
setNodeIds(node_ids);
for (unsigned int n=0; n<nnodes; ++n) {
weights[n] = FMath::Sqrt(weights_1d[node_ids[n][0]]*weights_1d[node_ids[n][1]]*weights_1d[node_ids[n][2]]);
}
}
/**
* Sets the interpolation points in the cluster with @p center and @p width
*
* @param[in] center of cluster
* @param[in] width of cluster
* @param[out] rootPositions coordinates of interpolation points
*/
static
void setRoots(const FPoint& center, const FReal width, FPoint rootPositions[nnodes])
{
unsigned int node_ids[nnodes][3];
setNodeIds(node_ids);
const map_loc_glob map(center, width);
FPoint localPosition;
enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FChebRoots<ORDER> BasisType;
typedef FInterpTensor<ORDER,BasisType> ParentTensor;
public:
/**
* Sets the roots of the Chebyshev quadrature weights defined as \f$w_i =
* \frac{\pi}{\ell}\sqrt{1-\bar x_i^2}\f$ with the Chebyshev roots \f$\bar
* x\f$.
*
* @param weights[out] the root of the weights \f$\sqrt{w_i}\f$
*/
static
void setRootOfWeights(FReal weights[nnodes])
{
// weights in 1d
FReal weights_1d[ORDER];
for (unsigned int o=0; o<ORDER; ++o)
weights_1d[o] = FMath::FPi/ORDER * FMath::Sqrt(FReal(1.)-FReal(BasisType::roots[o])*FReal(BasisType::roots[o]));
// weights in 3d (tensor structure)
unsigned int node_ids[nnodes][3];
ParentTensor::setNodeIds(node_ids);
for (unsigned int n=0; n<nnodes; ++n) {
localPosition.setX(FReal(BasisType::roots[node_ids[n][0]]));
localPosition.setY(FReal(BasisType::roots[node_ids[n][1]]));
localPosition.setZ(FReal(BasisType::roots[node_ids[n][2]]));
map(localPosition, rootPositions[n]);
}
}
weights[n] = FMath::Sqrt(weights_1d[node_ids[n][0]]*weights_1d[node_ids[n][1]]*weights_1d[node_ids[n][2]]);
}
}
/**