diff --git a/Src/Chebyshev/FChebKernels.hpp b/Src/Chebyshev/FChebKernels.hpp index 13641ae7e5b867fc1bb4dfed74bb46b996b6a701..aa65a221f77f68a3cdc8949f534f2765695a9eb3 100644 --- a/Src/Chebyshev/FChebKernels.hpp +++ b/Src/Chebyshev/FChebKernels.hpp @@ -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) diff --git a/Src/Chebyshev/FChebM2LHandler.hpp b/Src/Chebyshev/FChebM2LHandler.hpp index abfe6a222c80fd659686b0de6b808c7796eeaa82..7f17a4147717bda76928c5efa06117441f130a8a 100644 --- a/Src/Chebyshev/FChebM2LHandler.hpp +++ b/Src/Chebyshev/FChebM2LHandler.hpp @@ -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) @@ -237,13 +241,23 @@ 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--] diff --git a/Src/Chebyshev/FChebTensor.hpp b/Src/Chebyshev/FChebTensor.hpp index 89fb7ed8329b9adf8d91fa3003e0b319dbae6eb8..274ff2ab1fa6b2b6ffd4f70d2f00e37b56e76625 100644 --- a/Src/Chebyshev/FChebTensor.hpp +++ b/Src/Chebyshev/FChebTensor.hpp @@ -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; + } }