FChebSymKernel.hpp 19.2 KB
Newer Older
1 2
#ifndef FCHEBSYMKERNEL_HPP
#define FCHEBSYMKERNEL_HPP
BRAMAS Berenger's avatar
BRAMAS Berenger committed
3
// ===================================================================================
4
// Copyright ScalFmm 2011 INRIA,
5 6 7 8 9 10 11 12 13 14 15
// 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".
BRAMAS Berenger's avatar
BRAMAS Berenger committed
16
// ===================================================================================
17 18

#include "Utils/FGlobal.hpp"
19

20 21 22 23 24
#include "Utils/FSmartPointer.hpp"

#include "FAbstractChebKernel.hpp"
#include "FChebInterpolator.hpp"
#include "FChebSymM2LHandler.hpp"
25 26 27 28

class FTreeCoordinate;


29
// for verbosity only!!!
30 31
//#define COUNT_BLOCKED_INTERACTIONS

32 33
// if timings should be logged
//#define LOG_TIMINGS
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

/**
 * @author Matthias Messner(matthias.messner@inria.fr)
 * @class FChebSymKernel
 * @brief
 * Please read the license
 *
 * This kernels implement the Chebyshev interpolation based FMM operators
 * exploiting the symmetries in the far-field. 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
 */
51
template < class FReal, class CellClass, class ContainerClass,   class MatrixKernelClass, int ORDER, int NVALS = 1>
52
class FChebSymKernel
53
        : public FAbstractChebKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
