FChebTensorialM2LHandler.hpp 14.1 KB
Newer Older
1
// ===================================================================================
2 3 4 5
// Copyright ScalFmm 2016 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.
6
//
7
// This software is governed by the CeCILL-C and LGPL licenses and
8
// abiding by the rules of distribution of free software.
9 10 11
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
12
//
13 14 15
// 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
16 17 18
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
// ===================================================================================
#ifndef FCHEBTENSORIALM2LHANDLER_HPP
#define FCHEBTENSORIALM2LHANDLER_HPP

#include <numeric>
#include <stdexcept>
#include <string>
#include <sstream>
#include <fstream>
#include <typeinfo>

#include "../../Utils/FBlas.hpp"
#include "../../Utils/FTic.hpp"

#include "./FChebTensor.hpp"

/**
 * Computes and compresses all \f$K_t\f$.
 *
 * @param[in] epsilon accuracy
39
 * @param[out] C matrix of size \f$n\times 316 n\f$ storing \f$[C_1,\dots,C_{316}]\f$
40
 */
41
template <class FReal, int ORDER, class MatrixKernelClass>
42 43 44
unsigned int Compute(const MatrixKernelClass *const MatrixKernel, 
                                const FReal CellWidth,
                                FReal** &C);
45

46 47
//template <int ORDER>
//unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
48
//                                          FReal* &U,  FReal* &C, FReal* &B);
49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67


/**
 * @author Matthias Messner (matthias.messner@inria.fr)
 * @class FChebM2LHandler
 * Please read the license
 *
 * This class precomputes and compresses 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$.
 *
68 69 70
 * PB: BEWARE! Homogeneous matrix kernels do not support cell width extension
 * yet. Is it possible to find a reference width and a scale factor such that
 * only 1 set of M2L ops can be used for all levels?? 
71 72 73
 *
 * @tparam ORDER interpolation order \f$\ell\f$
 */
74
template <class FReal, int ORDER, class MatrixKernelClass, KERNEL_FUNCTION_TYPE TYPE> class FChebTensorialM2LHandler;
75

76 77
template <class FReal, int ORDER, class MatrixKernelClass>
class FChebTensorialM2LHandler<FReal, ORDER,MatrixKernelClass,HOMOGENEOUS> : FNoCopyable
78
{
79 80 81 82
    enum {order = ORDER,
          nnodes = TensorTraits<ORDER>::nnodes,
          ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field)
          ncmp = MatrixKernelClass::NCMP};
83

84
    FReal** C;
85

86
    const FReal CellWidthExtension; //<! extension of cells width
87

88
    unsigned int rank;   //<! truncation rank, satisfies @p epsilon
89 90


91
    static const std::string getFileName()
92 93 94 95
    {
        const char precision_type = (typeid(FReal)==typeid(double) ? 'd' : 'f');
        std::stringstream stream;
        stream << "m2l_k"<< MatrixKernelClass::getID() << "_" << precision_type
96
                     << "_o" << order << ".bin";
97 98
        return stream.str();
    }
99

100
    
101
public:
102 103 104
    FChebTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int, const FReal, const FReal inCellWidthExtension)
        : CellWidthExtension(inCellWidthExtension), 
          rank(0)
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
    {
        // measure time
        FTic time; time.tic();

        // allocate C
        C = new FReal*[ncmp];
        for (unsigned int d=0; d<ncmp; ++d) C[d]=nullptr;

        for (unsigned int d=0; d<ncmp; ++d)
            if (C[d]) throw std::runtime_error("Compressed M2L operator already set");

        // Compute matrix of interactions
        // reference cell width is arbitrarly set to 2. 
        // but it NEEDS to match the numerator of the scale factor in matrix kernel!
        // Therefore box width extension is not yet supported for homog kernels
        const FReal ReferenceCellWidth = FReal(2.);
121
        rank = Compute<order>(MatrixKernel, ReferenceCellWidth, 0., C);
122 123 124 125 126

        unsigned long sizeM2L = 343*ncmp*rank*rank*sizeof(FReal);


        // write info
127
        std::cout << "Compute and set full M2L operators (" << long(sizeM2L) << " B) in "
128 129 130 131 132 133 134 135
                  << time.tacAndElapsed() << "sec."   << std::endl;
    }

    ~FChebTensorialM2LHandler()
    {
        for (unsigned int d=0; d<ncmp; ++d)
            if (C[d] != nullptr) delete [] C[d];
    }
136

