8 #ifndef MUELU_LOCALPERMUTATIONSTRATEGY_DEF_HPP_ 9 #define MUELU_LOCALPERMUTATIONSTRATEGY_DEF_HPP_ 13 #include <Xpetra_MultiVector.hpp> 14 #include <Xpetra_Matrix.hpp> 15 #include <Xpetra_MatrixMatrix.hpp> 16 #include <Xpetra_CrsGraph.hpp> 17 #include <Xpetra_Vector.hpp> 18 #include <Xpetra_VectorFactory.hpp> 19 #include <Xpetra_CrsMatrixWrap.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>
54 SC SC_ZERO = Teuchos::ScalarTraits<SC>::zero();
56 size_t nDofsPerNode = 1;
57 if (A->IsView(
"stridedMaps")) {
58 Teuchos::RCP<const Map> permRowMapStrided = A->getRowMap(
"stridedMaps");
59 nDofsPerNode = Teuchos::rcp_dynamic_cast<
const StridedMap>(permRowMapStrided)->getFixedBlockSize();
62 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
65 std::vector<std::pair<GlobalOrdinal,GlobalOrdinal> > RowColPairs;
68 if(permWidth_ != nDofsPerNode)
69 BuildPermutations(nDofsPerNode);
72 LocalOrdinal lonDofsPerNode = Teuchos::as<LocalOrdinal> (nDofsPerNode);
73 Teuchos::ArrayView<const LocalOrdinal> indices;
74 Teuchos::ArrayView<const Scalar> vals;
75 Teuchos::SerialDenseMatrix<LocalOrdinal,Scalar> subBlockMatrix(nDofsPerNode, nDofsPerNode,
true);
76 std::vector<GlobalOrdinal> growIds(nDofsPerNode);
80 LocalOrdinal numLocalNodes = A->getRowMap()->getNodeNumElements()/nDofsPerNode;
81 for (
LocalOrdinal node = 0; node < numLocalNodes; ++node) {
84 subBlockMatrix.putScalar();
89 for (
LocalOrdinal lrdof = 0; lrdof < lonDofsPerNode; ++lrdof) {
91 growIds[lrdof] = grow;
94 A->getLocalRowView(A->getRowMap()->getLocalElement(grow), indices, vals);
97 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
99 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
100 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
101 maxVal = Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
107 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
108 GlobalOrdinal gcol = A->getColMap()->getGlobalElement(indices[j]);
110 if (grnodeid == gcnodeid) {
111 if(maxVal != Teuchos::ScalarTraits<MT>::zero ()) {
112 subBlockMatrix(lrdof, gcol % nDofsPerNode) = vals[j]/maxVal;
115 subBlockMatrix(lrdof, gcol % nDofsPerNode) = vals[j];
116 std::cout <<
"maxVal never should be zero!!!!" << std::endl;
135 std::vector<Scalar> performance_vector = std::vector<Scalar>(result_permvecs_.size());
136 for (
size_t t = 0; t < result_permvecs_.size(); t++) {
137 std::vector<int> vv = result_permvecs_[t];
139 for(
size_t j=0; j<vv.size(); j++) {
140 value = value * subBlockMatrix(j,vv[j]);
142 performance_vector[t] = value;
155 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
156 MT maxVal = Teuchos::ScalarTraits<MT>::zero();
157 size_t maxPerformancePermutationIdx = 0;
158 for (
size_t j = 0; j < Teuchos::as<size_t>(performance_vector.size()); j++) {
159 if(Teuchos::ScalarTraits<Scalar>::magnitude(performance_vector[j]) > maxVal) {
160 maxVal = Teuchos::ScalarTraits<Scalar>::magnitude(performance_vector[j]);
161 maxPerformancePermutationIdx = j;
166 std::vector<int> bestPerformancePermutation = result_permvecs_[maxPerformancePermutationIdx];
167 for(
size_t t = 0; t<nDofsPerNode; t++) {
168 RowColPairs.push_back(std::make_pair(growIds[t],growIds[bestPerformancePermutation[t]]));
175 Teuchos::RCP<Vector> Pperm = VectorFactory::Build(A->getRowMap());
176 Teuchos::RCP<Vector> Qperm = VectorFactory::Build(A->getDomainMap());
178 Pperm->putScalar(SC_ZERO);
179 Qperm->putScalar(SC_ZERO);
181 Teuchos::ArrayRCP<Scalar> PpermData = Pperm->getDataNonConst(0);
182 Teuchos::ArrayRCP<Scalar> QpermData = Qperm->getDataNonConst(0);
184 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = RowColPairs.begin();
185 while(p != RowColPairs.end() ) {
190 LocalOrdinal ljk = A->getDomainMap()->getLocalElement(jk);
192 Pperm->replaceLocalValue(lik,ik);
193 Qperm->replaceLocalValue(ljk,ik);
195 p = RowColPairs.erase(p);
198 if(RowColPairs.size()>0) GetOStream(
Warnings0) <<
"MueLu::LocalPermutationStrategy: There are Row/col pairs left!" << std::endl;
204 Teuchos::RCP<CrsMatrixWrap> permPTmatrix = Teuchos::rcp(
new CrsMatrixWrap(A->getRowMap(),1));
205 Teuchos::RCP<CrsMatrixWrap> permQTmatrix = Teuchos::rcp(
new CrsMatrixWrap(A->getRowMap(),1));
207 for(
size_t row=0; row<A->getNodeNumRows(); row++) {
208 Teuchos::ArrayRCP<GlobalOrdinal> indoutP(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(PpermData[row])));
209 Teuchos::ArrayRCP<GlobalOrdinal> indoutQ(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(QpermData[row])));
210 Teuchos::ArrayRCP<Scalar> valout(1,Teuchos::ScalarTraits<Scalar>::one());
211 permPTmatrix->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indoutP.view(0,indoutP.size()), valout.view(0,valout.size()));
212 permQTmatrix->insertGlobalValues (A->getRowMap()->getGlobalElement(row), indoutQ.view(0,indoutQ.size()), valout.view(0,valout.size()));
215 permPTmatrix->fillComplete();
216 permQTmatrix->fillComplete();
230 Teuchos::RCP<Matrix> ApermQt = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *permQTmatrix,
false, GetOStream(
Statistics2),
true,
true);
231 Teuchos::RCP<Matrix> permPApermQt = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*permPmatrix,
false, *ApermQt,
false, GetOStream(
Statistics2),
true,
true);
240 Teuchos::RCP<Vector> diagVec = VectorFactory::Build(permPApermQt->getRowMap(),
true);
241 Teuchos::RCP<Vector> invDiagVec = VectorFactory::Build(permPApermQt->getRowMap(),
true);
242 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
243 Teuchos::ArrayRCP< Scalar > invDiagVecData = invDiagVec->getDataNonConst(0);
245 LO lCntZeroDiagonals = 0;
246 permPApermQt->getLocalDiagCopy(*diagVec);
247 for(
size_t i = 0; i<diagVec->getMap()->getNodeNumElements(); ++i) {
248 if(diagVecData[i] != SC_ZERO)
249 invDiagVecData[i] = Teuchos::ScalarTraits<Scalar>::one()/diagVecData[i];
251 invDiagVecData[i] = Teuchos::ScalarTraits<Scalar>::one();
258 GO gCntZeroDiagonals = 0;
259 GO glCntZeroDiagonals = Teuchos::as<GlobalOrdinal>(lCntZeroDiagonals);
260 MueLu_sumAll(comm,glCntZeroDiagonals,gCntZeroDiagonals);
261 GetOStream(
Statistics0) <<
"MueLu::LocalPermutationStrategy: found " << gCntZeroDiagonals <<
" zeros on diagonal" << std::endl;
264 Teuchos::RCP<CrsMatrixWrap> diagScalingOp = Teuchos::rcp(
new CrsMatrixWrap(permPApermQt->getRowMap(),1));
266 for(
size_t row=0; row<A->getNodeNumRows(); row++) {
267 Teuchos::ArrayRCP<GlobalOrdinal> indout(1,permPApermQt->getRowMap()->getGlobalElement(row));
268 Teuchos::ArrayRCP<Scalar> valout(1,invDiagVecData[row]);
269 diagScalingOp->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
271 diagScalingOp->fillComplete();
273 Teuchos::RCP<Matrix> scaledA = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*diagScalingOp,
false, *permPApermQt,
false, GetOStream(
Statistics2),
true,
true);
274 currentLevel.
Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(scaledA), genFactory);
276 currentLevel.
Set(
"permA", Teuchos::rcp_dynamic_cast<Matrix>(permPApermQt), genFactory);
277 currentLevel.
Set(
"permP", Teuchos::rcp_dynamic_cast<Matrix>(permPmatrix), genFactory);
278 currentLevel.
Set(
"permQT", Teuchos::rcp_dynamic_cast<Matrix>(permQTmatrix), genFactory);
279 currentLevel.
Set(
"permScaling", Teuchos::rcp_dynamic_cast<Matrix>(diagScalingOp), genFactory);
283 Teuchos::RCP<Vector> diagPVec = VectorFactory::Build(permPmatrix->getRowMap(),
true);
284 permPmatrix->getLocalDiagCopy(*diagPVec);
285 Teuchos::ArrayRCP< const Scalar > diagPVecData = diagPVec->getData(0);
288 for(
size_t i = 0; i<diagPVec->getMap()->getNodeNumElements(); ++i) {
289 if(diagPVecData[i] == SC_ZERO) {
290 lNumRowPermutations++;
295 MueLu_sumAll(diagPVec->getMap()->getComm(), lNumRowPermutations, gNumRowPermutations);
299 Teuchos::RCP<Vector> diagQTVec = VectorFactory::Build(permQTmatrix->getRowMap(),
true);
300 permQTmatrix->getLocalDiagCopy(*diagQTVec);
301 Teuchos::ArrayRCP< const Scalar > diagQTVecData = diagQTVec->getData(0);
304 for(
size_t i = 0; i<diagQTVec->getMap()->getNodeNumElements(); ++i) {
305 if(diagQTVecData[i] == SC_ZERO) {
306 lNumColPermutations++;
311 MueLu_sumAll(diagQTVec->getMap()->getComm(), lNumColPermutations, gNumColPermutations);
313 currentLevel.
Set(
"#RowPermutations", gNumRowPermutations, genFactory);
314 currentLevel.
Set(
"#ColPermutations", gNumColPermutations, genFactory);
315 currentLevel.
Set(
"#WideRangeRowPermutations", 0, genFactory);
316 currentLevel.
Set(
"#WideRangeColPermutations", 0, genFactory);
318 GetOStream(
Statistics0) <<
"#Row permutations/max possible permutations: " << gNumRowPermutations <<
"/" << diagPVec->getMap()->getGlobalNumElements() << std::endl;
319 GetOStream(
Statistics0) <<
"#Column permutations/max possible permutations: " << gNumColPermutations <<
"/" << diagQTVec->getMap()->getGlobalNumElements() << std::endl;
322 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
324 size_t nDofsPerNode = 1;
325 if (A->IsView(
"stridedMaps")) {
326 Teuchos::RCP<const Map> permRowMapStrided = A->getRowMap(
"stridedMaps");
327 nDofsPerNode = Teuchos::rcp_dynamic_cast<
const StridedMap>(permRowMapStrided)->getFixedBlockSize();
330 LocalOrdinal localDofId = localNodeId * nDofsPerNode + localDof;
332 return A->getRowMap()->getGlobalElement(localDofId);
335 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
337 size_t nDofsPerNode = 1;
338 if (A->IsView(
"stridedMaps")) {
339 Teuchos::RCP<const Map> permRowMapStrided = A->getRowMap(
"stridedMaps");
340 nDofsPerNode = Teuchos::rcp_dynamic_cast<
const StridedMap>(permRowMapStrided)->getFixedBlockSize();
Important warning messages (one line)
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
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
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
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(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
GlobalOrdinal globalDofId2globalNodeId(const Teuchos::RCP< Matrix > &A, GlobalOrdinal grid) const
void BuildPermutation(const Teuchos::RCP< Matrix > &A, const Teuchos::RCP< const Map > permRowMap, Level ¤tLevel, const FactoryBase *genFactory) const
build permutation operators