Commit addf94a5 authored by BRAMAS Berenger's avatar BRAMAS Berenger
parents 0ab55fd3 cb425050
......@@ -154,6 +154,57 @@ public:
void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
const int NumSourceCells,
const int TreeLevel)
{
// permute and copy multipole expansion
memset(countExp, 0, sizeof(int) * 343);
for (unsigned int idx=0; idx<343; ++idx) {
if (SourceCells[idx]) {
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[idx]->getMultipole();
for (unsigned int n=0; n<nnodes; ++n)
mul[pvec[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);
for (unsigned int idx=0; idx<343; ++idx) {
if (SourceCells[idx]) {
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];
for (unsigned int n=0; n<nnodes; ++n)
LocalExpansion[n] += loc[pvec[n]];
}
}
}
/*
void M2L(CellClass* const FRestrict TargetCell,
const CellClass* SourceCells[343],
const int NumSourceCells,
......@@ -198,8 +249,8 @@ public:
}
}
}
*/
/*
void M2L(CellClass* const FRestrict TargetCell,
......
......@@ -99,8 +99,8 @@ class FChebSymmetries
const unsigned int getPermutationArrayAndIndex(const int i, const int j, const int k,
unsigned int permutation[nnodes]) 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 unsigned int qidx = getQuadIdx(i,j,k);
......
......@@ -79,11 +79,17 @@ FReal computeINFnorm(unsigned int N, FReal *const u, FReal *const v)
// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
const unsigned int ORDER = 5;
const FReal epsilon = FParameters::getValue(argc, argv, "-eps", FReal(1e-5));
const long NbPart = FParameters::getValue(argc, argv, "-num", 100000);
const unsigned int TreeHeight = FParameters::getValue(argc, argv, "-h", 5);
const unsigned int ORDER = 3;
const FReal epsilon = FParameters::getValue(argc, argv, "-eps", FReal(1e-3));
const long NbPart = FParameters::getValue(argc, argv, "-num", 4000000);
const unsigned int TreeHeight = FParameters::getValue(argc, argv, "-h", 7);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, "-sh", 2);
const unsigned int NbThreads = FParameters::getValue(argc, argv, "-t", 1);
omp_set_num_threads(NbThreads);
std::cout << "Using " << omp_get_max_threads() << " threads." << std::endl;
const FReal Width = 10.;
// init random fun
......
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