137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
    /**
     * @return rank of the SVD compressed M2L operators
     */
    unsigned int getRank() const {return rank;}

    /**
     * Compressed 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.
     *
     * @param[in] transfer transfer vector
     * @param[in] Y compressed multipole expansion
     * @param[out] X compressed local expansion
     * @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 unsigned int idx, const unsigned int , 
                const FReal scale, const unsigned int d,
                const FReal *const Y, FReal *const X) const
    {
        FBlas::gemva(rank, rank, scale, C[d] + idx*rank*rank, const_cast<FReal*>(Y), X);
    }

160 161 162
};


163 164
template <class FReal, int ORDER, class MatrixKernelClass>
class FChebTensorialM2LHandler<FReal,ORDER,MatrixKernelClass,NON_HOMOGENEOUS> : FNoCopyable
165
{
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
    enum {order = ORDER,
          nnodes = TensorTraits<ORDER>::nnodes,
          ninteractions = 316,// 7^3 - 3^3 (max num cells in far-field)
          ncmp = MatrixKernelClass::NCMP};

    // Tensorial MatrixKernel and homogeneity specific
    FReal*** C;

    const unsigned int TreeHeight; //<! number of levels
    const FReal RootCellWidth; //<! width of root cell
    const FReal CellWidthExtension; //<! extension of cells width

    unsigned int *rank;   //<! truncation rank, satisfies @p epsilon


181
    static const std::string getFileName()
182 183 184 185
    {
        const char precision_type = (typeid(FReal)==typeid(double) ? 'd' : 'f');
        std::stringstream stream;
        stream << "m2l_k"<< MatrixKernelClass::getID() << "_" << precision_type
186
                     << "_o" << order << ".bin";
187 188
        return stream.str();
    }
189

190 191
    
public:
192
    FChebTensorialM2LHandler(const MatrixKernelClass *const MatrixKernel, const unsigned int inTreeHeight, const FReal inRootCellWidth, const FReal inCellWidthExtension)
193 194
        : TreeHeight(inTreeHeight),
          RootCellWidth(inRootCellWidth),
195
          CellWidthExtension(inCellWidthExtension)
196 197 198 199 200 201 202 203 204 205 206 207 208 209
    {
        // measure time
        FTic time; time.tic();

        // allocate rank
        rank = new unsigned int[TreeHeight];

        // allocate C
        C = new FReal**[TreeHeight];
        for (unsigned int l=0; l<TreeHeight; ++l){ 
            C[l] = new FReal*[ncmp];
            for (unsigned int d=0; d<ncmp; ++d)
                C[l][d]=nullptr;
        }
210

211 212 213 214
        for (unsigned int l=0; l<TreeHeight; ++l) {
            for (unsigned int d=0; d<ncmp; ++d)
                if (C[l][d]) throw std::runtime_error("Compressed M2L operator already set");
        }
215

216 217 218 219 220 221
        // Compute matrix of interactions at each level !! (since non homog)
        FReal CellWidth = RootCellWidth / FReal(2.); // at level 1
        CellWidth /= FReal(2.);                      // at level 2
        rank[0]=rank[1]=0;
        for (unsigned int l=2; l<TreeHeight; ++l) {
            // compute m2l operator on extended cell
222
            rank[l] = Compute<FReal, order>(MatrixKernel, CellWidth, CellWidthExtension, C[l]);
223 224 225 226
            // update cell width
            CellWidth /= FReal(2.);                    // at level l+1 
        }
        unsigned long sizeM2L = (TreeHeight-2)*343*ncmp*rank[2]*rank[2]*sizeof(FReal);
227

228
        // write info
229
        std::cout << "Compute and Set full M2L operators of " << TreeHeight-2 << " levels ("<< long(sizeM2L/**1e-6*/) <<" Bytes) in "
230
                                << time.tacAndElapsed() << "sec."   << std::endl;
231 232
    }

233 234 235 236 237 238 239
    ~FChebTensorialM2LHandler()
    {
        if (rank != nullptr) delete [] rank;
        for (unsigned int l=0; l<TreeHeight; ++l) {
            for (unsigned int d=0; d<ncmp; ++d)
                if (C[l][d] != nullptr) delete [] C[l][d];
        }
240 241
    }