54
{
55 56
    typedef FAbstractChebKernel<FReal, CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> AbstractBaseClass;
    typedef SymmetryHandler<FReal, ORDER, MatrixKernelClass::Type> SymmetryHandlerClass;
57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
    enum {nnodes = AbstractBaseClass::nnodes};

    /// Needed for P2P and M2L operators
    const MatrixKernelClass *const MatrixKernel;

    /// Needed for handling all symmetries
    const FSmartPointer<SymmetryHandlerClass,FSmartPointerMemory> SymHandler;

    // permuted local and multipole expansions
    FReal** Loc;
    FReal** Mul;
    unsigned int* countExp;


    
    /**
     * Allocate memory for storing locally permuted mulipole and local expansions
     */
    void allocateMemoryForPermutedExpansions()
    {
        assert(Loc==nullptr && Mul==nullptr && countExp==nullptr);
        Loc = new FReal* [343];
        Mul = new FReal* [343];
        countExp = new unsigned int [343];

        // set all 343 to nullptr
        for (unsigned int idx=0; idx<343; ++idx) {
            Mul[idx] = Loc[idx] = nullptr;
        }

        // init only 16 of 343 possible translations due to symmetries
        for (int i=2; i<=3; ++i)
            for (int j=0; j<=i; ++j)
                for (int k=0; k<=j; ++k) {
                    const unsigned int idx = (i+3)*7*7 + (j+3)*7 + (k+3);
                    assert(Mul[idx]==nullptr || Loc[idx]==nullptr);
                    Mul[idx] = new FReal [24 * nnodes];
                    Loc[idx] = new FReal [24 * nnodes];
                }
    }
    
98

99
#ifdef LOG_TIMINGS
100 101
    FTic time;
    FReal t_m2l_1, t_m2l_2, t_m2l_3;
102
#endif
103 104

public:
105 106 107 108 109 110 111
    /**
     * 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).
     */
    FChebSymKernel(const int inTreeHeight,
                   const FReal inBoxWidth,
112
                   const FPoint<FReal>& inBoxCenter,
113 114
                   const MatrixKernelClass *const inMatrixKernel,
                   const FReal Epsilon)
115 116 117 118
        : AbstractBaseClass(inTreeHeight, inBoxWidth, inBoxCenter),
          MatrixKernel(inMatrixKernel),
          SymHandler(new SymmetryHandlerClass(MatrixKernel, Epsilon, inBoxWidth, inTreeHeight)),
          Loc(nullptr), Mul(nullptr), countExp(nullptr)
119 120
    {
        this->allocateMemoryForPermutedExpansions();
121
#ifdef LOG_TIMINGS
122 123 124
        t_m2l_1 = FReal(0.);
        t_m2l_2 = FReal(0.);
        t_m2l_3 = FReal(0.);
125
#endif
126 127 128 129 130 131 132 133 134
    }
    
    /**
     * 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).
     */
    FChebSymKernel(const int inTreeHeight,
                   const FReal inBoxWidth,
135
                   const FPoint<FReal>& inBoxCenter,
136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162
                   const MatrixKernelClass *const inMatrixKernel) :FChebSymKernel(inTreeHeight,inBoxWidth, inBoxCenter,inMatrixKernel,FMath::pow(10.0,static_cast<FReal>(-ORDER)))
    {}

    

    /** Copy constructor */
    FChebSymKernel(const FChebSymKernel& other)
        : AbstractBaseClass(other),
          MatrixKernel(other.MatrixKernel),
          SymHandler(other.SymHandler),
          Loc(nullptr), Mul(nullptr), countExp(nullptr)
    {
        this->allocateMemoryForPermutedExpansions();
    }



    /** Destructor */
    ~FChebSymKernel()
    {
        for (unsigned int t=0; t<343; ++t) {
            if (Loc[t]!=nullptr) delete [] Loc[t];
            if (Mul[t]!=nullptr) delete [] Mul[t];
        }
        if (Loc!=nullptr)      delete [] Loc;
        if (Mul!=nullptr)      delete [] Mul;
        if (countExp!=nullptr) delete [] countExp;
163 164

#ifdef LOG_TIMINGS
165 166 167 168
        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;
169
#endif
170
    }
171 172


173
    const SymmetryHandlerClass * getPtrToSymHandler() const
174 175
    {   return SymHandler.getPtr(); }
    
176

177

178
    void P2M(CellClass* const LeafCell,
179
             const ContainerClass* const SourceParticles/*, const int level = AbstractBaseClass::TreeHeight*/)
180 181
    {
        // apply Sy
182
        const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
183 184
        AbstractBaseClass::Interpolator->applyP2M(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                  LeafCell->getMultipole(0), SourceParticles);
185
    }
186 187 188



189 190 191 192
    void M2M(CellClass* const FRestrict ParentCell,
             const CellClass*const FRestrict *const FRestrict ChildCells,
             const int /*TreeLevel*/)
    {
193
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
194 195
            // Reset the Parent expansion to zero
            //   FBlas::scal(nnodes*2, FReal(0.), ParentCell->getMultipole(idxRhs));
196
            for (unsigned int ChildIndex=0; ChildIndex < 8; ++ChildIndex){
197
                // apply Sy
198
                if (ChildCells[ChildIndex]){
199 200 201 202
                    AbstractBaseClass::Interpolator->applyM2M(ChildIndex, ChildCells[ChildIndex]->getMultipole(idxRhs), ParentCell->getMultipole(idxRhs));
                }
            }
        }
203
    }
204 205 206



207 208
    void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
             const int neighborPositions[], const int inSize, const int TreeLevel)  override {
209
#ifdef LOG_TIMINGS
210
        time.tic();
211
#endif
212 213 214
        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
            // permute and copy multipole expansion
            memset(countExp, 0, sizeof(int) * 343);
215 216 217 218 219 220 221 222 223
            for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
                const int idx = neighborPositions[idxExistingNeigh];
                const unsigned int pidx = SymHandler->pindices[idx];
                const unsigned int count = (countExp[pidx])++;
                FReal *const mul = Mul[pidx] + count*nnodes;
                const unsigned int *const pvec = SymHandler->pvectors[idx];
                const FReal *const MultiExp = SourceCells[idxExistingNeigh]->getMultipole(idxRhs);

                /*
224 225 226 227 228
                    // no loop unrolling
                    for (unsigned int n=0; n<nnodes; ++n)
                        mul[pvec[n]] = MultiExp[n];
                    */

229 230 231 232 233 234 235 236 237
                // explicit loop unrolling
                for (unsigned int n=0; n<nnodes/4 * 4; n+=4) {
                    mul[pvec[n  ]] = MultiExp[n];
                    mul[pvec[n+1]] = MultiExp[n+1];
                    mul[pvec[n+2]] = MultiExp[n+2];
                    mul[pvec[n+3]] = MultiExp[n+3];
                }
                for (unsigned int n=nnodes/4 * 4; n<nnodes; ++n){
                    mul[pvec[n]] = MultiExp[n];
238 239 240
                }
            }

241
#ifdef LOG_TIMINGS
242
            t_m2l_1 += time.tacAndElapsed();
243
#endif
244 245


246
#ifdef COUNT_BLOCKED_INTERACTIONS ////////////////////////////////////
247 248 249 250 251 252 253 254
            unsigned int count_lidx = 0;
            unsigned int count_interactions = 0;
            for (unsigned int idx=0; idx<343; ++idx)
                count_interactions += countExp[idx];
            if (count_interactions==189) {
                for (unsigned int idx=0; idx<343; ++idx) {
                    if (countExp[idx])
                        std::cout << "gidx = " << idx << " gives lidx = " << count_lidx++ << " and has "
255
                                  << countExp[idx] << " interactions" << std::endl;
256 257 258
                }
                std::cout << std::endl;
            }
259
#endif ///////////////////////////////////////////////////////////////
260 261


262
#ifdef LOG_TIMINGS
263
            time.tic();
264
#endif
265 266 267

            // multiply (mat-mat-mul)
            FReal Compressed [nnodes * 24];
268
            const FReal scale = MatrixKernel->getScaleFactor(AbstractBaseClass::BoxWidth, TreeLevel);
269 270 271 272 273 274
            for (unsigned int pidx=0; pidx<343; ++pidx) {
                const unsigned int count = countExp[pidx];
                if (count) {
                    const unsigned int rank = SymHandler->getLowRank(TreeLevel, pidx);
                    // rank * count * (2*nnodes-1) flops
                    FBlas::gemtm(nnodes, rank, count, FReal(1.),
275 276
                                 const_cast<FReal*>(SymHandler->getK(TreeLevel, pidx))+rank*nnodes,
                                 nnodes, Mul[pidx], nnodes, Compressed, rank);
277 278
                    // nnodes *count * (2*rank-1) flops
                    FBlas::gemm( nnodes, rank, count, scale,
279 280
                                 const_cast<FReal*>(SymHandler->getK(TreeLevel, pidx)),
                                 nnodes, Compressed, rank, Loc[pidx], nnodes);
281 282 283
                }
            }

284
#ifdef LOG_TIMINGS
285
            t_m2l_2 += time.tacAndElapsed();
286
#endif
287

288
#ifdef LOG_TIMINGS
289
            time.tic();
290
#endif
291 292 293 294

            // permute and add contribution to local expansions
            FReal *const LocalExpansion = TargetCell->getLocal(idxRhs);
            memset(countExp, 0, sizeof(int) * 343);
295 296 297 298 299 300
            for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
                const int idx = neighborPositions[idxExistingNeigh];
                const unsigned int pidx = SymHandler->pindices[idx];
                const unsigned int count = (countExp[pidx])++;
                const FReal *const loc = Loc[pidx] + count*nnodes;
                const unsigned int *const pvec = SymHandler->pvectors[idx];
301

302 303 304 305 306 307 308 309 310 311 312 313
                /*
                // no loop unrolling
                for (unsigned int n=0; n<nnodes; ++n)
                    LocalExpansion[n] += loc[pvec[n]];
                */

                // explicit loop unrolling
                for (unsigned int n=0; n<nnodes/4 * 4; n+=4) {
                    LocalExpansion[n  ] += loc[pvec[n  ]];
                    LocalExpansion[n+1] += loc[pvec[n+1]];
                    LocalExpansion[n+2] += loc[pvec[n+2]];
                    LocalExpansion[n+3] += loc[pvec[n+3]];
314
                }
315 316 317
                for (unsigned int n=nnodes/4 * 4; n<nnodes; ++n)
                    LocalExpansion[n] += loc[pvec[n]];

318 319
            }

320
#ifdef LOG_TIMINGS
321
            t_m2l_3 += time.tacAndElapsed();
322
#endif
323
        }
324

325 326 327 328 329
    }



    /*
330 331
    void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
             const int neighborPositions[], const int inSize, const int TreeLevel)  override {
332 333
        // permute and copy multipole expansion
        memset(countExp, 0, sizeof(int) * 343);
334 335
            for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
                const int idx = neighborPositions[idxExistingNeigh];
336 337
                const unsigned int pidx = SymHandler->pindices[idx];
                const unsigned int count = (countExp[pidx])++;
338
                const FReal *const MultiExp = SourceCells[idxExistingNeigh]->getMultipole();
339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360
                for (unsigned int n=0; n<nnodes; ++n)
                    Mul[pidx][count*nnodes + SymHandler->pvectors[idx][n]] = MultiExp[n];
            }

        // multiply (mat-mat-mul)
        FReal Compressed [nnodes * 30];
        const FReal CellWidth(AbstractBaseClass::BoxWidth / FReal(FMath::pow(2, TreeLevel)));
        const FReal scale(AbstractBaseClass::MatrixKernel->getScaleFactor(CellWidth));
        for (unsigned int pidx=0; pidx<343; ++pidx) {
            const unsigned int count = countExp[pidx];
            if (count) {
                const unsigned int rank = SymHandler->LowRank[pidx];
                FBlas::gemtm(nnodes, rank, count, FReal(1.),
                                         SymHandler->K[pidx]+rank*nnodes, nnodes, Mul[pidx], nnodes, Compressed, rank);
                FBlas::gemm( nnodes, rank, count, scale,
                                         SymHandler->K[pidx], nnodes, Compressed, rank, Loc[pidx], nnodes);
            }
        }

        // permute and add contribution to local expansions
        FReal *const LocalExpansion = TargetCell->getLocal();
        memset(countExp, 0, sizeof(int) * 343);
361 362
            for(int idxExistingNeigh = 0 ; idxExistingNeigh < inSize ; ++idxExistingNeigh){
                const int idx = neighborPositions[idxExistingNeigh];
363 364 365 366 367 368 369 370 371 372
                const unsigned int pidx = SymHandler->pindices[idx];
                const unsigned int count = (countExp[pidx])++;
                for (unsigned int n=0; n<nnodes; ++n)
                    LocalExpansion[n] += Loc[pidx][count*nnodes + SymHandler->pvectors[idx][n]];
            }

    }
    */

    /*
373 374
    void M2L(CellClass* const FRestrict TargetCell, const CellClass* SourceCells[],
             const int neighborPositions[], const int inSize, const int TreeLevel)  override {
375 376 377 378 379 380 381 382 383 384
        FReal *const LocalExpansion = TargetCell->getLocal();
        const FReal CellWidth(BoxWidth / FReal(FMath::pow(2, TreeLevel)));
        const FReal scale(MatrixKernel->getScaleFactor(CellWidth));

        FReal PermLocalExp[nnodes];
        FReal PermMultiExp[nnodes];
        FReal Compressed[nnodes];
        for (int i=-3; i<=3; ++i) {
        for (int j=-3; j<=3; ++j) {
        for (int k=-3; k<=3; ++k) {
385

386
        const unsigned int idx = ((i+3) * 7 + (j+3)) * 7 + (k+3);
387

388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416
        if (SourceCells[idx]) {
        const FReal *const MultiExp = SourceCells[idx]->getMultipole();

        // permute
        for (unsigned int n=0; n<nnodes; ++n) PermMultiExp[pvectors[idx][n]] = MultiExp[n];

        // mat-vec-mult (svd)
        assert(K[pindices[idx]]!=nullptr);
        const int rank = LowRank[pindices[idx]];
        FBlas::gemtv(nnodes, rank, FReal(1.), K[pindices[idx]]+rank*nnodes, PermMultiExp, Compressed);
        FBlas::gemv( nnodes, rank, scale, K[pindices[idx]], Compressed, PermLocalExp);

        // permute
        for (unsigned int n=0; n<nnodes; ++n) LocalExpansion[n] += PermLocalExp[pvectors[idx][n]];
        }

        }
        }
        }

        }
    */



    void L2L(const CellClass* const FRestrict ParentCell,
             CellClass* FRestrict *const FRestrict ChildCells,
             const int /*TreeLevel*/)
    {
417 418 419 420 421 422 423 424
        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));
                }
            }
        }
