testChebSymmetries.cpp 5.46 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
// ===================================================================================
// Ce LOGICIEL "ScalFmm" est couvert par le copyright Inria 20xx-2012.
// Inria détient tous les droits de propriété sur le LOGICIEL, et souhaite que
// la communauté scientifique l'utilise afin de le tester et de l'évaluer.
// Inria donne gracieusement le droit d'utiliser ce LOGICIEL. Toute utilisation
// dans un but lucratif ou à des fins commerciales est interdite sauf autorisation
// expresse et préalable d'Inria.
// Toute utilisation hors des limites précisées ci-dessus et réalisée sans l'accord
// expresse préalable d'Inria constituerait donc le délit de contrefaçon.
// Le LOGICIEL étant un produit en cours de développement, Inria ne saurait assurer
// aucune responsabilité et notamment en aucune manière et en aucun cas, être tenu
// de répondre d'éventuels dommages directs ou indirects subits par l'utilisateur.
// Tout utilisateur du LOGICIEL s'engage à communiquer à Inria ses remarques
// relatives à l'usage du LOGICIEL
// ===================================================================================

// ==== CMAKE =====
// @FUSE_BLAS
// ================


// system includes
#include <iostream>
#include <stdexcept>
#include <climits>



29 30
#include "Utils/FTic.hpp"
#include "Utils/FBlas.hpp"
31

32 33 34
#include "Kernels/Chebyshev/FChebTensor.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp"
#include "Kernels/Interpolation/FInterpSymmetries.hpp"
35

36
#include "Utils/FParameterNames.hpp"
37

38
template < class FReal, int ORDER>
Matthias Messner's avatar
Matthias Messner committed
39 40
void permuteMatrix(const unsigned int perm[ORDER*ORDER*ORDER], FReal *const Matrix)
{
41 42 43 44 45 46 47 48 49
    const unsigned int nnodes = ORDER*ORDER*ORDER;
    // allocate temporary memory for matrix
    FReal temp[nnodes*nnodes];
    // permute rows
    for (unsigned int r=0; r<nnodes; ++r)
        FBlas::copy(nnodes, Matrix+r, nnodes, temp+perm[r], nnodes);
    // permute columns
    for (unsigned int c=0; c<nnodes; ++c)
        FBlas::copy(nnodes, temp + c * nnodes, Matrix + perm[c] * nnodes);
Matthias Messner's avatar
Matthias Messner committed
50 51
}

52 53 54 55 56 57




int main(int argc, char* argv[])
{
58 59 60 61 62 63
    FHelpDescribeAndExit(argc, argv, "Test Chebyshev symmetries.");

    // start timer /////////////////////////////////
    FTic time;

    // define set matrix kernel
64 65
    typedef double FReal;
    typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
66 67 68 69 70 71 72 73 74 75
    MatrixKernelClass MatrixKernel;

    // constants
    //const FReal epsilon     = FReal(atof(argv[1]));
    const unsigned int order = 5;

    // number of interpolation points per cell
    const unsigned int nnodes = TensorTraits<order>::nnodes;

    // interpolation points of source (Y) and target (X) cell
76
    FPoint<FReal> X[nnodes], Y[nnodes];
77
    // set roots of target cell (X)
78
    FChebTensor<FReal,order>::setRoots(FPoint<FReal>(0.,0.,0.), FReal(2.), X);
79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95

    // allocate 343 pointers to K, but only 16 are actually filled
    FReal** K = new FReal* [343];
    for (unsigned int t=0; t<343; ++t) K[t] = nullptr;

    {
        unsigned int counter = 0;

        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);
                    K[idx] = new FReal [nnodes*nnodes];

                    //std::cout << i << "," << j << "," << k << "\t" << idx << std::endl;

96 97
                    const FPoint<FReal> cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
                    FChebTensor<FReal,order>::setRoots(cy, FReal(2.), Y);
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
                    for (unsigned int n=0; n<nnodes; ++n)
                        for (unsigned int m=0; m<nnodes; ++m)
                            K[idx][n*nnodes + m] = MatrixKernel.evaluate(X[m], Y[n]);

                    counter++;
                }
            }
        }

        std::cout << "num interactions = " << counter << std::endl;
    }

    // max difference
    FReal maxdiff(0.);

    // permuter
    FInterpSymmetries<order> permuter;

    // permutation vector
    unsigned int perm[nnodes];

    FReal* K0 = new FReal [nnodes*nnodes];

    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) {

128 129
                    const FPoint<FReal> cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
                    FChebTensor<FReal,order>::setRoots(cy, FReal(2.), Y);
130 131 132 133 134 135
                    for (unsigned int n=0; n<nnodes; ++n)
                        for (unsigned int m=0; m<nnodes; ++m)
                            K0[n*nnodes + m] = MatrixKernel.evaluate(X[m], Y[n]);

                    // permute
                    const unsigned int pidx = permuter.getPermutationArrayAndIndex(i, j, k, perm);
136
                    permuteMatrix<FReal,order>(perm, K0);
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

                    if (K[pidx]==NULL) std::cout << " - not existing index " << pidx << std::endl;

                    FReal mdiff(0.);
                    for (unsigned int n=0; n<nnodes*nnodes; ++n) {
                        FReal diff(K0[n] - K[pidx][n]);
                        if (FMath::Abs(diff)>mdiff) mdiff = diff;
                    }
                    if (FMath::Abs(mdiff)>maxdiff) maxdiff = FMath::Abs(mdiff);

                    if (mdiff > 1e-15) exit(-1);
                    counter++;

                }
            }
        }
    }

    std::cout << "Max error = " << maxdiff << " of counter = " << counter << std::endl;

    for (unsigned int t=0; t<343; ++t) if (K[t]!=NULL) delete [] K[t];
    delete [] K;
    delete [] K0;

    return 0;
162 163 164 165
}