242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262
    /**
     * @return rank of the SVD compressed M2L operators
     */
    unsigned int getRank(unsigned int l = 2) const {return rank[l];}

    /**
     * Compressed 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.
     *
     * @param[in] transfer transfer vector
     * @param[in] Y compressed multipole expansion
     * @param[out] X compressed local expansion
     * @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 unsigned int idx, const unsigned int l, 
                const FReal, const unsigned int d,
                const FReal *const Y, FReal *const X) const
    {
        FBlas::gemva(rank[l], rank[l], 1., C[l][d] + idx*rank[l]*rank[l], const_cast<FReal*>(Y), X);
263 264 265 266
    }

};

267 268 269 270 271 272 273 274 275 276 277 278 279




//////////////////////////////////////////////////////////////////////
// definition ////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////






280
template <class FReal, int ORDER, class MatrixKernelClass>
281
unsigned int Compute(const MatrixKernelClass *const MatrixKernel, 
282
                                const FReal CellWidth, 
283 284
                                const FReal CellWidthExtension,
                                FReal** &C)
285
{
286 287 288 289 290 291 292 293
    // PB: need to redefine some constant since not function from m2lhandler class
    const unsigned int order = ORDER;
    const unsigned int nnodes = TensorTraits<ORDER>::nnodes;
    const unsigned int ninteractions = 316;
    const unsigned int ncmp = MatrixKernelClass::NCMP;

    // allocate memory and store compressed M2L operators
    for (unsigned int d=0; d<ncmp; ++d) 
294 295
    if (C[d]) throw std::runtime_error("Compressed M2L operators are already set");

296
    // interpolation points of source (Y) and target (X) cell
297
    FPoint<FReal> X[nnodes], Y[nnodes];
298 299
    // set roots of target cell (X)
    const FReal ExtendedCellWidth(CellWidth+CellWidthExtension);
300
    FChebTensor<FReal,order>::setRoots(FPoint<FReal>(0.,0.,0.), ExtendedCellWidth, X);
301

302 303 304 305
    // allocate memory and compute 316 m2l operators
    FReal** _C; 
    _C  = new FReal* [ncmp];
    for (unsigned int d=0; d<ncmp; ++d) 
306 307
    _C[d] = new FReal [nnodes*nnodes * ninteractions];

308 309 310 311 312 313
    unsigned int counter = 0;
    for (int i=-3; i<=3; ++i) {
        for (int j=-3; j<=3; ++j) {
            for (int k=-3; k<=3; ++k) {
                if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
                    // set roots of source cell (Y)
314 315
                    const FPoint<FReal> cy(CellWidth*FReal(i), CellWidth*FReal(j), CellWidth*FReal(k));
                    FChebTensor<FReal,order>::setRoots(cy, ExtendedCellWidth, Y);
316 317 318 319 320 321 322 323 324 325 326 327 328 329 330

                    // evaluate m2l operator
                    for (unsigned int n=0; n<nnodes; ++n)
                        for (unsigned int m=0; m<nnodes; ++m){

                            // Compute current M2L interaction (block matrix)
                            FReal* block; 
                            block = new FReal[ncmp]; 
                            MatrixKernel->evaluateBlock(X[m], Y[n], block);

                            // Copy block in C
                            for (unsigned int d=0; d<ncmp; ++d) 
                                _C[d][counter*nnodes*nnodes + n*nnodes + m] = block[d];

                            delete [] block;
331
              
332 333 334 335 336
                        }

                    // increment interaction counter
                    counter++;
                }
337 338
            }
        }
339 340 341
    }
    if (counter != ninteractions)
        throw std::runtime_error("Number of interactions must correspond to 316");
342
   
343

344
    // Copy M2L operators
345
    const unsigned int rank   = nnodes; //PB: dense Chebyshev
346
    if (!(rank>0)) throw std::runtime_error("Size must be larger than 0!");
347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367

    // store C
    counter = 0;
    for (unsigned int d=0; d<ncmp; ++d) 
    C[d] = new FReal [343 * rank*rank];
    for (int i=-3; i<=3; ++i)
        for (int j=-3; j<=3; ++j)
            for (int k=-3; k<=3; ++k) {
                const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
                if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
                    for (unsigned int d=0; d<ncmp; ++d) 
                        FBlas::copy(rank*rank, _C[d] + counter*rank*rank, C[d] + idx*rank*rank);
                    counter++;
                } else {
                    for (unsigned int d=0; d<ncmp; ++d) 
                        FBlas::setzero(rank*rank, C[d] + idx*rank*rank);
                }
            }
    if (counter != ninteractions)
        throw std::runtime_error("Number of interactions must correspond to 316");
    for (unsigned int d=0; d<ncmp; ++d) 
368
    delete [] _C[d];
369

370 371
    // return low rank
    return rank;
372 373 374 375 376 377 378 379
}




#endif // FCHEBTENSORIALM2LHANDLER_HPP

// [--END--]