MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_LocalPermutationStrategy_def.hpp
Go to the documentation of this file.
1/*
2 * MueLu_LocalPermutationStrategy_def.hpp
3 *
4 * Created on: Feb 19, 2013
5 * Author: tobias
6 */
7
8#ifndef MUELU_LOCALPERMUTATIONSTRATEGY_DEF_HPP_
9#define MUELU_LOCALPERMUTATIONSTRATEGY_DEF_HPP_
10
11#include <algorithm>
12
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>
20
21#include "MueLu_Utilities.hpp"
23
24namespace MueLu {
25
26 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
28
29 permWidth_ = nDofsPerNode;
30
31 result_permvecs_.clear();
32
33 // build permutation string
34 std::stringstream ss;
35 for(size_t t = 0; t<nDofsPerNode; t++)
36 ss << t;
37 std::string cs = ss.str();
38 //std::vector<std::string> result_perms;
39 do {
40 //result_perms.push_back(cs);
41
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');
45 }
46 result_permvecs_.push_back(newPerm);
47
48 } while (std::next_permutation(cs.begin(),cs.end()));
49 }
50
51 template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
52 void LocalPermutationStrategy<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildPermutation(const Teuchos::RCP<Matrix> & A, const Teuchos::RCP<const Map> /* permRowMap */, Level & currentLevel, const FactoryBase* genFactory) const {
53
54 SC SC_ZERO = Teuchos::ScalarTraits<SC>::zero();
55
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();
60 }
61
62 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
63
65 std::vector<std::pair<GlobalOrdinal,GlobalOrdinal> > RowColPairs;
66
67 // check whether we have to (re)build the permutation vector
68 if(permWidth_ != nDofsPerNode)
69 BuildPermutations(nDofsPerNode);
70
71 // declare local variables used inside loop
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);
77
78 // loop over local nodes
79 // TODO what about nOffset?
80 LocalOrdinal numLocalNodes = A->getRowMap()->getLocalNumElements()/nDofsPerNode;
81 for (LocalOrdinal node = 0; node < numLocalNodes; ++node) {
82
83 // zero out block matrix
84 subBlockMatrix.putScalar();
85
86 // loop over all DOFs in current node
87 // Note: were assuming constant number of Dofs per node here!
88 // TODO This is more complicated for variable dofs per node
89 for (LocalOrdinal lrdof = 0; lrdof < lonDofsPerNode; ++lrdof) {
90 GlobalOrdinal grow = getGlobalDofId(A, node, lrdof);
91 growIds[lrdof] = grow;
92
93 // extract local row information from matrix
94 A->getLocalRowView(A->getRowMap()->getLocalElement(grow), indices, vals);
95
96 // find column entry with max absolute value
97 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
98 MT maxVal = 0.0;
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]);
102 }
103 }
104
105 GlobalOrdinal grnodeid = globalDofId2globalNodeId(A,grow);
106
107 for (size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
108 GlobalOrdinal gcol = A->getColMap()->getGlobalElement(indices[j]);
109 GlobalOrdinal gcnodeid = globalDofId2globalNodeId(A,gcol); // -> global node id
110 if (grnodeid == gcnodeid) {
111 if(maxVal != Teuchos::ScalarTraits<MT>::zero ()) {
112 subBlockMatrix(lrdof, gcol % nDofsPerNode) = vals[j]/maxVal;
113 } else
114 {
115 subBlockMatrix(lrdof, gcol % nDofsPerNode) = vals[j]; // there is a problem
116 std::cout << "maxVal never should be zero!!!!" << std::endl;
117 }
118 }
119 }
120 }
121
122 // now we have the sub block matrix
123
124 // build permutation string
125 /*std::stringstream ss;
126 for(size_t t = 0; t<nDofsPerNode; t++)
127 ss << t;
128 std::string cs = ss.str();
129 std::vector<std::string> result_perms;
130 do {
131 result_perms.push_back(cs);
132 //std::cout << result_perms.back() << std::endl;
133 } while (std::next_permutation(cs.begin(),cs.end()));*/
134
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];
138 Scalar value = 1.0;
139 for(size_t j=0; j<vv.size(); j++) {
140 value = value * subBlockMatrix(j,vv[j]);
141 }
142 performance_vector[t] = value;
143 }
144 /*for(size_t t = 0; t < result_perms.size(); t++) {
145 std::string s = result_perms[t];
146 Scalar value = 1.0;
147 for(size_t len=0; len<s.length(); len++) {
148 int col = Teuchos::as<int>(s[len]-'0');
149 value = value * subBlockMatrix(len,col);
150 }
151 performance_vector[t] = value;
152 }*/
153
154 // find permutation with maximum performance 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;
162 }
163 }
164
165 // build RowColPairs for best permutation
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]]));
169 }
170
171 } // end loop over local nodes
172
173
174 // build Pperm and Qperm vectors
175 Teuchos::RCP<Vector> Pperm = VectorFactory::Build(A->getRowMap());
176 Teuchos::RCP<Vector> Qperm = VectorFactory::Build(A->getDomainMap());
177
178 Pperm->putScalar(SC_ZERO);
179 Qperm->putScalar(SC_ZERO);
180
181 Teuchos::ArrayRCP<Scalar> PpermData = Pperm->getDataNonConst(0);
182 Teuchos::ArrayRCP<Scalar> QpermData = Qperm->getDataNonConst(0);
183
184 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = RowColPairs.begin();
185 while(p != RowColPairs.end() ) {
186 GlobalOrdinal ik = (*p).first;
187 GlobalOrdinal jk = (*p).second;
188
189 LocalOrdinal lik = A->getRowMap()->getLocalElement(ik);
190 LocalOrdinal ljk = A->getDomainMap()->getLocalElement(jk);
191
192 Pperm->replaceLocalValue(lik,ik);
193 Qperm->replaceLocalValue(ljk,ik);
194
195 p = RowColPairs.erase(p);
196 }
197
198 if(RowColPairs.size()>0) GetOStream(Warnings0) << "MueLu::LocalPermutationStrategy: There are Row/col pairs left!" << std::endl;
199
200 // Qperm should be fine
201 // build matrices
202
203 // create new empty Matrix
204 Teuchos::RCP<CrsMatrixWrap> permPTmatrix = Teuchos::rcp(new CrsMatrixWrap(A->getRowMap(),1));
205 Teuchos::RCP<CrsMatrixWrap> permQTmatrix = Teuchos::rcp(new CrsMatrixWrap(A->getRowMap(),1));
206
207 for(size_t row=0; row<A->getLocalNumRows(); row++) {
208 Teuchos::ArrayRCP<GlobalOrdinal> indoutP(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(PpermData[row]))); // column idx for Perm^T
209 Teuchos::ArrayRCP<GlobalOrdinal> indoutQ(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(QpermData[row]))); // column idx for Qperm
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()));
213 }
214
215 permPTmatrix->fillComplete();
216 permQTmatrix->fillComplete();
217
218 Teuchos::RCP<Matrix> permPmatrix = Utilities::Transpose(*permPTmatrix,true);
219
220 /*for(size_t row=0; row<permPTmatrix->getLocalNumRows(); row++) {
221 if(permPTmatrix->getNumEntriesInLocalRow(row) != 1)
222 GetOStream(Warnings0) <<"#entries in row " << row << " of permPTmatrix is " << permPTmatrix->getNumEntriesInLocalRow(row) << std::endl;
223 if(permPmatrix->getNumEntriesInLocalRow(row) != 1)
224 GetOStream(Warnings0) <<"#entries in row " << row << " of permPmatrix is " << permPmatrix->getNumEntriesInLocalRow(row) << std::endl;
225 if(permQTmatrix->getNumEntriesInLocalRow(row) != 1)
226 GetOStream(Warnings0) <<"#entries in row " << row << " of permQmatrix is " << permQTmatrix->getNumEntriesInLocalRow(row) << std::endl;
227 }*/
228
229 // build permP * A * permQT
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);
232
233 /*
234 MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("A.mat", *A);
235 MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permP.mat", *permPmatrix);
236 MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permQt.mat", *permQTmatrix);
237 MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permPApermQt.mat", *permPApermQt);
238 */
239 // build scaling matrix
240 Teuchos::RCP<Vector> diagVec = VectorFactory::Build(permPApermQt->getRowMap(),true);
241 Teuchos::RCP<Vector> invDiagVec = VectorFactory::Build(permPApermQt->getRowMap(),true);
242 Teuchos::ArrayRCP< Scalar > invDiagVecData = invDiagVec->getDataNonConst(0);
243
244 LO lCntZeroDiagonals = 0;
245 permPApermQt->getLocalDiagCopy(*diagVec);
246 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
247 for(size_t i = 0; i<diagVec->getMap()->getLocalNumElements(); ++i) {
248 if(diagVecData[i] != SC_ZERO)
249 invDiagVecData[i] = Teuchos::ScalarTraits<Scalar>::one()/diagVecData[i];
250 else {
251 invDiagVecData[i] = Teuchos::ScalarTraits<Scalar>::one();
252 lCntZeroDiagonals++;
253 //GetOStream(Statistics0) << "MueLu::LocalPermutationStrategy: found zero on diagonal in row " << i << std::endl;
254 }
255 }
256
257#if 0
258 GO gCntZeroDiagonals = 0;
259 GO glCntZeroDiagonals = Teuchos::as<GlobalOrdinal>(lCntZeroDiagonals); /* LO->GO conversion */
260 MueLu_sumAll(comm,glCntZeroDiagonals,gCntZeroDiagonals);
261 GetOStream(Statistics0) << "MueLu::LocalPermutationStrategy: found " << gCntZeroDiagonals << " zeros on diagonal" << std::endl;
262#endif
263
264 Teuchos::RCP<CrsMatrixWrap> diagScalingOp = Teuchos::rcp(new CrsMatrixWrap(permPApermQt->getRowMap(),1));
265
266 for(size_t row=0; row<A->getLocalNumRows(); row++) {
267 Teuchos::ArrayRCP<GlobalOrdinal> indout(1,permPApermQt->getRowMap()->getGlobalElement(row)); // column idx for Perm^T
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()));
270 }
271 diagScalingOp->fillComplete();
272
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);
275
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);
280
282 // count zeros on diagonal in P -> number of row permutations
283 Teuchos::RCP<Vector> diagPVec = VectorFactory::Build(permPmatrix->getRowMap(),true);
284 permPmatrix->getLocalDiagCopy(*diagPVec);
285 Teuchos::ArrayRCP< const Scalar > diagPVecData = diagPVec->getData(0);
286 GlobalOrdinal lNumRowPermutations = 0;
287 GlobalOrdinal gNumRowPermutations = 0;
288 for(size_t i = 0; i<diagPVec->getMap()->getLocalNumElements(); ++i) {
289 if(diagPVecData[i] == SC_ZERO) {
290 lNumRowPermutations++;
291 }
292 }
293
294 // sum up all entries in multipleColRequests over all processors
295 MueLu_sumAll(diagPVec->getMap()->getComm(), lNumRowPermutations, gNumRowPermutations);
296
298 // count zeros on diagonal in Q^T -> number of column permutations
299 Teuchos::RCP<Vector> diagQTVec = VectorFactory::Build(permQTmatrix->getRowMap(),true);
300 permQTmatrix->getLocalDiagCopy(*diagQTVec);
301 Teuchos::ArrayRCP< const Scalar > diagQTVecData = diagQTVec->getData(0);
302 GlobalOrdinal lNumColPermutations = 0;
303 GlobalOrdinal gNumColPermutations = 0;
304 for(size_t i = 0; i<diagQTVec->getMap()->getLocalNumElements(); ++i) {
305 if(diagQTVecData[i] == SC_ZERO) {
306 lNumColPermutations++;
307 }
308 }
309
310 // sum up all entries in multipleColRequests over all processors
311 MueLu_sumAll(diagQTVec->getMap()->getComm(), lNumColPermutations, gNumColPermutations);
312
313 currentLevel.Set("#RowPermutations", gNumRowPermutations, genFactory/*this*/);
314 currentLevel.Set("#ColPermutations", gNumColPermutations, genFactory/*this*/);
315 currentLevel.Set("#WideRangeRowPermutations", 0, genFactory/*this*/);
316 currentLevel.Set("#WideRangeColPermutations", 0, genFactory/*this*/);
317
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;
320 }
321
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();
328 }
329
330 LocalOrdinal localDofId = localNodeId * nDofsPerNode + localDof;
331
332 return A->getRowMap()->getGlobalElement(localDofId);
333 }
334
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();
341 }
342
343 return (GlobalOrdinal) grid / (GlobalOrdinal)nDofsPerNode; // TODO what about nOffset???
344 }
345
346} // namespace MueLu
347
348
349#endif /* MUELU_LOCALPERMUTATIONSTRATEGY_DEF_HPP_ */
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Base class for factories (e.g., R, P, and A_coarse).
Class that holds all level-specific information.
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 &currentLevel, const FactoryBase *genFactory) const
build permutation operators
GlobalOrdinal getGlobalDofId(const Teuchos::RCP< Matrix > &A, LocalOrdinal localNodeId, LocalOrdinal localDof) const
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 > &params=Teuchos::null)
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Statistics2
Print even more statistics.
@ Statistics0
Print statistics that do not involve significant additional computation.