425
    }
426 427 428



429 430 431
    void L2P(const CellClass* const LeafCell,
             ContainerClass* const TargetParticles)
    {
432
        const FPoint<FReal> LeafCellCenter(AbstractBaseClass::getLeafCellCenter(LeafCell->getCoordinate()));
433 434 435 436 437 438 439 440 441 442
        //      // a) apply Sx
        //      AbstractBaseClass::Interpolator->applyL2P(LeafCellCenter,
        //                                                AbstractBaseClass::BoxWidthLeaf,
        //                                                LeafCell->getLocal(),
        //                                                TargetParticles);
        //      // b) apply Px (grad Sx)
        //      AbstractBaseClass::Interpolator->applyL2PGradient(LeafCellCenter,
        //                                                        AbstractBaseClass::BoxWidthLeaf,
        //                                                        LeafCell->getLocal(),
        //                                                        TargetParticles);
443

444 445 446 447
        // c) apply Sx and Px (grad Sx)
        AbstractBaseClass::Interpolator->applyL2PTotal(LeafCellCenter, AbstractBaseClass::BoxWidthLeaf,
                                                       LeafCell->getLocal(0), TargetParticles);

448
    }
449

450
    void P2P(const FTreeCoordinate& inPosition,
451
             ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict inSources,
452 453
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
454 455 456 457 458 459 460 461 462
        if(inTargets == inSources){
            P2POuter(inPosition, inTargets, inNeighbors, neighborPositions, inSize);
            DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PInner(inTargets,MatrixKernel);
        }
        else{
            const ContainerClass* const srcPtr[1] = {inSources};
            DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,srcPtr,1,MatrixKernel);
            DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,inSize,MatrixKernel);
        }
463 464 465 466 467 468
    }

    void P2POuter(const FTreeCoordinate& /*inLeafPosition*/,
             ContainerClass* const FRestrict inTargets,
             ContainerClass* const inNeighbors[], const int neighborPositions[],
             const int inSize) override {
469 470 471 472 473 474
        int nbNeighborsToCompute = 0;
        while(nbNeighborsToCompute < inSize
              && neighborPositions[nbNeighborsToCompute] < 14){
            nbNeighborsToCompute += 1;
        }
        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2P(inTargets,inNeighbors,nbNeighborsToCompute,MatrixKernel);
475 476 477 478 479
    }


    void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                   ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
480
                   const ContainerClass* const inNeighbors[], const int /*neighborPositions*/[],
481 482
                   const int inSize) override {
        DirectInteractionComputer<FReal, MatrixKernelClass::NCMP, NVALS>::P2PRemote(inTargets,inNeighbors,inSize,MatrixKernel);
483
    }
484 485 486 487 488 489 490 491 492 493 494 495 496

};








#endif //FCHEBSYMKERNELS_HPP

// [--END--]