MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
24 namespace 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>
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);
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()->getNodeNumElements()/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
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) {
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
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->getNodeNumRows(); 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
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->getNodeNumRows(); 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< const Scalar > diagVecData = diagVec->getData(0);
243  Teuchos::ArrayRCP< Scalar > invDiagVecData = invDiagVec->getDataNonConst(0);
244 
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];
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->getNodeNumRows(); 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()->getNodeNumElements(); ++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()->getNodeNumElements(); ++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_ */
Important warning messages (one line)
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
void BuildPermutation(const Teuchos::RCP< Matrix > &A, const Teuchos::RCP< const Map > permRowMap, Level &currentLevel, const FactoryBase *genFactory) const
build permutation operators
size_type size() const
Print even more statistics.
GlobalOrdinal getGlobalDofId(const Teuchos::RCP< Matrix > &A, LocalOrdinal localNodeId, LocalOrdinal localDof) const
Base class for factories (e.g., R, P, and A_coarse).
Print statistics that do not involve significant additional computation.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void BuildPermutations(size_t nDofsPerNode) const
int putScalar(const ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
GlobalOrdinal globalDofId2globalNodeId(const Teuchos::RCP< Matrix > &A, GlobalOrdinal grid) 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)
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static magnitudeType magnitude(T a)
iterator begin() const
ArrayView< T > view(size_type lowerOffset, size_type size) const