MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_AlgebraicPermutationStrategy_def.hpp
Go to the documentation of this file.
1 /*
2  * MueLu_AlgebraicPermutationStrategy_def.hpp
3  *
4  * Created on: Feb 20, 2013
5  * Author: tobias
6  */
7 
8 #ifndef MUELU_ALGEBRAICPERMUTATIONSTRATEGY_DEF_HPP_
9 #define MUELU_ALGEBRAICPERMUTATIONSTRATEGY_DEF_HPP_
10 
11 #include <queue>
12 
13 #include <Teuchos_ScalarTraits.hpp>
14 
15 #include <Xpetra_MultiVector.hpp>
16 #include <Xpetra_Matrix.hpp>
17 #include <Xpetra_CrsGraph.hpp>
18 #include <Xpetra_Vector.hpp>
20 #include <Xpetra_VectorFactory.hpp>
21 #include <Xpetra_CrsMatrixWrap.hpp>
22 #include <Xpetra_Export.hpp>
23 #include <Xpetra_ExportFactory.hpp>
24 #include <Xpetra_Import.hpp>
25 #include <Xpetra_ImportFactory.hpp>
26 #include <Xpetra_MatrixMatrix.hpp>
27 
28 #include "MueLu_Utilities.hpp"
30 
31 namespace MueLu {
32 
33  template<class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
36  Level & currentLevel, const FactoryBase* genFactory) const {
37 
38  const Scalar SC_ZERO = Teuchos::ScalarTraits<Scalar>::zero();
39  const Scalar SC_ONE = Teuchos::ScalarTraits<Scalar>::one();
40 
42  const MT MT_ZERO = Teuchos::ScalarTraits<MT>::zero();
43  const MT MT_ONE = Teuchos::ScalarTraits<MT>::one();
44 
45  const Teuchos::RCP< const Teuchos::Comm< int > > comm = A->getRowMap()->getComm();
46  int numProcs = comm->getSize();
47  int myRank = comm->getRank();
48 
49  size_t nDofsPerNode = 1;
50  if (A->IsView("stridedMaps")) {
51  Teuchos::RCP<const Map> permRowMapStrided = A->getRowMap("stridedMaps");
52  nDofsPerNode = Teuchos::rcp_dynamic_cast<const StridedMap>(permRowMapStrided)->getFixedBlockSize();
53  }
54 
55  std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > permutedDiagCandidates;
56  std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > keepDiagonalEntries;
57  std::vector<MT> Weights;
58 
59  // loop over all local rows in matrix A and keep diagonal entries if corresponding
60  // matrix rows are not contained in permRowMap
61  for (size_t row = 0; row < A->getRowMap()->getNodeNumElements(); row++) {
62  GlobalOrdinal grow = A->getRowMap()->getGlobalElement(row);
63 
64  if(permRowMap->isNodeGlobalElement(grow) == true) continue;
65 
66  size_t nnz = A->getNumEntriesInLocalRow(row);
67 
68  // extract local row information from matrix
71  A->getLocalRowView(row, indices, vals);
72 
73  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz, Exceptions::RuntimeError,
74  "MueLu::PermutationFactory::Build: number of nonzeros not equal to number of indices? Error.");
75 
76  // find column entry with max absolute value
77  GlobalOrdinal gMaxValIdx = 0;
78  MT norm1 = MT_ZERO;
79  MT maxVal = MT_ZERO;
80  for (size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
82  if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
84  gMaxValIdx = A->getColMap()->getGlobalElement(indices[j]);
85  }
86  }
87 
88  if(grow == gMaxValIdx) // only keep row/col pair if it's diagonal dominant!!!
89  keepDiagonalEntries.push_back(std::make_pair(grow,grow));
90  }
91 
93  // handle rows that are marked to be relevant for permutations
94  for (size_t row = 0; row < permRowMap->getNodeNumElements(); row++) {
95  GlobalOrdinal grow = permRowMap->getGlobalElement(row);
96  LocalOrdinal lArow = A->getRowMap()->getLocalElement(grow);
97  size_t nnz = A->getNumEntriesInLocalRow(lArow);
98 
99  // extract local row information from matrix
102  A->getLocalRowView(lArow, indices, vals);
103 
104  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz, Exceptions::RuntimeError,
105  "MueLu::PermutationFactory::Build: number of nonzeros not equal to number of indices? Error.");
106 
107  // find column entry with max absolute value
108  GlobalOrdinal gMaxValIdx = 0;
109  MT norm1 = MT_ZERO;
110  MT maxVal = MT_ZERO;
111  for (size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
113  if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
115  gMaxValIdx = A->getColMap()->getGlobalElement(indices[j]);
116  }
117  }
118 
119  if(maxVal > Teuchos::ScalarTraits<MT>::zero ()) { // keep only max Entries \neq 0.0
120  permutedDiagCandidates.push_back(std::make_pair(grow,gMaxValIdx));
121  Weights.push_back(maxVal/(norm1*Teuchos::as<MT>(nnz)));
122  } else {
123  std::cout << "ATTENTION: row " << grow << " has only zero entries -> singular matrix!" << std::endl;
124  }
125 
126  }
127 
128  // sort Weights in descending order
129  std::vector<int> permutation;
130  sortingPermutation(Weights,permutation);
131 
132  // create new vector with exactly one possible entry for each column
133 
134  // each processor which requests the global column id gcid adds 1 to gColVec
135  // gColVec will be summed up over all processors and communicated to gDomVec
136  // which is based on the non-overlapping domain map of A.
137 
138  Teuchos::RCP<Vector> gColVec = VectorFactory::Build(A->getColMap());
139  Teuchos::RCP<Vector> gDomVec = VectorFactory::Build(A->getDomainMap());
140  gColVec->putScalar(SC_ZERO);
141  gDomVec->putScalar(SC_ZERO);
142 
143  // put in all keep diagonal entries
144  for (typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::const_iterator p = keepDiagonalEntries.begin(); p != keepDiagonalEntries.end(); ++p) {
145  gColVec->sumIntoGlobalValue((*p).second,1.0);
146  }
147 
148  Teuchos::RCP<Export> exporter = ExportFactory::Build(gColVec->getMap(), gDomVec->getMap());
149  gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD); // communicate blocked gcolids to all procs
150  gColVec->doImport(*gDomVec,*exporter,Xpetra::INSERT);
151 
152  std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > permutedDiagCandidatesFiltered; // TODO reserve memory
153  std::map<GlobalOrdinal, MT> gColId2Weight;
154 
155  Teuchos::ArrayRCP< Scalar > ddata = gColVec->getDataNonConst(0);
156  for(size_t i = 0; i < permutedDiagCandidates.size(); ++i) {
157  // loop over all candidates
158  std::pair<GlobalOrdinal, GlobalOrdinal> pp = permutedDiagCandidates[permutation[i]];
159  GlobalOrdinal grow = pp.first;
160  GlobalOrdinal gcol = pp.second;
161 
162  LocalOrdinal lcol = A->getColMap()->getLocalElement(gcol);
163  if(Teuchos::ScalarTraits<Scalar>::real (ddata[lcol]) > MT_ZERO){
164  continue; // skip lcol: column already handled by another row
165  }
166 
167  // mark column as already taken
168  ddata[lcol] += SC_ONE;
169 
170  permutedDiagCandidatesFiltered.push_back(std::make_pair(grow,gcol));
171  gColId2Weight[gcol] = Weights[permutation[i]];
172  }
173 
174  // communicate how often each column index is requested by the different procs
175  gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD);
176  gColVec->doImport(*gDomVec,*exporter,Xpetra::INSERT); // probably not needed // TODO check me
177 
178  //*****************************************************************************************
179  // first communicate ALL global ids of column indices which are requested by more
180  // than one proc to all other procs
181  // detect which global col indices are requested by more than one proc
182  // and store them in the multipleColRequests vector
183  std::vector<GlobalOrdinal> multipleColRequests; // store all global column indices from current processor that are also
184  // requested by another processor. This is possible, since they are stored
185  // in gDomVec which is based on the nonoverlapping domain map. That is, each
186  // global col id is handled by exactly one proc.
187  std::queue<GlobalOrdinal> unusedColIdx; // unused column indices on current processor
188 
189  for(size_t sz = 0; sz < gDomVec->getLocalLength(); ++sz) {
190  Teuchos::ArrayRCP< const Scalar > arrDomVec = gDomVec->getData(0);
191  //
192  // FIXME (mfh 30 Oct 2015) I _think_ it's OK to check just the
193  // real part, because this is a count. (Shouldn't MueLu use
194  // mag_type instead of Scalar here to save space?)
195  //
196  if(Teuchos::ScalarTraits<Scalar>::real(arrDomVec[sz]) > MT_ONE) {
197  multipleColRequests.push_back(gDomVec->getMap()->getGlobalElement(sz));
198  } else if(Teuchos::ScalarTraits<Scalar>::real (arrDomVec[sz]) == MT_ZERO) {
199  unusedColIdx.push(gDomVec->getMap()->getGlobalElement(sz));
200  }
201  }
202 
203  // communicate the global number of column indices which are requested by more than one proc
204  LocalOrdinal localMultColRequests = Teuchos::as<LocalOrdinal>(multipleColRequests.size());
205  LocalOrdinal globalMultColRequests = 0;
206 
207  // sum up all entries in multipleColRequests over all processors
208  //
209  // FIXME (lbv 11 Feb 2019) why cast localMultColRequests as LocalOrdinal
210  // when it was properly declared as such just a line above?
211  //
212  MueLu_sumAll(gDomVec->getMap()->getComm(), (LocalOrdinal)localMultColRequests, globalMultColRequests);
213 
214  if(globalMultColRequests > 0) {
215  // special handling: two processors request the same global column id.
216  // decide which processor gets it
217 
218  // distribute number of multipleColRequests to all processors
219  // each processor stores how many column ids for exchange are handled by the cur proc
220  std::vector<GlobalOrdinal> numMyMultColRequests(numProcs, 0);
221  std::vector<GlobalOrdinal> numGlobalMultColRequests(numProcs, 0);
222  numMyMultColRequests[myRank] = localMultColRequests;
223  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
224  numMyMultColRequests.data(),
225  numGlobalMultColRequests.data());
226 
227  // communicate multipleColRequests entries to all processors
228  int nMyOffset = 0;
229  for (int i=0; i<myRank-1; i++)
230  nMyOffset += numGlobalMultColRequests[i]; // calculate offset to store the weights on the corresponding place in procOverlappingWeights
231 
232  const GlobalOrdinal zero = 0;
233  std::vector<GlobalOrdinal> procMultRequestedColIds(globalMultColRequests, zero);
234  std::vector<GlobalOrdinal> global_procMultRequestedColIds(globalMultColRequests, zero);
235 
236  // loop over all local column GIDs that are also requested by other procs
237  for(size_t i = 0; i < multipleColRequests.size(); i++) {
238  procMultRequestedColIds[nMyOffset + i] = multipleColRequests[i]; // all weights are > 0 ?
239  }
240 
241  // template ordinal, package (double)
242  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(globalMultColRequests),
243  procMultRequestedColIds.data(),
244  global_procMultRequestedColIds.data());
245 
246  // loop over global_procOverlappingWeights and eliminate wrong entries...
247  for (size_t k = 0; k<global_procMultRequestedColIds.size(); k++) {
248  GlobalOrdinal globColId = global_procMultRequestedColIds[k];
249 
250  std::vector<MT> MyWeightForColId(numProcs, MT_ZERO);
251  std::vector<MT> GlobalWeightForColId(numProcs, MT_ZERO);
252 
253  if(gColVec->getMap()->isNodeGlobalElement(globColId)) {
254  MyWeightForColId[myRank] = gColId2Weight[globColId];
255  } else {
256  MyWeightForColId[myRank] = MT_ZERO;
257  }
258 
259  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
260  MyWeightForColId.data(),
261  GlobalWeightForColId.data());
262 
263  if(gColVec->getMap()->isNodeGlobalElement(globColId)) {
264  // note: 2 procs could have the same weight for a column index.
265  // pick the first one.
266  MT winnerValue = MT_ZERO;
267  int winnerProcRank = 0;
268  for (int proc = 0; proc < numProcs; proc++) {
269  if(GlobalWeightForColId[proc] > winnerValue) {
270  winnerValue = GlobalWeightForColId[proc];
271  winnerProcRank = proc;
272  }
273  }
274 
275  // winnerProcRank is the winner for handling globColId.
276  // winnerProcRank is unique (even if two procs have the same weight for a column index)
277 
278  if(myRank != winnerProcRank) {
279  // remove corresponding entry from permutedDiagCandidatesFiltered
280  typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = permutedDiagCandidatesFiltered.begin();
281  while(p != permutedDiagCandidatesFiltered.end() )
282  {
283  if((*p).second == globColId)
284  p = permutedDiagCandidatesFiltered.erase(p);
285  else
286  p++;
287  }
288  }
289 
290  } // end if isNodeGlobalElement
291  } // end loop over global_procOverlappingWeights and eliminate wrong entries...
292  } // end if globalMultColRequests > 0
293 
294  // put together all pairs:
295  //size_t sizeRowColPairs = keepDiagonalEntries.size() + permutedDiagCandidatesFiltered.size();
296  std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > RowColPairs;
297  RowColPairs.insert( RowColPairs.end(), keepDiagonalEntries.begin(), keepDiagonalEntries.end());
298  RowColPairs.insert( RowColPairs.end(), permutedDiagCandidatesFiltered.begin(), permutedDiagCandidatesFiltered.end());
299 
300 #ifdef DEBUG_OUTPUT
301  //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
302  // plausibility check
303  gColVec->putScalar(SC_ZERO);
304  gDomVec->putScalar(SC_ZERO);
305  typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator pl = RowColPairs.begin();
306  while(pl != RowColPairs.end() )
307  {
308  //GlobalOrdinal ik = (*pl).first;
309  GlobalOrdinal jk = (*pl).second;
310 
311  gColVec->sumIntoGlobalValue(jk,1.0);
312  pl++;
313  }
314  gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD);
315  for(size_t sz = 0; sz < gDomVec->getLocalLength(); ++sz) {
316  Teuchos::ArrayRCP< const Scalar > arrDomVec = gDomVec->getData(0);
317  if(arrDomVec[sz] > 1.0) {
318  GetOStream(Runtime0) << "RowColPairs has multiple column [" << sz << "]=" << arrDomVec[sz] << std::endl;
319  } else if(arrDomVec[sz] == SC_ZERO) {
320  GetOStream(Runtime0) << "RowColPairs has empty column [" << sz << "]=" << arrDomVec[sz] << std::endl;
321  }
322  }
323  //&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
324 #endif
325 
327  // assumption: on each processor RowColPairs now contains
328  // a valid set of (row,column) pairs, where the row entries
329  // are a subset of the processor's rows and the column entries
330  // are unique throughout all processors.
331  // Note: the RowColPairs are only defined for a subset of all rows,
332  // so there might be rows without an entry in RowColPairs.
333  // It can be, that some rows seem to be missing in RowColPairs, since
334  // the entry in that row with maximum absolute value has been reserved
335  // by another row already (e.g. as already diagonal dominant row outside
336  // of perRowMap).
337  // In fact, the RowColPairs vector only defines the (row,column) pairs
338  // that will be definitely moved to the diagonal after permutation.
339 
340 #ifdef DEBUG_OUTPUT
341  // for (typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::const_iterator p = RowColPairs.begin(); p != RowColPairs.end(); ++p) {
342  // std::cout << "proc: " << myRank << " r/c: " << (*p).first << "/" << (*p).second << std::endl;
343  // }
344  // for (typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::const_iterator p = RowColPairs.begin(); p != RowColPairs.end(); ++p)
345  // {
347  // std::cout << (*p).first +1 << " " << (*p).second+1 << std::endl;
348  // }
349  // std::cout << "\n";
350 #endif
351 
352  // vectors to store permutation information
353  Teuchos::RCP<Vector> Pperm = VectorFactory::Build(A->getRowMap());
354  Teuchos::RCP<Vector> Qperm = VectorFactory::Build(A->getDomainMap()); // global variant (based on domain map)
355  Teuchos::RCP<Vector> lQperm = VectorFactory::Build(A->getColMap()); // local variant (based on column map)
356 
357  Teuchos::ArrayRCP< Scalar > PpermData = Pperm->getDataNonConst(0);
358  Teuchos::ArrayRCP< Scalar > QpermData = Qperm->getDataNonConst(0);
359 
360  Pperm->putScalar(SC_ZERO);
361  Qperm->putScalar(SC_ZERO);
362  lQperm->putScalar(SC_ZERO);
363 
364  // setup exporter for Qperm
365  Teuchos::RCP<Export> QpermExporter = ExportFactory::Build(lQperm->getMap(), Qperm->getMap());
366 
367  Teuchos::RCP<Vector> RowIdStatus = VectorFactory::Build(A->getRowMap());
368  Teuchos::RCP<Vector> ColIdStatus = VectorFactory::Build(A->getDomainMap()); // global variant (based on domain map)
369  Teuchos::RCP<Vector> lColIdStatus = VectorFactory::Build(A->getColMap()); // local variant (based on column map)
370  Teuchos::RCP<Vector> ColIdUsed = VectorFactory::Build(A->getDomainMap()); // mark column ids to be already in use
371  Teuchos::ArrayRCP< Scalar > RowIdStatusArray = RowIdStatus->getDataNonConst(0);
372  Teuchos::ArrayRCP< Scalar > ColIdStatusArray = ColIdStatus->getDataNonConst(0);
373  Teuchos::ArrayRCP< Scalar > lColIdStatusArray = lColIdStatus->getDataNonConst(0);
374  Teuchos::ArrayRCP< Scalar > ColIdUsedArray = ColIdUsed->getDataNonConst(0); // not sure about this
375  RowIdStatus->putScalar(SC_ZERO);
376  ColIdStatus->putScalar(SC_ZERO);
377  lColIdStatus->putScalar(SC_ZERO);
378  ColIdUsed->putScalar(SC_ZERO); // no column ids are used
379 
380  // count wide-range permutations
381  // a wide-range permutation is defined as a permutation of rows/columns which do not
382  // belong to the same node
383  LocalOrdinal lWideRangeRowPermutations = 0;
384  GlobalOrdinal gWideRangeRowPermutations = 0;
385  LocalOrdinal lWideRangeColPermutations = 0;
386  GlobalOrdinal gWideRangeColPermutations = 0;
387 
388  // run 1: mark all "identity" permutations
389  typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = RowColPairs.begin();
390  while(p != RowColPairs.end() )
391  {
392  GlobalOrdinal ik = (*p).first;
393  GlobalOrdinal jk = (*p).second;
394 
395  LocalOrdinal lik = A->getRowMap()->getLocalElement(ik);
396  LocalOrdinal ljk = A->getColMap()->getLocalElement(jk);
397 
398  if(RowIdStatusArray[lik] == SC_ZERO) {
399  RowIdStatusArray[lik] = SC_ONE; // use this row id
400  lColIdStatusArray[ljk] = SC_ONE; // use this column id
401  Pperm->replaceLocalValue(lik, ik);
402  lQperm->replaceLocalValue(ljk, ik); // use column map
403  ColIdUsed->replaceGlobalValue(ik,SC_ONE); // ik is now used
404  p = RowColPairs.erase(p);
405 
406  // detect wide range permutations
407  if(floor(ik/nDofsPerNode) != floor(jk/nDofsPerNode)) {
408  lWideRangeColPermutations++;
409  }
410  }
411  else
412  p++;
413  }
414 
415  // communicate column map -> domain map
416  Qperm->doExport(*lQperm, *QpermExporter, Xpetra::ABSMAX);
417  ColIdStatus->doExport(*lColIdStatus, *QpermExporter, Xpetra::ABSMAX);
418 
419  // plausibility check
420  if(RowColPairs.size()>0) GetOStream(Warnings0) << "MueLu::PermutationFactory: There are Row/Col pairs left!!!" << std::endl; // TODO fix me
421 
422  // close Pperm
423 
424  // count, how many row permutations are missing on current proc
425  size_t cntFreeRowIdx = 0;
426  std::queue<GlobalOrdinal> qFreeGRowIdx; // store global row ids of "free" rows
427  for (size_t lik = 0; lik < RowIdStatus->getLocalLength(); ++lik) {
428  if(RowIdStatusArray[lik] == SC_ZERO) {
429  cntFreeRowIdx++;
430  qFreeGRowIdx.push(RowIdStatus->getMap()->getGlobalElement(lik));
431  }
432  }
433 
434  // fix Pperm
435  for (size_t lik = 0; lik < RowIdStatus->getLocalLength(); ++lik) {
436  if(RowIdStatusArray[lik] == SC_ZERO) {
437  RowIdStatusArray[lik] = SC_ONE; // use this row id
438  Pperm->replaceLocalValue(lik, qFreeGRowIdx.front());
439  // detect wide range permutations
440  if(floor(qFreeGRowIdx.front()/nDofsPerNode) != floor(RowIdStatus->getMap()->getGlobalElement(lik)/nDofsPerNode)) {
441  lWideRangeRowPermutations++;
442  }
443  qFreeGRowIdx.pop();
444  }
445  }
446 
447  // close Qperm (free permutation entries in Qperm)
448  size_t cntFreeColIdx = 0;
449  std::queue<GlobalOrdinal> qFreeGColIdx; // store global column ids of "free" available columns
450  for (size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
451  if(ColIdStatusArray[ljk] == SC_ZERO) {
452  cntFreeColIdx++;
453  qFreeGColIdx.push(ColIdStatus->getMap()->getGlobalElement(ljk));
454  }
455  }
456 
457  size_t cntUnusedColIdx = 0;
458  std::queue<GlobalOrdinal> qUnusedGColIdx; // store global column ids of "free" available columns
459  for (size_t ljk = 0; ljk < ColIdUsed->getLocalLength(); ++ljk) {
460  if(ColIdUsedArray[ljk] == SC_ZERO) {
461  cntUnusedColIdx++;
462  qUnusedGColIdx.push(ColIdUsed->getMap()->getGlobalElement(ljk));
463  }
464  }
465 
466  // fix Qperm with local entries
467  for (size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
468  // stop if no (local) unused column idx are left
469  if(cntUnusedColIdx == 0) break;
470 
471  if(ColIdStatusArray[ljk] == SC_ZERO) {
472  ColIdStatusArray[ljk] = SC_ONE; // use this row id
473  Qperm->replaceLocalValue(ljk, qUnusedGColIdx.front()); // loop over ColIdStatus (lives on domain map)
474  ColIdUsed->replaceGlobalValue(qUnusedGColIdx.front(),SC_ONE); // ljk is now used, too
475  // detect wide range permutations
476  if(floor(qUnusedGColIdx.front()/nDofsPerNode) != floor(ColIdStatus->getMap()->getGlobalElement(ljk)/nDofsPerNode)) {
477  lWideRangeColPermutations++;
478  }
479  qUnusedGColIdx.pop();
480  cntUnusedColIdx--;
481  cntFreeColIdx--;
482  }
483  }
484 
485  //Qperm->doExport(*lQperm,*QpermExporter,Xpetra::ABSMAX); // no export necessary, since changes only locally
486  //ColIdStatus->doExport(*lColIdStatus,*QpermExporter,Xpetra::ABSMAX);
487 
488  // count, how many unused column idx are needed on current processor
489  // to complete Qperm
490  cntFreeColIdx = 0;
491  for (size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) { // TODO avoid this loop
492  if(ColIdStatusArray[ljk] == SC_ZERO) {
493  cntFreeColIdx++;
494  }
495  }
496 
497  GlobalOrdinal global_cntFreeColIdx = 0;
498  LocalOrdinal local_cntFreeColIdx = cntFreeColIdx;
499  MueLu_sumAll(comm, Teuchos::as<GlobalOrdinal>(local_cntFreeColIdx), global_cntFreeColIdx);
500 #ifdef DEBUG_OUTPUT
501  std::cout << "global # of empty column idx entries in Qperm: " << global_cntFreeColIdx << std::endl;
502 #endif
503 
504  // avoid global communication if possible
505  if(global_cntFreeColIdx > 0) {
506 
507  // 1) count how many unused column ids are left
508  GlobalOrdinal global_cntUnusedColIdx = 0;
509  LocalOrdinal local_cntUnusedColIdx = cntUnusedColIdx;
510  MueLu_sumAll(comm, Teuchos::as<GlobalOrdinal>(local_cntUnusedColIdx), global_cntUnusedColIdx);
511 #ifdef DEBUG_OUTPUT
512  std::cout << "global # of unused column idx: " << global_cntUnusedColIdx << std::endl;
513 #endif
514 
515  // 2) communicate how many unused column ids are available on procs
516  std::vector<LocalOrdinal> local_UnusedColIdxOnProc (numProcs);
517  std::vector<LocalOrdinal> global_UnusedColIdxOnProc(numProcs);
518  local_UnusedColIdxOnProc[myRank] = local_cntUnusedColIdx;
519  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
520  local_UnusedColIdxOnProc.data(),
521  global_UnusedColIdxOnProc.data());
522 
523 #ifdef DEBUG_OUTPUT
524  std::cout << "PROC " << myRank << " global num unused indices per proc: ";
525  for (size_t ljk = 0; ljk < global_UnusedColIdxOnProc.size(); ++ljk) {
526  std::cout << " " << global_UnusedColIdxOnProc[ljk];
527  }
528  std::cout << std::endl;
529 #endif
530 
531  // 3) build array of length global_cntUnusedColIdx to globally replicate unused column idx
532  std::vector<GlobalOrdinal> local_UnusedColIdxVector(Teuchos::as<size_t>(global_cntUnusedColIdx));
533  std::vector<GlobalOrdinal> global_UnusedColIdxVector(Teuchos::as<size_t>(global_cntUnusedColIdx));
534  GlobalOrdinal global_cntUnusedColIdxStartIter = 0;
535  for(int proc=0; proc<myRank; proc++) {
536  global_cntUnusedColIdxStartIter += global_UnusedColIdxOnProc[proc];
537  }
538  for(GlobalOrdinal k = global_cntUnusedColIdxStartIter; k < global_cntUnusedColIdxStartIter+local_cntUnusedColIdx; k++) {
539  local_UnusedColIdxVector[k] = qUnusedGColIdx.front();
540  qUnusedGColIdx.pop();
541  }
542  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(global_cntUnusedColIdx),
543  local_UnusedColIdxVector.data(),
544  global_UnusedColIdxVector.data());
545 #ifdef DEBUG_OUTPUT
546  std::cout << "PROC " << myRank << " global UnusedGColIdx: ";
547  for (size_t ljk = 0; ljk < global_UnusedColIdxVector.size(); ++ljk) {
548  std::cout << " " << global_UnusedColIdxVector[ljk];
549  }
550  std::cout << std::endl;
551 #endif
552 
553 
554 
555  // 4) communicate, how many column idx are needed on each processor
556  // to complete Qperm
557  std::vector<LocalOrdinal> local_EmptyColIdxOnProc (numProcs);
558  std::vector<LocalOrdinal> global_EmptyColIdxOnProc(numProcs);
559  local_EmptyColIdxOnProc[myRank] = local_cntFreeColIdx;
560  Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
561  local_EmptyColIdxOnProc.data(),
562  global_EmptyColIdxOnProc.data());
563 
564 #ifdef DEBUG_OUTPUT
565  std::cout << "PROC " << myRank << " global num of needed column indices: ";
566  for (size_t ljk = 0; ljk < global_EmptyColIdxOnProc.size(); ++ljk) {
567  std::cout << " " << global_EmptyColIdxOnProc[ljk];
568  }
569  std::cout << std::endl;
570 #endif
571 
572  // 5) determine first index in global_UnusedColIdxVector for unused column indices,
573  // that are marked to be used by this processor
574  GlobalOrdinal global_UnusedColStartIdx = 0;
575  for(int proc=0; proc<myRank; proc++) {
576  global_UnusedColStartIdx += global_EmptyColIdxOnProc[proc];
577  }
578 
579 #ifdef DEBUG_OUTPUT
580  GetOStream(Statistics0) << "PROC " << myRank << " is allowd to use the following column gids: ";
581  for(GlobalOrdinal k = global_UnusedColStartIdx; k < global_UnusedColStartIdx + Teuchos::as<GlobalOrdinal>(cntFreeColIdx); k++) {
582  GetOStream(Statistics0) << global_UnusedColIdxVector[k] << " ";
583  }
584  GetOStream(Statistics0) << std::endl;
585 #endif
586 
587  // 6.) fix Qperm with global entries
588  GlobalOrdinal array_iter = 0;
589  for (size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
590 
591  if(ColIdStatusArray[ljk] == SC_ZERO) {
592  ColIdStatusArray[ljk] = SC_ONE; // use this row id
593  Qperm->replaceLocalValue(ljk, global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter]);
594  ColIdUsed->replaceGlobalValue(global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter],SC_ONE);
595  // detect wide range permutations
596  if(floor(global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter]/nDofsPerNode) != floor(ColIdStatus->getMap()->getGlobalElement(ljk)/nDofsPerNode)) {
597  lWideRangeColPermutations++;
598  }
599  array_iter++;
600  //cntUnusedColIdx--; // check me
601  }
602  }
603  } // end if global_cntFreeColIdx > 0
605 
606 
607  // create new empty Matrix
608  Teuchos::RCP<CrsMatrixWrap> permPTmatrix = Teuchos::rcp(new CrsMatrixWrap(A->getRowMap(),1,Xpetra::StaticProfile));
609  Teuchos::RCP<CrsMatrixWrap> permQTmatrix = Teuchos::rcp(new CrsMatrixWrap(A->getRowMap(),1,Xpetra::StaticProfile));
610 
611  for(size_t row=0; row<A->getNodeNumRows(); row++) {
612  // FIXME (mfh 30 Oct 2015): Teuchos::as doesn't know how to
613  // convert from complex Scalar to GO, so we have to take the real
614  // part first. I think that's the right thing to do in this case.
615  Teuchos::ArrayRCP<GlobalOrdinal> indoutP(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(PpermData[row]))); // column idx for Perm^T
616  Teuchos::ArrayRCP<GlobalOrdinal> indoutQ(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(QpermData[row]))); // column idx for Qperm
617  Teuchos::ArrayRCP<Scalar> valout(1,SC_ONE);
618  permPTmatrix->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indoutP.view(0,indoutP.size()), valout.view(0,valout.size()));
619  permQTmatrix->insertGlobalValues (A->getRowMap()->getGlobalElement(row), indoutQ.view(0,indoutQ.size()), valout.view(0,valout.size()));
620  }
621 
622  permPTmatrix->fillComplete();
623  permQTmatrix->fillComplete();
624 
625  Teuchos::RCP<Matrix> permPmatrix = Utilities::Transpose(*permPTmatrix, true);
626 
627  for(size_t row=0; row<permPTmatrix->getNodeNumRows(); row++) {
628  if(permPTmatrix->getNumEntriesInLocalRow(row) != 1)
629  GetOStream(Warnings0) <<"#entries in row " << row << " of permPTmatrix is " << permPTmatrix->getNumEntriesInLocalRow(row) << std::endl;
630  if(permPmatrix->getNumEntriesInLocalRow(row) != 1)
631  GetOStream(Warnings0) <<"#entries in row " << row << " of permPmatrix is " << permPmatrix->getNumEntriesInLocalRow(row) << std::endl;
632  if(permQTmatrix->getNumEntriesInLocalRow(row) != 1)
633  GetOStream(Warnings0) <<"#entries in row " << row << " of permQmatrix is " << permQTmatrix->getNumEntriesInLocalRow(row) << std::endl;
634  }
635 
636  // build permP * A * permQT
638  Teuchos::RCP<Matrix> permPApermQt = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*permPmatrix, false, *ApermQt, false, GetOStream(Statistics2));
639 
640  /*
641  MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("A.mat", *A);
642  MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permP.mat", *permPmatrix);
643  MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permQt.mat", *permQTmatrix);
644  MueLu::Utils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Write("permPApermQt.mat", *permPApermQt);
645  */
646  // build scaling matrix
647  Teuchos::RCP<Vector> diagVec = VectorFactory::Build(permPApermQt->getRowMap(),true);
648  Teuchos::RCP<Vector> invDiagVec = VectorFactory::Build(permPApermQt->getRowMap(),true);
649  Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
650  Teuchos::ArrayRCP< Scalar > invDiagVecData = invDiagVec->getDataNonConst(0);
651 
652  permPApermQt->getLocalDiagCopy(*diagVec);
653  for(size_t i = 0; i<diagVec->getMap()->getNodeNumElements(); ++i) {
654  if(diagVecData[i] != SC_ZERO)
655  invDiagVecData[i] = SC_ONE / diagVecData[i];
656  else {
657  invDiagVecData[i] = SC_ONE;
658  GetOStream(Statistics0) << "MueLu::PermutationFactory: found zero on diagonal in row " << i << std::endl;
659  }
660  }
661 
662  Teuchos::RCP<CrsMatrixWrap> diagScalingOp = Teuchos::rcp(new CrsMatrixWrap(permPApermQt->getRowMap(),1,Xpetra::StaticProfile));
663 
664  for(size_t row=0; row<A->getNodeNumRows(); row++) {
665  Teuchos::ArrayRCP<GlobalOrdinal> indout(1,permPApermQt->getRowMap()->getGlobalElement(row)); // column idx for Perm^T
666  Teuchos::ArrayRCP<Scalar> valout(1,invDiagVecData[row]);
667  diagScalingOp->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
668  }
669  diagScalingOp->fillComplete();
670 
671  Teuchos::RCP<Matrix> scaledA = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*diagScalingOp, false, *permPApermQt, false, GetOStream(Statistics2));
672  currentLevel.Set("A", Teuchos::rcp_dynamic_cast<Matrix>(scaledA), genFactory/*this*/);
673 
674  currentLevel.Set("permA", Teuchos::rcp_dynamic_cast<Matrix>(permPApermQt), genFactory/*this*/); // TODO careful with this!!!
675  currentLevel.Set("permP", Teuchos::rcp_dynamic_cast<Matrix>(permPmatrix), genFactory/*this*/);
676  currentLevel.Set("permQT", Teuchos::rcp_dynamic_cast<Matrix>(permQTmatrix), genFactory/*this*/);
677  currentLevel.Set("permScaling", Teuchos::rcp_dynamic_cast<Matrix>(diagScalingOp), genFactory/*this*/);
678 
680  // count zeros on diagonal in P -> number of row permutations
681  Teuchos::RCP<Vector> diagPVec = VectorFactory::Build(permPmatrix->getRowMap(),true);
682  permPmatrix->getLocalDiagCopy(*diagPVec);
683  Teuchos::ArrayRCP< const Scalar > diagPVecData = diagPVec->getData(0);
684  LocalOrdinal lNumRowPermutations = 0;
685  GlobalOrdinal gNumRowPermutations = 0;
686  for(size_t i = 0; i<diagPVec->getMap()->getNodeNumElements(); ++i) {
687  if(diagPVecData[i] == SC_ZERO) {
688  lNumRowPermutations++;
689  }
690  }
691 
692  // sum up all entries in multipleColRequests over all processors
693  MueLu_sumAll(diagPVec->getMap()->getComm(), Teuchos::as<GlobalOrdinal>(lNumRowPermutations), gNumRowPermutations);
694 
696  // count zeros on diagonal in Q^T -> number of column permutations
697  Teuchos::RCP<Vector> diagQTVec = VectorFactory::Build(permQTmatrix->getRowMap(),true);
698  permQTmatrix->getLocalDiagCopy(*diagQTVec);
699  Teuchos::ArrayRCP< const Scalar > diagQTVecData = diagQTVec->getData(0);
700  LocalOrdinal lNumColPermutations = 0;
701  GlobalOrdinal gNumColPermutations = 0;
702  for(size_t i = 0; i<diagQTVec->getMap()->getNodeNumElements(); ++i) {
703  if(diagQTVecData[i] == SC_ZERO) {
704  lNumColPermutations++;
705  }
706  }
707 
708  // sum up all entries in multipleColRequests over all processors
709  MueLu_sumAll(diagQTVec->getMap()->getComm(), Teuchos::as<GlobalOrdinal>(lNumColPermutations), gNumColPermutations);
710 
711  currentLevel.Set("#RowPermutations", gNumRowPermutations, genFactory/*this*/);
712  currentLevel.Set("#ColPermutations", gNumColPermutations, genFactory/*this*/);
713  currentLevel.Set("#WideRangeRowPermutations", gWideRangeRowPermutations, genFactory/*this*/);
714  currentLevel.Set("#WideRangeColPermutations", gWideRangeColPermutations, genFactory/*this*/);
715 
716  GetOStream(Statistics0) << "#Row permutations/max possible permutations: " << gNumRowPermutations << "/" << diagPVec->getMap()->getGlobalNumElements() << std::endl;
717  GetOStream(Statistics0) << "#Column permutations/max possible permutations: " << gNumColPermutations << "/" << diagQTVec->getMap()->getGlobalNumElements() << std::endl;
718  GetOStream(Runtime1) << "#wide range row permutations: " << gWideRangeRowPermutations << " #wide range column permutations: " << gWideRangeColPermutations << std::endl;
719  }
720 
721 } // namespace MueLu
722 
723 #endif /* MUELU_ALGEBRAICPERMUTATIONSTRATEGY_DEF_HPP_ */
Important warning messages (one line)
void BuildPermutation(const Teuchos::RCP< Matrix > &A, const Teuchos::RCP< const Map > &permRowMap, Level &currentLevel, const FactoryBase *genFactory) const
build permutation operators
#define MueLu_sumAll(rcpComm, in, out)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void sortingPermutation(const std::vector< Scalar > &values, std::vector< LocalOrdinal > &v)
size_type size() const
One-liner description of what is happening.
size_type size() 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.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
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)
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
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 > &params=null)
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
iterator begin() const
ArrayView< T > view(size_type lowerOffset, size_type size) const