Commit 655fc192 authored by messner's avatar messner

update of weighting functionality of M2L operators (WIP)


git-svn-id: svn+ssh://scm.gforge.inria.fr/svn/scalfmm/scalfmm/trunk@445 2616d619-271b-44dc-8df4-d4a8f33a7222
parent 4e1c99bf
......@@ -85,12 +85,10 @@ public:
Epsilon(inEpsilon)
{
// read precomputed compressed m2l operators from binary file
M2LHandler->ReadFromBinaryFileAndSet();
//M2LHandler->ReadFromBinaryFileAndSet();
M2LHandler->ComputeAndCompressAndSet();
}
// /** Default destructor */
// ~FChebKernels() {}
void P2M(CellClass* const LeafCell,
const ContainerClass* const SourceParticles)
......
......@@ -86,7 +86,8 @@ public:
if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
rank = ComputeAndCompress(epsilon, U, C, B);
// write info
std::cout << "Compressed and set M2L operators in " << time.tacAndElapsed() << "sec." << std::endl;
std::cout << "Compressed and set M2L operators (" << rank << ") in "
<< time.tacAndElapsed() << "sec." << std::endl;
}
/**
......@@ -205,6 +206,9 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
FReal* &C,
FReal* &B)
{
// allocate memory and store compressed M2L operators
if (U||C||B) throw std::runtime_error("Compressed M2L operators are already set");
// interpolation points of source (Y) and target (X) cell
F3DPosition X[nnodes], Y[nnodes];
// set roots of target cell (X)
......@@ -238,12 +242,22 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
if (counter != ninteractions)
throw std::runtime_error("Number of interactions must correspond to 316");
//////////////////////////////////////////////////////////
FReal weights[nnodes];
FChebTensor<order>::setRootOfWeights(weights);
for (unsigned int i=0; i<316; ++i)
for (unsigned int n=0; n<nnodes; ++n) {
FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n, nnodes); // scale rows
FBlas::scal(nnodes, weights[n], _C+i*nnodes*nnodes + n * nnodes); // scale cols
}
//////////////////////////////////////////////////////////
// svd compression of M2L
const unsigned int rank = Compress<ORDER>(epsilon, ninteractions, _U, _C, _B);
if (!(rank>0)) throw std::runtime_error("Low rank must be larger then 0!");
// allocate memory and store compressed M2L operators
if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
// store U
U = new FReal [nnodes * rank];
FBlas::copy(rank*nnodes, _U, U);
......@@ -262,12 +276,21 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompress(const FReal epsilo
if (abs(i)>1 || abs(j)>1 || abs(k)>1) {
FBlas::copy(rank*rank, _C + counter*rank*rank, C + idx*rank*rank);
counter++;
} else FBlas::scal(rank*rank, FReal(0.), C + idx*rank*rank);
} else FBlas::setzero(rank*rank, C + idx*rank*rank);
}
if (counter != ninteractions)
throw std::runtime_error("Number of interactions must correspond to 316");
delete [] _C;
//////////////////////////////////////////////////////////
for (unsigned int n=0; n<nnodes; ++n) {
FBlas::scal(rank, FReal(1.) / weights[n], U+n, nnodes); // scale rows
FBlas::scal(rank, FReal(1.) / weights[n], B+n, nnodes); // scale rows
}
//////////////////////////////////////////////////////////
// return low rank
return rank;
}
......@@ -312,8 +335,7 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ComputeAndCompressAndStoreInBinaryFil
if (B != NULL) delete [] B;
if (C != NULL) delete [] C;
// write info
std::cout << "Compressed M2L operators (rank="<< rank
<< ") stored in binary file " << filename
std::cout << "Compressed M2L operators ("<< rank << ") stored in binary file " << filename
<< " in " << time.tacAndElapsed() << "sec." << std::endl;
}
......@@ -356,9 +378,8 @@ FChebM2LHandler<ORDER, MatrixKernelClass>::ReadFromBinaryFileAndSet()
} else throw std::runtime_error("File could not be opened to read");
stream.close();
// write info
std::cout << "Compressed M2L operators read from binary file "
<< filename << " in " << time.tacAndElapsed() << "sec."
<< std::endl;
std::cout << "Compressed M2L operators (" << rank << ") read from binary file "
<< filename << " in " << time.tacAndElapsed() << "sec." << std::endl;
}
/*
......@@ -532,7 +553,6 @@ unsigned int Compress(const FReal epsilon, const unsigned int ninteractions,
#endif // FCHEBM2LHANDLER_HPP
// [--END--]
......@@ -76,8 +76,14 @@ public:
// weights in 3d (tensor structure)
unsigned int node_ids[nnodes][3];
setNodeIds(node_ids);
for (unsigned int n=0; n<nnodes; ++n)
weights[n] = sqrt(weights_1d[node_ids[n][0]]*weights_1d[node_ids[n][1]]*weights_1d[node_ids[n][2]]);
for (unsigned int n=0; n<nnodes; ++n) {
FReal weight = FReal(1.); // no weighting
weight *= sqrt(weights_1d[node_ids[n][0]]*weights_1d[node_ids[n][1]]*weights_1d[node_ids[n][2]]); // 1/2
//weight *= weight; // 1
//weight *= weight; // 2
//weight *= weight; // 4
weights[n] = weight;
}
}
......
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