Commit e62560e2 authored by BRAMAS Berenger's avatar BRAMAS Berenger
parents fca208c7 086c6a53
......@@ -46,7 +46,8 @@ SET(my_include_dirs "."
"Arranger" "Core" "Utils"
"Kernels/Chebyshev" "Components" "Extensions"
"Containers" "Files" "Kernels/Spherical"
"Kernels/P2P" "Kernels/Taylor" "Kernels/ROtation")
"Kernels/P2P" "Kernels/Taylor" "Kernels/ROtation"
"Kernels/Uniform" "Kernels/Interpolation")
FOREACH(my_dir ${my_include_dirs})
file(GLOB
......
......@@ -31,18 +31,18 @@
* This class defines a cell used in the Chebyshev based FMM.
* @param NVALS is the number of right hand side.
*/
template <int ORDER, int NVALS = 1>
template <int ORDER, int NRHS = 1, int NLHS = 1, int NVALS = 1>
class FChebCell : public FBasicCell
{
static const int VectorSize = TensorTraits<ORDER>::nnodes * 2;
FReal multipole_exp[NVALS * VectorSize]; //< Multipole expansion
FReal local_exp[NVALS * VectorSize]; //< Local expansion
FReal multipole_exp[NRHS * NVALS * VectorSize]; //< Multipole expansion
FReal local_exp[NLHS * NVALS * VectorSize]; //< Local expansion
public:
FChebCell(){
memset(multipole_exp, 0, sizeof(FReal) * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NVALS * VectorSize);
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
}
~FChebCell() {}
......@@ -72,8 +72,8 @@ public:
/** Make it like the begining */
void resetToInitialState(){
memset(multipole_exp, 0, sizeof(FReal) * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NVALS * VectorSize);
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
}
///////////////////////////////////////////////////////
......@@ -81,20 +81,20 @@ public:
///////////////////////////////////////////////////////
template <class BufferWriterClass>
void serializeUp(BufferWriterClass& buffer) const{
buffer.write(multipole_exp, VectorSize*NVALS);
buffer.write(multipole_exp, VectorSize*NVALS*NRHS);
}
template <class BufferReaderClass>
void deserializeUp(BufferReaderClass& buffer){
buffer.fillArray(multipole_exp, VectorSize*NVALS);
buffer.fillArray(multipole_exp, VectorSize*NVALS*NRHS);
}
template <class BufferWriterClass>
void serializeDown(BufferWriterClass& buffer) const{
buffer.write(local_exp, VectorSize*NVALS);
buffer.write(local_exp, VectorSize*NVALS*NLHS);
}
template <class BufferReaderClass>
void deserializeDown(BufferReaderClass& buffer){
buffer.fillArray(local_exp, VectorSize*NVALS);
buffer.fillArray(local_exp, VectorSize*NVALS*NLHS);
}
///////////////////////////////////////////////////////
......@@ -103,37 +103,37 @@ public:
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FBasicCell::save(buffer);
buffer.write(multipole_exp, VectorSize);
buffer.write(local_exp, VectorSize);
buffer.write(multipole_exp, VectorSize*NVALS*NRHS);
buffer.write(local_exp, VectorSize*NVALS*NLHS);
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FBasicCell::restore(buffer);
buffer.fillArray(multipole_exp, VectorSize);
buffer.fillArray(local_exp, VectorSize);
buffer.fillArray(multipole_exp, VectorSize*NVALS*NRHS);
buffer.fillArray(local_exp, VectorSize*NVALS*NLHS);
}
static int GetSize(){
return sizeof(FReal)*VectorSize*2*NVALS;
return sizeof(FReal) * VectorSize*(NRHS+NLHS)*NVALS;
}
};
template <int ORDER, int NVALS = 1>
class FTypedChebCell : public FChebCell<ORDER,NVALS>, public FExtendCellType {
template <int ORDER, int NRHS = 1, int NLHS = 1, int NVALS = 1>
class FTypedChebCell : public FChebCell<ORDER,NRHS,NLHS,NVALS>, public FExtendCellType {
public:
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FChebCell<ORDER,NVALS>::save(buffer);
FChebCell<ORDER,NRHS,NLHS,NVALS>::save(buffer);
FExtendCellType::save(buffer);
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FChebCell<ORDER,NVALS>::restore(buffer);
FChebCell<ORDER,NRHS,NLHS,NVALS>::restore(buffer);
FExtendCellType::restore(buffer);
}
void resetToInitialState(){
FChebCell<ORDER,NVALS>::resetToInitialState();
FChebCell<ORDER,NRHS,NLHS,NVALS>::resetToInitialState();
FExtendCellType::resetToInitialState();
}
};
......
......@@ -25,7 +25,7 @@
#include "../../Components/FAbstractKernels.hpp"
#include "./FChebInterpolator.hpp"
#include "./FChebSymmetries.hpp"
#include "../Interpolation/FInterpSymmetries.hpp"
class FTreeCoordinate;
......@@ -322,7 +322,7 @@ struct FChebFlopsSymKernel<CellClass, ContainerClass, MatrixKernelClass, ORDER>
}
// set permutation vector and indices
const FChebSymmetries<ORDER> Symmetries;
const FInterpSymmetries<ORDER> Symmetries;
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
for (int k=-3; k<=3; ++k)
......
......@@ -21,7 +21,7 @@
#include "../../Utils/FBlas.hpp"
#include "./FChebTensor.hpp"
#include "./FChebSymmetries.hpp"
#include "../Interpolation/FInterpSymmetries.hpp"
#include "./FChebM2LHandler.hpp"
/**
......@@ -539,7 +539,7 @@ public:
}
// set permutation vector and indices
const FChebSymmetries<ORDER> Symmetries;
const FInterpSymmetries<ORDER> Symmetries;
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
for (int k=-3; k<=3; ++k) {
......@@ -620,7 +620,7 @@ public:
// set permutation vector and indices
const FChebSymmetries<ORDER> Symmetries;
const FInterpSymmetries<ORDER> Symmetries;
for (int i=-3; i<=3; ++i)
for (int j=-3; j<=3; ++j)
for (int k=-3; k<=3; ++k) {
......
This diff is collapsed.
This diff is collapsed.
......@@ -48,6 +48,42 @@ struct DirectInteractionComputer<LENNARD_JONES_POTENTIAL, 1>
}
};
/*! Specialization for ID_OVER_R potential */
template <>
struct DirectInteractionComputer<ID_OVER_R, 1>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutualIOR(TargetParticles,NeighborSourceParticles,14);
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemoteIOR(inTargets,inNeighbors,inSize);
}
};
/*! Specialization for GradGradR potential */
template <>
struct DirectInteractionComputer<R_IJ, 1>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutualRIJ(TargetParticles,NeighborSourceParticles,14);
}
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemoteRIJ(inTargets,inNeighbors,inSize);
}
};
///////////////////////////////////////////////////////
// In case of multi right hand side
///////////////////////////////////////////////////////
......@@ -98,4 +134,48 @@ struct DirectInteractionComputer<LENNARD_JONES_POTENTIAL, NVALS>
}
};
/*! Specialization for ID_OVER_R potential */
template <int NVALS>
struct DirectInteractionComputer<ID_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::FullMutualIOR(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::FullRemoteIOR(inTargets,inNeighbors,inSize);
}
}
};
/*! Specialization for GradGradR potential */
template <int NVALS>
struct DirectInteractionComputer<R_IJ, NVALS>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FP2P::FullMutualRIJ(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::FullRemoteRIJ(inTargets,inNeighbors,inSize);
}
}
};
#endif // FINTERPP2PKERNELS_HPP
......@@ -13,8 +13,8 @@
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FCHEBSYMMETRIES_HPP
#define FCHEBSYMMETRIES_HPP
#ifndef FINTERPSYMMETRIES_HPP
#define FINTERPSYMMETRIES_HPP
#include <climits>
......@@ -25,12 +25,12 @@
*/
/**
* @class FChebSymmetries
* @class FInterpSymmetries
*
* The class @p FChebSymmetries exploits all symmetries
* The class @p FInterpSymmetries exploits all symmetries
*/
template <int ORDER>
class FChebSymmetries
class FInterpSymmetries
{
enum {nnodes = ORDER*ORDER*ORDER};
......@@ -49,11 +49,11 @@ class FChebSymmetries
return (sk | sj | si);
}
public:
/** Constructor */
FChebSymmetries()
FInterpSymmetries()
{
// permutations for 8 quadrants
unsigned int quads[8][nnodes];
......
This diff is collapsed.
......@@ -45,7 +45,7 @@ class FAbstractUnifKernel : public FAbstractKernels< CellClass, ContainerClass>
{
protected:
enum {nnodes = TensorTraits<ORDER>::nnodes};
typedef FUnifInterpolator<ORDER> InterpolatorClass;
typedef FUnifInterpolator<ORDER,MatrixKernelClass> InterpolatorClass;
/// Needed for P2M, M2M, L2L and L2P operators
const FSmartPointer<InterpolatorClass,FSmartPointerMemory> Interpolator;
......
......@@ -33,36 +33,31 @@
*
* This class defines a cell used in the Lagrange based FMM.
*
* PB: !!! exactly the same as FChebCell except the TensorTraits used is
* the one from FUnifTensor. This is a quick work around to avoid multiple
* definition of TensorTraits until it is factorized for all interpolations.
*
* PB: This class also contains the storage and accessors for the tranformed
* expansion (in Fourier space, i.e. complex valued).
* PB: This class also contains the storage and accessors for the transformed
* expansion (in Fourier space, i.e. complex valued).
*
* @param NVALS is the number of right hand side.
*/
template <int ORDER, int NVALS = 1>
template <int ORDER, int NRHS = 1, int NLHS = 1, int NVALS = 1>
class FUnifCell : public FBasicCell
{
static const int VectorSize = TensorTraits<ORDER>::nnodes;
static const int TransformedVectorSize = (2*ORDER-1)*(2*ORDER-1)*(2*ORDER-1);
FReal multipole_exp[NVALS * VectorSize]; //< Multipole expansion
FReal local_exp[NVALS * VectorSize]; //< Local expansion
FReal multipole_exp[NRHS * NVALS * VectorSize]; //< Multipole expansion
FReal local_exp[NLHS * NVALS * VectorSize]; //< Local expansion
// PB: Store multipole and local expansion in Fourier space
FComplexe transformed_multipole_exp[NVALS * TransformedVectorSize];
FComplexe transformed_local_exp[NVALS * TransformedVectorSize];
FComplexe transformed_multipole_exp[NRHS * NVALS * TransformedVectorSize];
FComplexe transformed_local_exp[NLHS * NVALS * TransformedVectorSize];
public:
FUnifCell(){
memset(multipole_exp, 0, sizeof(FReal) * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NVALS * VectorSize);
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
memset(transformed_multipole_exp, 0,
sizeof(FComplexe) * NVALS * TransformedVectorSize);
sizeof(FComplexe) * NRHS * NVALS * TransformedVectorSize);
memset(transformed_local_exp, 0,
sizeof(FComplexe) * NVALS * TransformedVectorSize);
sizeof(FComplexe) * NLHS * NVALS * TransformedVectorSize);
}
~FUnifCell() {}
......@@ -90,19 +85,9 @@ public:
return VectorSize;
}
/** Make it like the begining */
void resetToInitialState(){
memset(multipole_exp, 0, sizeof(FReal) * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NVALS * VectorSize);
memset(transformed_multipole_exp, 0,
sizeof(FComplexe) * NVALS * TransformedVectorSize);
memset(transformed_local_exp, 0,
sizeof(FComplexe) * NVALS * TransformedVectorSize);
}
/** Get Transformed Multipole */
const FComplexe* getTransformedMultipole(const int inRhs) const
{ return this->transformed_multipole_exp + inRhs*TransformedVectorSize;
const FComplexe* getTransformedMultipole(const int inRhs) const{
return this->transformed_multipole_exp + inRhs*TransformedVectorSize;
}
/** Get Transformed Local */
const FComplexe* getTransformedLocal(const int inRhs) const{
......@@ -118,31 +103,46 @@ public:
return this->transformed_local_exp + inRhs*TransformedVectorSize;
}
/** To get the leading dim of a vec */
int getTransformedVectorSize() const{
return TransformedVectorSize;
}
/** Make it like the begining */
void resetToInitialState(){
memset(multipole_exp, 0, sizeof(FReal) * NRHS * NVALS * VectorSize);
memset(local_exp, 0, sizeof(FReal) * NLHS * NVALS * VectorSize);
memset(transformed_multipole_exp, 0,
sizeof(FComplexe) * NRHS * NVALS * TransformedVectorSize);
memset(transformed_local_exp, 0,
sizeof(FComplexe) * NLHS * NVALS * TransformedVectorSize);
}
///////////////////////////////////////////////////////
// to extend FAbstractSendable
///////////////////////////////////////////////////////
template <class BufferWriterClass>
void serializeUp(BufferWriterClass& buffer) const{
buffer.write(multipole_exp, VectorSize * NVALS);
buffer.write(transformed_multipole_exp, VectorSize * NVALS);
buffer.write(multipole_exp, VectorSize*NVALS*NRHS);
buffer.write(transformed_multipole_exp, TransformedVectorSize*NVALS*NRHS);
}
template <class BufferReaderClass>
void deserializeUp(BufferReaderClass& buffer){
buffer.fillArray(multipole_exp, VectorSize*NVALS);
buffer.fillArray(transformed_multipole_exp, VectorSize*NVALS);
buffer.fillArray(multipole_exp, VectorSize*NVALS*NRHS);
buffer.fillArray(transformed_multipole_exp, TransformedVectorSize*NVALS*NRHS);
}
template <class BufferWriterClass>
void serializeDown(BufferWriterClass& buffer) const{
buffer.write(local_exp, VectorSize*NVALS);
buffer.write(transformed_local_exp, VectorSize * NVALS);
buffer.write(local_exp, VectorSize*NVALS*NLHS);
buffer.write(transformed_local_exp, TransformedVectorSize*NVALS*NLHS);
}
template <class BufferReaderClass>
void deserializeDown(BufferReaderClass& buffer){
buffer.fillArray(local_exp, VectorSize*NVALS);
buffer.fillArray(transformed_local_exp, VectorSize*NVALS);
buffer.fillArray(local_exp, VectorSize*NVALS*NLHS);
buffer.fillArray(transformed_local_exp, TransformedVectorSize*NVALS*NLHS);
}
///////////////////////////////////////////////////////
......@@ -151,43 +151,43 @@ public:
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FBasicCell::save(buffer);
buffer.write(multipole_exp, VectorSize*NVALS);
buffer.write(transformed_multipole_exp, VectorSize*NVALS);
buffer.write(local_exp, VectorSize*NVALS);
buffer.write(transformed_local_exp, VectorSize*NVALS);
buffer.write(multipole_exp, VectorSize*NVALS*NRHS);
buffer.write(transformed_multipole_exp, TransformedVectorSize*NVALS*NRHS);
buffer.write(local_exp, VectorSize*NVALS*NLHS);
buffer.write(transformed_local_exp, TransformedVectorSize*NVALS*NLHS);
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FBasicCell::restore(buffer);
buffer.fillArray(multipole_exp, VectorSize*NVALS);
buffer.fillArray(transformed_multipole_exp, VectorSize*NVALS);
buffer.fillArray(local_exp, VectorSize*NVALS);
buffer.fillArray(transformed_local_exp, VectorSize*NVALS);
buffer.fillArray(multipole_exp, VectorSize*NVALS*NRHS);
buffer.fillArray(transformed_multipole_exp, TransformedVectorSize*NVALS*NRHS);
buffer.fillArray(local_exp, VectorSize*NVALS*NLHS);
buffer.fillArray(transformed_local_exp, TransformedVectorSize*NVALS*NLHS);
}
static int GetSize(){
return 2 * (int) sizeof(FReal)*VectorSize*NVALS + 2*NVALS*TransformedVectorSize*(int) sizeof(FComplexe);
return (NRHS+NLHS)*NVALS*VectorSize * (int) sizeof(FReal) + (NRHS+NLHS)*NVALS*TransformedVectorSize * (int) sizeof(FComplexe);
}
};
template <int ORDER, int NVALS = 1>
class FTypedUnifCell : public FUnifCell<ORDER,NVALS>, public FExtendCellType {
template <int ORDER, int NRHS = 1, int NLHS = 1, int NVALS = 1>
class FTypedUnifCell : public FUnifCell<ORDER,NRHS,NLHS,NVALS>, public FExtendCellType {
public:
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FUnifCell<ORDER,NVALS>::save(buffer);
FUnifCell<ORDER,NRHS,NLHS,NVALS>::save(buffer);
FExtendCellType::save(buffer);
}
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FUnifCell<ORDER,NVALS>::restore(buffer);
FUnifCell<ORDER,NRHS,NLHS,NVALS>::restore(buffer);
FExtendCellType::restore(buffer);
}
void resetToInitialState(){
FUnifCell<ORDER,NVALS>::resetToInitialState();
FUnifCell<ORDER,NRHS,NLHS,NVALS>::resetToInitialState();
FExtendCellType::resetToInitialState();
}
};
......
This diff is collapsed.
......@@ -96,7 +96,7 @@ public:
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
FBlas::scal(AbstractBaseClass::nnodes*2, FReal(0.), ParentCell->getMultipole(idxRhs));
FBlas::scal(AbstractBaseClass::nnodes, FReal(0.), ParentCell->getMultipole(idxRhs));
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
......@@ -170,7 +170,6 @@ public:
const int /*TreeLevel*/)
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) Apply Inverse Discete Fourier Transform
M2LHandler->unapplyZeroPaddingAndDFT(ParentCell->getTransformedLocal(idxRhs),
const_cast<CellClass*>(ParentCell)->getLocal(idxRhs));
......
......@@ -54,20 +54,22 @@ class FUnifM2LHandler : FNoCopyable
{
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes,
ninteractions = 316}; // 7^3 - 3^3 (max num cells in far-field)
ninteractions = 316, // 7^3 - 3^3 (max num cells in far-field)
rc = (2*ORDER-1)*(2*ORDER-1)*(2*ORDER-1)};
const MatrixKernelClass MatrixKernel;
FComplexe *FC;
unsigned int rc;
// for real valued kernel only n/2+1 complex values are stored
// after performing the DFT (the rest is deduced by conjugation)
unsigned int opt_rc;
typedef FUnifTensor<ORDER> TensorType;
unsigned int node_diff[nnodes*nnodes]; // PB: used in applyC to get id=i-j
unsigned int node_diff[nnodes*nnodes];
// FDft Dft; // Direct Discrete Fourier Transformator
FFft Dft; // Fast Discrete Fourier Transformator
FFft<1> Dft; // Fast Discrete Fourier Transformator
static const std::string getFileName()
{
......@@ -81,17 +83,12 @@ class FUnifM2LHandler : FNoCopyable
public:
FUnifM2LHandler()
: MatrixKernel(), FC(NULL),
rc((2*ORDER-1)*(2*ORDER-1)*(2*ORDER-1)),
: MatrixKernel(), FC(NULL),
opt_rc(rc/2+1),
Dft(rc) // initialize Discrete Fourier Transformator
{
// initialize root node ids
TensorType::setNodeIdsDiff(node_diff);
// for real valued kernel only n/2+1 complex values are stored
// after performing the DFT (the rest is deduced by conjugation)
}
~FUnifM2LHandler()
......@@ -222,18 +219,12 @@ public:
FReal Py[rc];
FBlas::setzero(rc,Py);
FComplexe tmpFY[rc]; // not mandatory?
// Apply Zero Padding
for (unsigned int i=0; i<nnodes; ++i)
Py[node_diff[i*nnodes]]=y[i];
// Apply forward Discrete Fourier Transform
Dft.applyDFT(Py,tmpFY);
// not mandatory?
for (unsigned int j=0; j<rc; ++j) // could be opt_rc
FY[j]+=tmpFY[j];
Dft.applyDFT(Py,FY);
}
......@@ -277,22 +268,19 @@ FUnifM2LHandler<ORDER, MatrixKernelClass>::Compute(FComplexe* &FC)
_C = new FReal [rc];
_FC = new FComplexe [rc * ninteractions];
// initialize root node ids pairs
unsigned int node_ids_pairs[rc][2];
TensorType::setNodeIdsPairs(node_ids_pairs);
// init Discrete Fourier Transformator
// FDft Dft(rc);
FFft Dft(rc);
FFft<1> Dft(rc);
// get first column of K via permutation
unsigned int perm[rc];
for(unsigned int p=0; p<rc; ++p){
// permutation to order WHILE computing the entire 1st row
if(p<rc-1) perm[p]=p+1;
else perm[p]=p+1-rc;
// std::cout << "perm["<< p << "]="<< perm[p] << std::endl;
}
TensorType::setStoragePermutation(perm);
unsigned int counter = 0;
unsigned int li,lj, mi,mj, ni,nj;
unsigned int idi, idj, ido;
for (int i=-3; i<=3; ++i) {
for (int j=-3; j<=3; ++j) {
for (int k=-3; k<=3; ++k) {
......@@ -301,23 +289,11 @@ FUnifM2LHandler<ORDER, MatrixKernelClass>::Compute(FComplexe* &FC)
const FPoint cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
FUnifTensor<order>::setRoots(cy, FReal(2.), Y);
// evaluate m2l operator
ido=0;
unsigned int ido=0;
for(unsigned int l=0; l<2*order-1; ++l)
for(unsigned int m=0; m<2*order-1; ++m)
for(unsigned int n=0; n<2*order-1; ++n){
// l=0:(2*order-1) => li-lj=-(order-1):(order-1)
// Convention:
// lj=order-1 & li=0:order-1 => li-lj=1-order:0
// lj=1 & li=0:order-1 => li-lj=1:order-1
if(l<order-1) lj=order-1; else lj=0;
if(m<order-1) mj=order-1; else mj=0;
if(n<order-1) nj=order-1; else nj=0;
li=(l-(order-1))+lj; mi=(m-(order-1))+mj; ni=(n-(order-1))+nj;
// deduce corresponding index of K[nnodes x nnodes]
idi=li*order*order + mi*order + ni;
idj=lj*order*order + mj*order + nj;
// store value at current position in C
// use permutation if DFT is used because
// the storage of the first column is required
......@@ -325,7 +301,8 @@ FUnifM2LHandler<ORDER, MatrixKernelClass>::Compute(FComplexe* &FC)
// i.e. C[rc-1] C[0] C[1] .. C[rc-2] < RIGHT!
// _C[counter*rc + ido]
_C[perm[ido]]
= MatrixKernel.evaluate(X[idi], Y[idj]);
= MatrixKernel.evaluate(X[node_ids_pairs[ido][0]],
Y[node_ids_pairs[ido][1]]);
ido++;
}
......
......@@ -63,14 +63,14 @@ struct FUnifRoots : FNoCopyable
x = (x < FReal(-1.) ? FReal(-1.) : x);
}