Commit 49873018 authored by Pierre BLANCHARD's avatar Pierre BLANCHARD

Renamed FChebCmp to FCheb (original name), and FCheb to FChebDense.

parent e9d279f3
......@@ -17,8 +17,8 @@
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FCHEBCMPKERNEL_HPP
#define FCHEBCMPKERNEL_HPP
#ifndef FCHEBDENSEKERNEL_HPP
#define FCHEBDENSEKERNEL_HPP
#include "../../Utils/FGlobal.hpp"
......@@ -26,13 +26,13 @@
#include "./FAbstractChebKernel.hpp"
#include "./FChebCmpM2LHandler.hpp"
#include "./FChebDenseM2LHandler.hpp"
class FTreeCoordinate;
/**
* @author Matthias Messner(matthias.messner@inria.fr)
* @class FChebCmpKernel
* @class FChebDenseKernel
* @brief Chebyshev interpolation based FMM operators for general non oscillatory kernels..
*
* This class implements the Chebyshev interpolation based FMM operators. It
......@@ -52,11 +52,11 @@ class FTreeCoordinate;
* approximation of the M2L operators, then you can try to set EPSILON to 10 ^- (ORDER+{1,2,...}).
*/
template < class FReal, class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FChebCmpKernel
class FChebDenseKernel
: public FAbstractChebKernel< FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
{
// private types
typedef FChebCmpM2LHandler<FReal, ORDER,MatrixKernelClass> M2LHandlerClass;
typedef FChebDenseM2LHandler<FReal, ORDER,MatrixKernelClass> M2LHandlerClass;
// using from
typedef FAbstractChebKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
......@@ -81,7 +81,7 @@ public:
* but this will significantly slow down the computations.
*
*/
FChebCmpKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter, const MatrixKernelClass *const inMatrixKernel,
FChebDenseKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter, const MatrixKernelClass *const inMatrixKernel,
const FReal Epsilon)
: FAbstractChebKernel< FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>(inTreeHeight,inBoxWidth,inBoxCenter),
MatrixKernel(inMatrixKernel),
......@@ -100,8 +100,8 @@ public:
* Same as \see above constructor but the epsilon is automatically set to EPSILON=10^-ORDER
*/
FChebCmpKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter, const MatrixKernelClass *const inMatrixKernel)
: FChebCmpKernel(inTreeHeight, inBoxWidth,inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
FChebDenseKernel(const int inTreeHeight, const FReal inBoxWidth, const FPoint<FReal>& inBoxCenter, const MatrixKernelClass *const inMatrixKernel)
: FChebDenseKernel(inTreeHeight, inBoxWidth,inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
{
}
......@@ -111,15 +111,8 @@ public:
const ContainerClass* const SourceParticles)
{
const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
// 1) apply Sy
AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
LeafCell->getMultipole(0), SourceParticles);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 2) apply B
M2LHandler->applyB(LeafCell->getMultipole(idxRhs), LeafCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
}
}
......@@ -128,73 +121,34 @@ public:
const int /*TreeLevel*/)
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs),
ParentCell->getMultipole(idxRhs));
}
}
// 2) apply B
M2LHandler->applyB(ParentCell->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs) + AbstractBaseClass::nnodes);
}
}
// void M2L(CellClass* const FRestrict TargetCell,
// const CellClass* SourceCells[],
// const int NumSourceCells,
// const int TreeLevel) const
// {
// const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
// const FTreeCoordinate& cx = TargetCell->getCoordinate();
// for (int idx=0; idx<NumSourceCells; ++idx) {
// const FTreeCoordinate& cy = SourceCells[idx]->getCoordinate();
// const int transfer[3] = {cy.getX()-cx.getX(),
// cy.getY()-cx.getY(),
// cy.getZ()-cx.getZ()};
// M2LHandler->applyC(transfer, CellWidth,
// SourceCells[idx]->getMultipole() + AbstractBaseClass::nnodes,
// TargetCell->getLocal() + AbstractBaseClass::nnodes);
// }
// }
void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
const int neighborPositions[], const int inSize, const int TreeLevel) override {
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
FReal *const CompressedLocalExpansion = TargetCell->getLocal(idxRhs) + AbstractBaseClass::nnodes;
FReal *const localExpansion = TargetCell->getLocal(idxRhs);
const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
const int idx = neighborPositions[idxExistingNeigh];
M2LHandler->applyC(idx, CellWidth, SourceCells[idxExistingNeigh]->getMultipole(idxRhs) + AbstractBaseClass::nnodes,
CompressedLocalExpansion);
M2LHandler->applyC(idx, CellWidth, SourceCells[idxExistingNeigh]->getMultipole(idxRhs),
localExpansion);
}
}
}
// void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
// const int neighborPositions[], const int inSize, const int TreeLevel) override {
// const unsigned int rank = M2LHandler.getRank();
// FBlas::scal(343*rank, FReal(0.), MultipoleExpansion);
// const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
// for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
// const int idx = neighborPositions[idxExistingNeigh];
// FBlas::copy(rank, const_cast<FReal *const>(SourceCells[idxExistingNeigh]->getMultipole())+AbstractBaseClass::nnodes,
// MultipoleExpansion+idx*rank);
//
// M2LHandler->applyC(CellWidth, MultipoleExpansion, TargetCell->getLocal() + AbstractBaseClass::nnodes);
// }
void L2L(const CellClass* const FRestrict ParentCell,
CellClass* FRestrict *const FRestrict ChildCells,
const int /*TreeLevel*/)
{
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply U
M2LHandler->applyU(ParentCell->getLocal(idxRhs) + AbstractBaseClass::nnodes,
const_cast<CellClass*>(ParentCell)->getLocal(idxRhs));
// 2) apply Sx
for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
if (ChildCells[ChildIndex]){
AbstractBaseClass::Interpolator->applyL2L(ChildIndex, ParentCell->getLocal(idxRhs), ChildCells[ChildIndex]->getLocal(idxRhs));
......@@ -208,11 +162,6 @@ public:
{
const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply U
M2LHandler->applyU(LeafCell->getLocal(idxRhs) + AbstractBaseClass::nnodes, const_cast<CellClass*>(LeafCell)->getLocal(idxRhs));
}
//// 2.a) apply Sx
//AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
// AbstractBaseClass::BoxWidthLeaf,
......
......@@ -17,8 +17,8 @@
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
#ifndef FCHEBCMPM2LHANDLER_HPP
#define FCHEBCMPM2LHANDLER_HPP
#ifndef FCHEBDENSEM2LHANDLER_HPP
#define FCHEBDENSEM2LHANDLER_HPP
#include <numeric>
#include <stdexcept>
......@@ -32,32 +32,20 @@
#include "FChebTensor.hpp"
#include "Utils/FSvd.hpp"
template <class FReal, int ORDER>
unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
FReal* &U, FReal* &C, FReal* &B);
/**
* @author Matthias Messner (matthias.messner@inria.fr)
* @class FChebCmpM2LHandler
* @author Pierre Blanchard (pierre.blanchard@inria.fr)
* @class FChebDenseM2LHandler
* Please read the license
*
* This class precomputes and compresses the M2L operators
* This class precomputes the M2L operators
* \f$[K_1,\dots,K_{316}]\f$ for all (\f$7^3-3^3 = 316\f$ possible interacting
* cells in the far-field) interactions for the Chebyshev interpolation
* approach. The class uses the compression via a truncated SVD and represents
* the compressed M2L operator as \f$K_t \sim U C_t B^\top\f$ with
* \f$t=1,\dots,316\f$. The truncation rank is denoted by \f$r\f$ and is
* determined by the prescribed accuracy \f$\varepsilon\f$. Hence, the
* originally \f$K_t\f$ of size \f$\ell^3\times\ell^3\f$ times \f$316\f$ for
* all interactions is reduced to only one \f$U\f$ and one \f$B\f$, each of
* size \f$\ell^3\times r\f$, and \f$316\f$ \f$C_t\f$, each of size \f$r\times
* r\f$.
* approach.
*
* PB: FChebCmpM2LHandler does not seem to support non_homogeneous kernels!
* In fact nothing appears to handle this here (i.e. adapt scaling and storage
* PB: FChebDenseM2LHandler does not seem to support non_homogeneous kernels!
* In fact nothing appears to handle this here (i.e., adapt scaling and storage
* to MatrixKernelClass::Type). Given the relatively important cost of the
* Chebyshev variant, it is probably a choice not to have implemented this
* feature here but instead in the ChebyshevSym variant. But what if the
......@@ -69,7 +57,7 @@ unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
* @tparam ORDER interpolation order \f$\ell\f$
*/
template <class FReal, int ORDER, class MatrixKernelClass>
class FChebCmpM2LHandler : FNoCopyable
class FChebDenseM2LHandler : FNoCopyable
{
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes,
......@@ -77,102 +65,85 @@ class FChebCmpM2LHandler : FNoCopyable
const MatrixKernelClass *const MatrixKernel;
FReal *U, *C, *B;
const FReal epsilon; //<! accuracy which determines trucation of SVD
FReal *C;
unsigned int rank; //<! truncation rank, satisfies @p epsilon
static const std::string getFileName(FReal epsilon)
static const std::string getFileName()
{
const char precision_type = (typeid(FReal)==typeid(double) ? 'd' : 'f');
std::stringstream stream;
stream << "m2l_k"<< MatrixKernelClass::getID() << "_" << precision_type
<< "_o" << order << "_e" << epsilon << ".bin";
<< "_o" << order << ".bin";
return stream.str();
}
public:
FChebCmpM2LHandler(const MatrixKernelClass *const inMatrixKernel, const FReal _epsilon)
: MatrixKernel(inMatrixKernel), U(nullptr), C(nullptr), B(nullptr), epsilon(_epsilon), rank(0)
FChebDenseM2LHandler(const MatrixKernelClass *const inMatrixKernel, const FReal dummy)
: MatrixKernel(inMatrixKernel), C(nullptr), rank(0)
{}
~FChebCmpM2LHandler()
~FChebDenseM2LHandler()
{
if (U != nullptr) delete [] U;
if (B != nullptr) delete [] B;
if (C != nullptr) delete [] C;
}
/**
* Computes, compresses and sets the matrices \f$Y, C_t, B\f$
* Computes and sets the matrices \f$C_t\f$
*/
void ComputeAndCompressAndSet()
{
// measure time
FTic time; time.tic();
// check if aready set
if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);
if (C) throw std::runtime_error("Full M2L operator already set");
rank = Compute(MatrixKernel, C);
unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);
// write info
std::cout << "Compressed and set M2L operators (" << long(sizeM2L) << " B) in "
std::cout << "Compute and set full M2L operators (" << long(sizeM2L) << " B) in "
<< time.tacAndElapsed() << "sec." << std::endl;
}
/**
* Computes, compresses, writes to binary file, reads it and sets the matrices \f$Y, C_t, B\f$
*/
void ComputeAndCompressAndStoreInBinaryFileAndReadFromFileAndSet()
void ComputeAndStoreInBinaryFileAndReadFromFileAndSet()
{
FChebCmpM2LHandler<FReal, ORDER,MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(epsilon);
FChebDenseM2LHandler<FReal, ORDER,MatrixKernelClass>::ComputeAndStoreInBinaryFile();
this->ReadFromBinaryFileAndSet();
}
/**
* Computes and compressed all \f$K_t\f$.
*
* @param[in] epsilon accuracy
* @param[out] U matrix of size \f$\ell^3\times r\f$
* @param[in] MatrixKernel kernel function evaluator
* @param[out] C matrix of size \f$r\times 316 r\f$ storing \f$[C_1,\dots,C_{316}]\f$
* @param[out] B matrix of size \f$\ell^3\times r\f$
*/
static unsigned int ComputeAndCompress(const MatrixKernelClass *const MatrixKernel, const FReal epsilon, FReal* &U, FReal* &C, FReal* &B);
static unsigned int Compute(const MatrixKernelClass *const MatrixKernel, FReal* &C);
/**
* Computes, compresses and stores the matrices \f$Y, C_t, B\f$ in a binary
* Computes and stores the matrices \f$C_t\f$ in a binary
* file
*/
static void ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon);
static void ComputeAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel);
/**
* Reads the matrices \f$Y, C_t, B\f$ from the respective binary file
* Reads the matrices \f$C_t\f$ from the respective binary file
*/
void ReadFromBinaryFileAndSet();
/**
* @return rank of the SVD compressed M2L operators
* @return rank of the M2L operators
*/
unsigned int getRank() const {return rank;}
/**
* Expands potentials \f$x+=UX\f$ of a target cell. This operation can be
* seen as part of the L2L operation.
*
* @param[in] X compressed local expansion of size \f$r\f$
* @param[out] x local expansion of size \f$\ell^3\f$
*/
void applyU(const FReal *const X, FReal *const x) const
{
FBlas::gemva(nnodes, rank, 1., U, const_cast<FReal*>(X), x);
}
/**
* Compressed M2L operation \f$X+=C_tY\f$, where \f$Y\f$ is the compressed
* Full M2L operation \f$X+=C_tY\f$, where \f$Y\f$ is the compressed
* multipole expansion and \f$X\f$ is the compressed local expansion, both
* of size \f$r\f$. The index \f$t\f$ denotes the transfer vector of the
* target cell to the source cell.
......@@ -183,39 +154,26 @@ public:
* @param[in] CellWidth needed for the scaling of the compressed M2L operators which are based on a homogeneous matrix kernel computed for the reference cell width \f$w=2\f$, ie in \f$[-1,1]^3\f$.
*/
void applyC(const int transfer[3], FReal CellWidth,
const FReal *const Y, FReal *const X) const
const FReal *const Y, FReal *const X) const
{
const unsigned int idx
= (transfer[0]+3)*7*7 + (transfer[1]+3)*7 + (transfer[2]+3);
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
const unsigned int idx
= (transfer[0]+3)*7*7 + (transfer[1]+3)*7 + (transfer[2]+3);
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
}
void applyC(const unsigned int idx, FReal CellWidth,
const FReal *const Y, FReal *const X) const
const FReal *const Y, FReal *const X) const
{
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank, scale, C + idx*rank*rank, const_cast<FReal*>(Y), X);
}
void applyC(FReal CellWidth,
const FReal *const Y, FReal *const X) const
const FReal *const Y, FReal *const X) const
{
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
const FReal scale(MatrixKernel->getScaleFactor(CellWidth));
FBlas::gemva(rank, rank * 343, scale, C, const_cast<FReal*>(Y), X);
}
/**
* Compresses densities \f$Y=B^\top y\f$ of a source cell. This operation
* can be seen as part of the M2M operation.
*
* @param[in] y multipole expansion of size \f$\ell^3\f$
* @param[out] Y compressed multipole expansion of size \f$r\f$
*/
void applyB(FReal *const y, FReal *const Y) const
{
FBlas::gemtv(nnodes, rank, 1., B, y, Y);
}
};
......@@ -234,14 +192,10 @@ public:
template <class FReal, int ORDER, class MatrixKernelClass>
unsigned int
FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompress(const MatrixKernelClass *const MatrixKernel,
const FReal epsilon,
FReal* &U,
FReal* &C,
FReal* &B)
FChebDenseM2LHandler<FReal, ORDER, MatrixKernelClass>::Compute(const MatrixKernelClass *const MatrixKernel, FReal* &C)
{
// allocate memory and store compressed M2L operators
if (U||C||B) throw std::runtime_error("Compressed M2L operators are already set");
if (C) throw std::runtime_error("Full M2L operators are already set");
// interpolation points of source (Y) and target (X) cell
FPoint<FReal> X[nnodes], Y[nnodes];
......@@ -249,8 +203,7 @@ FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompress(const Ma
FChebTensor<FReal, order>::setRoots(FPoint<FReal>(0.,0.,0.), FReal(2.), X);
// allocate memory and compute 316 m2l operators
FReal *_U, *_C, *_B;
_U = _B = nullptr;
FReal *_C;
_C = new FReal [nnodes*nnodes * ninteractions];
unsigned int counter = 0;
for (int i=-3; i<=3; ++i) {
......@@ -272,32 +225,11 @@ FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompress(const Ma
}
}
if (counter != ninteractions)
throw std::runtime_error("Number of interactions must correspond to 316");
//////////////////////////////////////////////////////////
FReal weights[nnodes];
FChebTensor<FReal, order>::setRootOfWeights(weights);
for (unsigned int i=0; i<316; ++i)
for (unsigned int n=0; n<nnodes; ++n) {
FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n, nnodes); // scale rows
FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n * nnodes); // scale cols
}
//////////////////////////////////////////////////////////
throw std::runtime_error("Number of interactions must correspond to 316");
// svd compression of M2L
const unsigned int rank = Compress<FReal, ORDER>(epsilon, ninteractions, _U, _C, _B);
if (!(rank>0)) throw std::runtime_error("Low rank must be larger then 0!");
// store U
U = new FReal [nnodes * rank];
FBlas::copy(rank*nnodes, _U, U);
delete [] _U;
// store B
B = new FReal [nnodes * rank];
FBlas::copy(rank*nnodes, _B, B);
delete [] _B;
const unsigned int rank = nnodes;
// store C
counter = 0;
C = new FReal [343 * rank*rank];
......@@ -312,16 +244,7 @@ FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompress(const Ma
}
if (counter != ninteractions)
throw std::runtime_error("Number of interactions must correspond to 316");
delete [] _C;
//////////////////////////////////////////////////////////
for (unsigned int n=0; n<nnodes; ++n) {
FBlas::scal(rank, FReal(1.) / weights[n], U+n, nnodes); // scale rows
FBlas::scal(rank, FReal(1.) / weights[n], B+n, nnodes); // scale rows
}
//////////////////////////////////////////////////////////
delete [] _C;
// return low rank
return rank;
......@@ -334,18 +257,17 @@ FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompress(const Ma
template <class FReal, int ORDER, class MatrixKernelClass>
void
FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel, const FReal epsilon)
FChebDenseM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndStoreInBinaryFile(const MatrixKernelClass *const MatrixKernel)
{
// measure time
FTic time; time.tic();
// start computing process
FReal *U, *C, *B;
U = C = B = nullptr;
const unsigned int rank = ComputeAndCompress(MatrixKernel, epsilon, U, C, B);
FReal *C;
C = nullptr;
const unsigned int rank = nnodes;
// store into binary file
const std::string filename(getFileName(epsilon));
std::ofstream stream(filename.c_str(),
std::ios::out | std::ios::binary | std::ios::trunc);
const std::string filename(getFileName());
std::ofstream stream(filename.c_str(), std::ios::out | std::ios::binary | std::ios::trunc);
if (stream.good()) {
stream.seekp(0);
// 1) write number of interpolation points (int)
......@@ -354,40 +276,34 @@ FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreI
// 2) write low rank (int)
int _rank = rank;
stream.write(reinterpret_cast<char*>(&_rank), sizeof(int));
// 3) write U (rank*nnodes * FReal)
stream.write(reinterpret_cast<char*>(U), sizeof(FReal)*rank*nnodes);
// 4) write B (rank*nnodes * FReal)
stream.write(reinterpret_cast<char*>(B), sizeof(FReal)*rank*nnodes);
// 5) write 343 C (343 * rank*rank * FReal)
stream.write(reinterpret_cast<char*>(C), sizeof(FReal)*rank*rank*343);
} else throw std::runtime_error("File could not be opened to write");
stream.close();
// free memory
if (U != nullptr) delete [] U;
if (B != nullptr) delete [] B;
if (C != nullptr) delete [] C;
// write info
std::cout << "Compressed M2L operators ("<< rank << ") stored in binary file " << filename
std::cout << "Full M2L operators ("<< rank << ") stored in binary file " << filename
<< " in " << time.tacAndElapsed() << "sec." << std::endl;
}
template <class FReal, int ORDER, class MatrixKernelClass>
void
FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ReadFromBinaryFileAndSet()
FChebDenseM2LHandler<FReal, ORDER, MatrixKernelClass>::ReadFromBinaryFileAndSet()
{
// measure time
FTic time; time.tic();
// start reading process
if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
const std::string filename(getFileName(epsilon));
if (C) throw std::runtime_error("Full M2L operator already set");
const std::string filename(getFileName());
std::ifstream stream(filename.c_str(),
std::ios::in | std::ios::binary | std::ios::ate);
const std::ifstream::pos_type size = stream.tellg();
if (size<=0) {
std::cout << "Info: The requested binary file " << filename
<< " does not yet exist. Compute it now ... " << std::endl;
this->ComputeAndCompressAndStoreInBinaryFileAndReadFromFileAndSet();
this->ComputeAndStoreInBinaryFileAndReadFromFileAndSet();
return;
}
if (stream.good()) {
......@@ -398,147 +314,19 @@ FChebCmpM2LHandler<FReal, ORDER, MatrixKernelClass>::ReadFromBinaryFileAndSet()
if (npts!=nnodes) throw std::runtime_error("nnodes and npts do not correspond");
// 2) read low rank (int)
stream.read(reinterpret_cast<char*>(&rank), sizeof(int));
// 3) write U (rank*nnodes * FReal)
U = new FReal [rank*nnodes];
stream.read(reinterpret_cast<char*>(U), sizeof(FReal)*rank*nnodes);
// 4) write B (rank*nnodes * FReal)
B = new FReal [rank*nnodes];
stream.read(reinterpret_cast<char*>(B), sizeof(FReal)*rank*nnodes);
// 5) write 343 C (343 * rank*rank * FReal)
C = new FReal [343 * rank*rank];
stream.read(reinterpret_cast<char*>(C), sizeof(FReal)*rank*rank*343);
} else throw std::runtime_error("File could not be opened to read");
stream.close();
// write info
std::cout << "Compressed M2L operators (" << rank << ") read from binary file "
std::cout << "Full M2L operators (" << rank << ") read from binary file "
<< filename << " in " << time.tacAndElapsed() << "sec." << std::endl;
}
/*
unsigned int ReadRankFromBinaryFile(const std::string& filename)
{
// start reading process
std::ifstream stream(filename.c_str(), std::ios::in | std::ios::binary | std::ios::ate);
const std::ifstream::pos_type size = stream.tellg();
if (size<=0) throw std::runtime_error("The requested binary file does not exist.");
unsigned int rank = -1;
if (stream.good()) {
stream.seekg(0);
// 1) read number of interpolation points (int)
int npts;
stream.read(reinterpret_cast<char*>(&npts), sizeof(int));
// 2) read low rank (int)
stream.read(reinterpret_cast<char*>(&rank), sizeof(int));
return rank;
} else throw std::runtime_error("File could not be opened to read");
stream.close();
return rank;
}
*/
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
/**
* Compresses \f$[K_1,\dots,K_{316}]\f$ in \f$C\f$. Attention: the matrices
* \f$U,B\f$ are not initialized, no memory is allocated as input, as output
* they store the respective matrices. The matrix \f$C\f$ stores
* \f$[K_1,\dots,K_{316}]\f$ as input and \f$[C_1,\dots,C_{316}]\f$ as output.
*
* @param[in] epsilon accuracy
* @param[out] U matrix of size \f$\ell^3\times r\f$
* @param[in] C matrix of size \f$\ell^3\times 316 \ell^e\f$ storing \f$[K_1,\dots,K_{316}]\f$
* @param[out] C matrix of size \f$r\times 316 r\f$ storing \f$[C_1,\dots,C_{316}]\f$
* @param[out] B matrix of size \f$\ell^3\times r\f$
*/
template <class FReal, int ORDER>
unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
FReal* &U, FReal* &C, FReal* &B)
{
// compile time constants
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes};
// init SVD
const unsigned int LWORK = 2 * (3*nnodes + ninteractions*nnodes);
FReal *const WORK = new FReal [LWORK];
// K_col ///////////////////////////////////////////////////////////
FReal *const K_col = new FReal [ninteractions * nnodes*nnodes];
for (unsigned int i=0; i<ninteractions; ++i)