Commit 06ad1f29 authored by messner's avatar messner
Browse files

added use of symmetries in the far-field; added swap function to FMemUtils


git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@474 2616d619-271b-44dc-8df4-d4a8f33a7222
parent d759f799
#ifndef FCHEBSYMMETRIES_HPP
#define FCHEBSYMMETRIES_HPP
#include "../Utils/FBlas.hpp"
#include "./FChebTensor.hpp"
/**
* @author Matthias Messner (matthias.matthias@inria.fr)
* Please read the license
*/
/**
* @class FChebSymmetries
*
* The class @p FChebSymmetries exploits all symmetries
*/
template <int ORDER>
class FChebSymmetries
{
enum {nnodes = TensorTraits<ORDER>::nnodes};
unsigned int quads[8][nnodes];
unsigned int cones[8][nnodes];
unsigned int perms[8][3];
public:
FChebSymmetries()
{
unsigned int evn[ORDER]; unsigned int odd[ORDER];
for (unsigned int o=0; o<ORDER; ++o) {evn[o] = o; odd[o] = ORDER-1 - o;}
for (unsigned int i=0; i<ORDER; ++i) {
for (unsigned int j=0; j<ORDER; ++j) {
for (unsigned int k=0; k<ORDER; ++k) {
const unsigned int index = k*ORDER*ORDER + j*ORDER + i;
// global axis parallel symmetries (8 quads) ///////////
quads[0][index] = evn[k]*ORDER*ORDER + evn[j]*ORDER + evn[i]; // - - -
quads[1][index] = evn[k]*ORDER*ORDER + evn[j]*ORDER + odd[i]; // - - +
quads[2][index] = evn[k]*ORDER*ORDER + odd[j]*ORDER + evn[i]; // - + -
quads[3][index] = evn[k]*ORDER*ORDER + odd[j]*ORDER + odd[i]; // - + +
quads[4][index] = odd[k]*ORDER*ORDER + evn[j]*ORDER + evn[i]; // + - -
quads[5][index] = odd[k]*ORDER*ORDER + evn[j]*ORDER + odd[i]; // + - +
quads[6][index] = odd[k]*ORDER*ORDER + odd[j]*ORDER + evn[i]; // + + -
quads[7][index] = odd[k]*ORDER*ORDER + odd[j]*ORDER + odd[i]; // + + +
// diagonal symmetries (j>i)(k>i)(k>j) /////////////////////
cones[0][index] = evn[k]*ORDER*ORDER + evn[j]*ORDER + evn[i]; // (0) 000
cones[1][index] = evn[j]*ORDER*ORDER + evn[k]*ORDER + evn[i]; // (1) 001
// cones[2] does not exist
cones[3][index] = evn[j]*ORDER*ORDER + evn[i]*ORDER + evn[k]; // (3) 011
cones[4][index] = evn[k]*ORDER*ORDER + evn[i]*ORDER + evn[j]; // (4) 100
// cones[5] does not exist
cones[6][index] = evn[i]*ORDER*ORDER + evn[k]*ORDER + evn[j]; // (6) 110
cones[7][index] = evn[i]*ORDER*ORDER + evn[j]*ORDER + evn[k]; // (7) 111
}
}
}
// permutation of interaction indices (already absolute value)
perms[0][0] = 0; perms[0][1] = 1; perms[0][2] = 2;
perms[1][0] = 0; perms[1][1] = 2; perms[1][2] = 1;
// perms[2] does not exist
perms[3][0] = 2; perms[3][1] = 0; perms[3][2] = 1;
perms[4][0] = 1; perms[4][1] = 0; perms[4][2] = 2;
// perms[5] does not exist
perms[6][0] = 1; perms[6][1] = 2; perms[6][2] = 0;
perms[7][0] = 2; perms[7][1] = 1; perms[7][2] = 0;
}
const unsigned int getPermutationArrayAndIndex(const int i, const int j, const int k,
unsigned int permutation[nnodes]) const
{
// find right quadrant index (if < 0 then 0, else 1)
const int si = ((unsigned int)i >> (sizeof(int)*CHAR_BIT - 1));
const int sj = ((unsigned int)j >> (sizeof(int)*CHAR_BIT - 2)) & 2;
const int sk = ((unsigned int)k >> (sizeof(int)*CHAR_BIT - 3)) & 4;
const int qidx = (sk | sj | si);
// get absolut values of i,j,k
const int imask = i >> (sizeof(int) * CHAR_BIT - 1);
const int jmask = j >> (sizeof(int) * CHAR_BIT - 1);
const int kmask = k >> (sizeof(int) * CHAR_BIT - 1);
const unsigned int u[3] = {(i + imask) ^ imask,
(j + jmask) ^ jmask,
(k + kmask) ^ kmask};
// find right cone index
const int q0 = (u[1]>u[0])<<2;
const int q1 = (u[2]>u[0])<<1;
const int q2 = (u[2]>u[1]);
const int cidx = (q2 | q1 | q0);
// set permutation array (combination of quadrant and cone) //////
for (unsigned int n=0; n<nnodes; ++n) permutation[n] = cones[cidx][quads[qidx][n]];
// set permutation index /////////////////////////////////////////
return (u[perms[cidx][0]]+3)*7*7 + (u[perms[cidx][1]]+3)*7 + (u[perms[cidx][2]]+3);
}
static void permuteMatrix(const unsigned int perm[nnodes], FReal *const Matrix)
{
// 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);
}
};
#endif
//// if (uj>ui) FMemUtils::swap(uj,ui);
//// if (uk>ui) FMemUtils::swap(uk,ui);
//// if (uk>uj) FMemUtils::swap(uk,uj);
//// const unsigned int *const getPermQuadrant(const int i, const int j, const int k) const
//// {
//// // if < 0 then 0, else 1
//// const int si = ((unsigned int)i >> (sizeof(int)*CHAR_BIT - 1));
//// const int sj = ((unsigned int)j >> (sizeof(int)*CHAR_BIT - 2)) & 2;
//// const int sk = ((unsigned int)k >> (sizeof(int)*CHAR_BIT - 3)) & 4;
//// const int qidx = sk | sj | si;
//// return quads[qidx];
//// }
//// const void permuteMatrix(const int i, const int j, const int k, FReal *const Matrix) const
//// { permuteMatrix(getPermQuadrant(i, j, k), Matrix); }
......@@ -25,7 +25,7 @@ namespace FMemUtils {
static const FSize MaxSize_t = UINT_MAX; //std::numeric_limits<std::size_t>::max();
/** memcpy */
static void* memcpy(void* const dest, const void* const source, const FSize nbBytes){
void* memcpy(void* const dest, const void* const source, const FSize nbBytes){
if( nbBytes < MaxSize_t){
return ::memcpy(dest, source, size_t(nbBytes));
}
......@@ -45,7 +45,7 @@ namespace FMemUtils {
}
/** memset */
static void* memset(void* const dest, const int val, const FSize nbBytes){
void* memset(void* const dest, const int val, const FSize nbBytes){
if( nbBytes < MaxSize_t){
return ::memset(dest, val, size_t(nbBytes));
}
......@@ -85,6 +85,14 @@ namespace FMemUtils {
(*dest++) = source;
}
}
/** swap values from a and b*/
template <class TypeClass>
void swap(TypeClass& a, TypeClass& b){
TypeClass c(a);
a=b;
b=c;
}
/** Delete all */
template <class TypeClass>
......
......@@ -50,7 +50,7 @@ int main(int argc, char* argv[])
// constants
const FReal epsilon = FReal(atof(argv[1]));
const unsigned int order = 4;
const unsigned int order = 9;
// number of interpolation points per cell
const unsigned int nnodes = TensorTraits<order>::nnodes;
......@@ -120,9 +120,9 @@ int main(int argc, char* argv[])
time.tic();
// compute 316 m2l operators
C1 = new FReal [nnodes*nnodes];
const unsigned int i = 2;
const unsigned int j = 0;
const unsigned int k = 0;
const unsigned int i = 3;
const unsigned int j = 3;
const unsigned int k = 3;
const F3DPosition cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
FChebTensor<order>::setRoots(cy, FReal(2.), Y);
// evaluate m2l operator
......
// ===================================================================================
// 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>
#include "../Src/Utils/FTic.hpp"
#include "../Src/Chebyshev/FChebTensor.hpp"
#include "../Src/Chebyshev/FChebMatrixKernel.hpp"
#include "../Src/Chebyshev/FChebSymmetries.hpp"
int main(int argc, char* argv[])
{
// start timer /////////////////////////////////
FTic time;
// define set matrix kernel
typedef FChebMatrixKernelR MatrixKernelClass;
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
F3DPosition X[nnodes], Y[nnodes];
// set roots of target cell (X)
FChebTensor<order>::setRoots(F3DPosition(0.,0.,0.), FReal(2.), X);
// 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] = NULL;
{
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;
const F3DPosition cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
FChebTensor<order>::setRoots(cy, FReal(2.), Y);
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
FChebSymmetries<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) {
const F3DPosition cy(FReal(2.*i), FReal(2.*j), FReal(2.*k));
FChebTensor<order>::setRoots(cy, FReal(2.), Y);
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);
permuter.permuteMatrix(perm, K0);
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;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment