8 #ifndef MUELU_LOCALPERMUTATIONSTRATEGY_DEF_HPP_ 9 #define MUELU_LOCALPERMUTATIONSTRATEGY_DEF_HPP_ 21 #include "MueLu_Utilities.hpp" 26 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
29 permWidth_ = nDofsPerNode;
31 result_permvecs_.clear();
35 for(
size_t t = 0; t<nDofsPerNode; t++)
37 std::string cs = ss.str();
42 std::vector<int> newPerm(cs.length(),-1);
43 for(
size_t len=0; len<cs.length(); len++) {
44 newPerm[len] = Teuchos::as<int>(cs[len]-
'0');
46 result_permvecs_.push_back(newPerm);
48 }
while (std::next_permutation(cs.begin(),cs.end()));
51 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
53 size_t nDofsPerNode = 1;
54 if (A->IsView(
"stridedMaps")) {
56 nDofsPerNode = Teuchos::rcp_dynamic_cast<
const StridedMap>(permRowMapStrided)->getFixedBlockSize();
62 std::vector<std::pair<GlobalOrdinal,GlobalOrdinal> > RowColPairs;
65 if(permWidth_ != nDofsPerNode)
66 BuildPermutations(nDofsPerNode);
69 LocalOrdinal lonDofsPerNode = Teuchos::as<LocalOrdinal> (nDofsPerNode);
73 std::vector<GlobalOrdinal> growIds(nDofsPerNode);
77 LocalOrdinal numLocalNodes = A->getRowMap()->getNodeNumElements()/nDofsPerNode;
78 for (LocalOrdinal node = 0; node < numLocalNodes; ++node) {
86 for (LocalOrdinal lrdof = 0; lrdof < lonDofsPerNode; ++lrdof) {
87 GlobalOrdinal grow = getGlobalDofId(A, node, lrdof);
88 growIds[lrdof] = grow;
91 A->getLocalRowView(A->getRowMap()->getLocalElement(grow), indices, vals);
96 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
102 GlobalOrdinal grnodeid = globalDofId2globalNodeId(A,grow);
104 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
105 GlobalOrdinal gcol = A->getColMap()->getGlobalElement(indices[j]);
106 GlobalOrdinal gcnodeid = globalDofId2globalNodeId(A,gcol);
107 if (grnodeid == gcnodeid) {
109 subBlockMatrix(lrdof, gcol % nDofsPerNode) = vals[j]/maxVal;
112 subBlockMatrix(lrdof, gcol % nDofsPerNode) = vals[j];
113 std::cout <<
"maxVal never should be zero!!!!" << std::endl;
132 std::vector<Scalar> performance_vector = std::vector<Scalar>(result_permvecs_.size());
133 for (
size_t t = 0; t < result_permvecs_.size(); t++) {
134 std::vector<int> vv = result_permvecs_[t];
136 for(
size_t j=0; j<vv.size(); j++) {
137 value = value * subBlockMatrix(j,vv[j]);
139 performance_vector[t] = value;
154 size_t maxPerformancePermutationIdx = 0;
155 for (
size_t j = 0; j < Teuchos::as<size_t>(performance_vector.size()); j++) {
158 maxPerformancePermutationIdx = j;
163 std::vector<int> bestPerformancePermutation = result_permvecs_[maxPerformancePermutationIdx];
164 for(
size_t t = 0; t<nDofsPerNode; t++) {
165 RowColPairs.push_back(std::make_pair(growIds[t],growIds[bestPerformancePermutation[t]]));
175 Pperm->putScalar(0.0);
176 Qperm->putScalar(0.0);
181 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = RowColPairs.
begin();
182 while(p != RowColPairs.end() ) {
183 GlobalOrdinal ik = (*p).first;
184 GlobalOrdinal jk = (*p).second;
186 LocalOrdinal lik = A->getRowMap()->getLocalElement(ik);
187 LocalOrdinal ljk = A->getDomainMap()->getLocalElement(jk);
189 Pperm->replaceLocalValue(lik,ik);
190 Qperm->replaceLocalValue(ljk,ik);
192 p = RowColPairs.erase(p);
195 if(RowColPairs.size()>0) GetOStream(
Warnings0) <<
"MueLu::LocalPermutationStrategy: There are Row/col pairs left!" << std::endl;
204 for(
size_t row=0; row<A->getNodeNumRows(); row++) {
208 permPTmatrix->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indoutP.
view(0,indoutP.
size()), valout.
view(0,valout.
size()));
209 permQTmatrix->insertGlobalValues (A->getRowMap()->getGlobalElement(row), indoutQ.
view(0,indoutQ.
size()), valout.
view(0,valout.
size()));
212 permPTmatrix->fillComplete();
213 permQTmatrix->fillComplete();
242 LO lCntZeroDiagonals = 0;
243 permPApermQt->getLocalDiagCopy(*diagVec);
244 for(
size_t i = 0; i<diagVec->getMap()->getNodeNumElements(); ++i) {
245 if(diagVecData[i] != 0.0)
255 GO gCntZeroDiagonals = 0;
256 GO glCntZeroDiagonals = Teuchos::as<GlobalOrdinal>(lCntZeroDiagonals);
257 MueLu_sumAll(comm,glCntZeroDiagonals,gCntZeroDiagonals);
258 GetOStream(
Statistics0) <<
"MueLu::LocalPermutationStrategy: found " << gCntZeroDiagonals <<
" zeros on diagonal" << std::endl;
263 for(
size_t row=0; row<A->getNodeNumRows(); row++) {
266 diagScalingOp->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.
view(0,valout.
size()));
268 diagScalingOp->fillComplete();
271 currentLevel.
Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(scaledA), genFactory);
273 currentLevel.
Set(
"permA", Teuchos::rcp_dynamic_cast<Matrix>(permPApermQt), genFactory);
274 currentLevel.
Set(
"permP", Teuchos::rcp_dynamic_cast<Matrix>(permPmatrix), genFactory);
275 currentLevel.
Set(
"permQT", Teuchos::rcp_dynamic_cast<Matrix>(permQTmatrix), genFactory);
276 currentLevel.
Set(
"permScaling", Teuchos::rcp_dynamic_cast<Matrix>(diagScalingOp), genFactory);
281 permPmatrix->getLocalDiagCopy(*diagPVec);
283 GlobalOrdinal lNumRowPermutations = 0;
284 GlobalOrdinal gNumRowPermutations = 0;
285 for(
size_t i = 0; i<diagPVec->getMap()->getNodeNumElements(); ++i) {
286 if(diagPVecData[i] == 0.0) {
287 lNumRowPermutations++;
292 MueLu_sumAll(diagPVec->getMap()->getComm(), lNumRowPermutations, gNumRowPermutations);
297 permQTmatrix->getLocalDiagCopy(*diagQTVec);
299 GlobalOrdinal lNumColPermutations = 0;
300 GlobalOrdinal gNumColPermutations = 0;
301 for(
size_t i = 0; i<diagQTVec->getMap()->getNodeNumElements(); ++i) {
302 if(diagQTVecData[i] == 0.0) {
303 lNumColPermutations++;
308 MueLu_sumAll(diagQTVec->getMap()->getComm(), lNumColPermutations, gNumColPermutations);
310 currentLevel.
Set(
"#RowPermutations", gNumRowPermutations, genFactory);
311 currentLevel.
Set(
"#ColPermutations", gNumColPermutations, genFactory);
312 currentLevel.
Set(
"#WideRangeRowPermutations", 0, genFactory);
313 currentLevel.
Set(
"#WideRangeColPermutations", 0, genFactory);
315 GetOStream(
Statistics0) <<
"#Row permutations/max possible permutations: " << gNumRowPermutations <<
"/" << diagPVec->getMap()->getGlobalNumElements() << std::endl;
316 GetOStream(
Statistics0) <<
"#Column permutations/max possible permutations: " << gNumColPermutations <<
"/" << diagQTVec->getMap()->getGlobalNumElements() << std::endl;
319 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
321 size_t nDofsPerNode = 1;
322 if (A->IsView(
"stridedMaps")) {
324 nDofsPerNode = Teuchos::rcp_dynamic_cast<
const StridedMap>(permRowMapStrided)->getFixedBlockSize();
327 LocalOrdinal localDofId = localNodeId * nDofsPerNode + localDof;
329 return A->getRowMap()->getGlobalElement(localDofId);
332 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
334 size_t nDofsPerNode = 1;
335 if (A->IsView(
"stridedMaps")) {
337 nDofsPerNode = Teuchos::rcp_dynamic_cast<
const StridedMap>(permRowMapStrided)->getFixedBlockSize();
340 return (GlobalOrdinal) grid / (GlobalOrdinal)nDofsPerNode;
Important warning messages (one line)
#define MueLu_sumAll(rcpComm, in, out)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string())
Namespace for MueLu classes and methods.
GlobalOrdinal getGlobalDofId(const Teuchos::RCP< Matrix > &A, LocalOrdinal localNodeId, LocalOrdinal localDof) const
Print even more statistics.
Base class for factories (e.g., R, P, and A_coarse).
Print statistics that do not involve significant additional computation.
void BuildPermutations(size_t nDofsPerNode) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static magnitudeType magnitude(T a)
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
GlobalOrdinal globalDofId2globalNodeId(const Teuchos::RCP< Matrix > &A, GlobalOrdinal grid) const
ArrayView< T > view(size_type lowerOffset, size_type size) const
void BuildPermutation(const Teuchos::RCP< Matrix > &A, const Teuchos::RCP< const Map > permRowMap, Level ¤tLevel, const FactoryBase *genFactory) const
build permutation operators