8 #ifndef MUELU_ALGEBRAICPERMUTATIONSTRATEGY_DEF_HPP_ 9 #define MUELU_ALGEBRAICPERMUTATIONSTRATEGY_DEF_HPP_ 13 #include <Xpetra_MultiVector.hpp> 14 #include <Xpetra_Matrix.hpp> 15 #include <Xpetra_CrsGraph.hpp> 16 #include <Xpetra_Vector.hpp> 17 #include <Xpetra_VectorFactory.hpp> 18 #include <Xpetra_CrsMatrixWrap.hpp> 19 #include <Xpetra_Export.hpp> 20 #include <Xpetra_ExportFactory.hpp> 21 #include <Xpetra_Import.hpp> 22 #include <Xpetra_ImportFactory.hpp> 23 #include <Xpetra_MatrixMatrix.hpp> 25 #include "MueLu_Utilities.hpp" 30 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
33 const Teuchos::RCP< const Teuchos::Comm< int > > comm = A->getRowMap()->getComm();
34 int numProcs = comm->getSize();
35 int myRank = comm->getRank();
37 size_t nDofsPerNode = 1;
38 if (A->IsView(
"stridedMaps")) {
39 Teuchos::RCP<const Map> permRowMapStrided = A->getRowMap(
"stridedMaps");
40 nDofsPerNode = Teuchos::rcp_dynamic_cast<
const StridedMap>(permRowMapStrided)->getFixedBlockSize();
43 std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > permutedDiagCandidates;
44 std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > keepDiagonalEntries;
45 std::vector<Scalar> Weights;
49 for (
size_t row = 0; row < A->getRowMap()->getNodeNumElements(); row++) {
50 GlobalOrdinal grow = A->getRowMap()->getGlobalElement(row);
52 if(permRowMap->isNodeGlobalElement(grow) ==
true)
continue;
54 size_t nnz = A->getNumEntriesInLocalRow(row);
57 Teuchos::ArrayView<const LocalOrdinal> indices;
58 Teuchos::ArrayView<const Scalar> vals;
59 A->getLocalRowView(row, indices, vals);
61 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz,
Exceptions::RuntimeError,
"MueLu::PermutationFactory::Build: number of nonzeros not equal to number of indices? Error.");
64 GlobalOrdinal gMaxValIdx = 0;
65 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
66 MT norm1 = Teuchos::ScalarTraits<MT>::zero ();
67 MT maxVal = Teuchos::ScalarTraits<MT>::zero ();
68 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
69 norm1 += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
70 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
71 maxVal = Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
72 gMaxValIdx = A->getColMap()->getGlobalElement(indices[j]);
76 if(grow == gMaxValIdx)
77 keepDiagonalEntries.push_back(std::make_pair(grow,grow));
82 for (
size_t row = 0; row < permRowMap->getNodeNumElements(); row++) {
83 GlobalOrdinal grow = permRowMap->getGlobalElement(row);
84 LocalOrdinal lArow = A->getRowMap()->getLocalElement(grow);
85 size_t nnz = A->getNumEntriesInLocalRow(lArow);
88 Teuchos::ArrayView<const LocalOrdinal> indices;
89 Teuchos::ArrayView<const Scalar> vals;
90 A->getLocalRowView(lArow, indices, vals);
92 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz,
Exceptions::RuntimeError,
"MueLu::PermutationFactory::Build: number of nonzeros not equal to number of indices? Error.");
95 GlobalOrdinal gMaxValIdx = 0;
96 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
97 MT norm1 = Teuchos::ScalarTraits<MT>::zero ();
98 MT maxVal = Teuchos::ScalarTraits<MT>::zero ();
99 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
100 norm1 += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
101 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
102 maxVal = Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
103 gMaxValIdx = A->getColMap()->getGlobalElement(indices[j]);
107 if(maxVal > Teuchos::ScalarTraits<MT>::zero ()) {
108 permutedDiagCandidates.push_back(std::make_pair(grow,gMaxValIdx));
109 Weights.push_back(maxVal/(norm1*Teuchos::as<Scalar>(nnz)));
111 std::cout <<
"ATTENTION: row " << grow <<
" has only zero entries -> singular matrix!" << std::endl;
117 std::vector<int> permutation;
126 Teuchos::RCP<Vector> gColVec = VectorFactory::Build(A->getColMap());
127 Teuchos::RCP<Vector> gDomVec = VectorFactory::Build(A->getDomainMap());
128 gColVec->putScalar(0.0);
129 gDomVec->putScalar(0.0);
132 for (
typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::const_iterator p = keepDiagonalEntries.begin(); p != keepDiagonalEntries.end(); ++p) {
133 gColVec->sumIntoGlobalValue((*p).second,1.0);
136 Teuchos::RCP<Export> exporter = ExportFactory::Build(gColVec->getMap(), gDomVec->getMap());
137 gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD);
138 gColVec->doImport(*gDomVec,*exporter,Xpetra::INSERT);
140 std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > permutedDiagCandidatesFiltered;
141 std::map<GlobalOrdinal, Scalar> gColId2Weight;
143 Teuchos::ArrayRCP< Scalar > ddata = gColVec->getDataNonConst(0);
144 for(
size_t i = 0; i < permutedDiagCandidates.size(); ++i) {
146 std::pair<GlobalOrdinal, GlobalOrdinal> pp = permutedDiagCandidates[permutation[i]];
147 GlobalOrdinal grow = pp.first;
148 GlobalOrdinal gcol = pp.second;
150 LocalOrdinal lcol = A->getColMap()->getLocalElement(gcol);
151 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
152 if(Teuchos::ScalarTraits<Scalar>::real (ddata[lcol]) > Teuchos::ScalarTraits<MT>::zero ()){
157 ddata[lcol] += Teuchos::ScalarTraits<Scalar>::one ();
159 permutedDiagCandidatesFiltered.push_back(std::make_pair(grow,gcol));
160 gColId2Weight[gcol] = Weights[permutation[i]];
164 gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD);
165 gColVec->doImport(*gDomVec,*exporter,Xpetra::INSERT);
172 std::vector<GlobalOrdinal> multipleColRequests;
176 std::queue<GlobalOrdinal> unusedColIdx;
178 for(
size_t sz = 0; sz<gDomVec->getLocalLength(); ++sz) {
179 Teuchos::ArrayRCP< const Scalar > arrDomVec = gDomVec->getData(0);
185 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
186 if(Teuchos::ScalarTraits<Scalar>::real (arrDomVec[sz]) > Teuchos::ScalarTraits<MT>::one ()) {
187 multipleColRequests.push_back(gDomVec->getMap()->getGlobalElement(sz));
188 }
else if(Teuchos::ScalarTraits<Scalar>::real (arrDomVec[sz]) ==
189 Teuchos::ScalarTraits<MT>::zero ()) {
190 unusedColIdx.push(gDomVec->getMap()->getGlobalElement(sz));
195 LocalOrdinal localMultColRequests = Teuchos::as<LocalOrdinal>(multipleColRequests.size());
196 LocalOrdinal globalMultColRequests = 0;
199 MueLu_sumAll(gDomVec->getMap()->getComm(), (LocalOrdinal)localMultColRequests, globalMultColRequests);
201 if(globalMultColRequests > 0) {
207 std::vector<GlobalOrdinal> numMyMultColRequests(numProcs,0);
208 std::vector<GlobalOrdinal> numGlobalMultColRequests(numProcs,0);
209 numMyMultColRequests[myRank] = localMultColRequests;
210 Teuchos::reduceAll(*comm,Teuchos::REDUCE_MAX,numProcs,&numMyMultColRequests[0],&numGlobalMultColRequests[0]);
214 for (
int i=0; i<myRank-1; i++)
215 nMyOffset += numGlobalMultColRequests[i];
217 GlobalOrdinal zero=0;
218 std::vector<GlobalOrdinal> procMultRequestedColIds(globalMultColRequests,zero);
219 std::vector<GlobalOrdinal> global_procMultRequestedColIds(globalMultColRequests,zero);
222 for(
size_t i = 0; i < multipleColRequests.size(); i++) {
223 procMultRequestedColIds[nMyOffset + i] = multipleColRequests[i];
227 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(globalMultColRequests), &procMultRequestedColIds[0], &global_procMultRequestedColIds[0]);
230 for (
size_t k = 0; k<global_procMultRequestedColIds.size(); k++) {
231 GlobalOrdinal globColId = global_procMultRequestedColIds[k];
233 std::vector<Scalar> MyWeightForColId(numProcs,0);
234 std::vector<Scalar> GlobalWeightForColId(numProcs,0);
236 if(gColVec->getMap()->isNodeGlobalElement(globColId)) {
237 MyWeightForColId[myRank] = gColId2Weight[globColId];
239 MyWeightForColId[myRank] = 0.0;
242 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs, &MyWeightForColId[0], &GlobalWeightForColId[0]);
244 if(gColVec->getMap()->isNodeGlobalElement(globColId)) {
247 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
248 MT winnerValue = 0.0;
249 int winnerProcRank = 0;
250 for (
int proc = 0; proc < numProcs; proc++) {
251 if(Teuchos::ScalarTraits<Scalar>::real (GlobalWeightForColId[proc]) > winnerValue) {
252 winnerValue = Teuchos::ScalarTraits<Scalar>::real (GlobalWeightForColId[proc]);
253 winnerProcRank = proc;
260 if(myRank != winnerProcRank) {
262 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = permutedDiagCandidatesFiltered.begin();
263 while(p != permutedDiagCandidatesFiltered.end() )
265 if((*p).second == globColId)
266 p = permutedDiagCandidatesFiltered.erase(p);
278 std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > RowColPairs;
279 RowColPairs.insert( RowColPairs.end(), keepDiagonalEntries.begin(), keepDiagonalEntries.end());
280 RowColPairs.insert( RowColPairs.end(), permutedDiagCandidatesFiltered.begin(), permutedDiagCandidatesFiltered.end());
285 gColVec->putScalar(0.0);
286 gDomVec->putScalar(0.0);
287 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator pl = RowColPairs.begin();
288 while(pl != RowColPairs.end() )
291 GlobalOrdinal jk = (*pl).second;
293 gColVec->sumIntoGlobalValue(jk,1.0);
296 gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD);
297 for(
size_t sz = 0; sz<gDomVec->getLocalLength(); ++sz) {
298 Teuchos::ArrayRCP< const Scalar > arrDomVec = gDomVec->getData(0);
299 if(arrDomVec[sz] > 1.0) {
300 GetOStream(
Runtime0) <<
"RowColPairs has multiple column [" << sz <<
"]=" << arrDomVec[sz] << std::endl;
301 }
else if(arrDomVec[sz] == 0.0) {
302 GetOStream(
Runtime0) <<
"RowColPairs has empty column [" << sz <<
"]=" << arrDomVec[sz] << std::endl;
335 Teuchos::RCP<Vector> Pperm = VectorFactory::Build(A->getRowMap());
336 Teuchos::RCP<Vector> Qperm = VectorFactory::Build(A->getDomainMap());
337 Teuchos::RCP<Vector> lQperm = VectorFactory::Build(A->getColMap());
339 Teuchos::ArrayRCP< Scalar > PpermData = Pperm->getDataNonConst(0);
340 Teuchos::ArrayRCP< Scalar > QpermData = Qperm->getDataNonConst(0);
342 Pperm->putScalar(0.0);
343 Qperm->putScalar(0.0);
344 lQperm->putScalar(0.0);
347 Teuchos::RCP<Export> QpermExporter = ExportFactory::Build(lQperm->getMap(), Qperm->getMap());
349 Teuchos::RCP<Vector> RowIdStatus = VectorFactory::Build(A->getRowMap());
350 Teuchos::RCP<Vector> ColIdStatus = VectorFactory::Build(A->getDomainMap());
351 Teuchos::RCP<Vector> lColIdStatus = VectorFactory::Build(A->getColMap());
352 Teuchos::RCP<Vector> ColIdUsed = VectorFactory::Build(A->getDomainMap());
353 Teuchos::ArrayRCP< Scalar > RowIdStatusArray = RowIdStatus->getDataNonConst(0);
354 Teuchos::ArrayRCP< Scalar > ColIdStatusArray = ColIdStatus->getDataNonConst(0);
355 Teuchos::ArrayRCP< Scalar > lColIdStatusArray = lColIdStatus->getDataNonConst(0);
356 Teuchos::ArrayRCP< Scalar > ColIdUsedArray = ColIdUsed->getDataNonConst(0);
357 RowIdStatus->putScalar(0.0);
358 ColIdStatus->putScalar(0.0);
359 lColIdStatus->putScalar(0.0);
360 ColIdUsed->putScalar(0.0);
365 LocalOrdinal lWideRangeRowPermutations = 0;
366 GlobalOrdinal gWideRangeRowPermutations = 0;
367 LocalOrdinal lWideRangeColPermutations = 0;
368 GlobalOrdinal gWideRangeColPermutations = 0;
371 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = RowColPairs.begin();
372 while(p != RowColPairs.end() )
374 GlobalOrdinal ik = (*p).first;
375 GlobalOrdinal jk = (*p).second;
377 LocalOrdinal lik = A->getRowMap()->getLocalElement(ik);
378 LocalOrdinal ljk = A->getColMap()->getLocalElement(jk);
380 if(RowIdStatusArray[lik] == 0.0) {
381 RowIdStatusArray[lik] = 1.0;
382 lColIdStatusArray[ljk] = 1.0;
383 Pperm->replaceLocalValue(lik, ik);
384 lQperm->replaceLocalValue(ljk, ik);
385 ColIdUsed->replaceGlobalValue(ik,1.0);
386 p = RowColPairs.erase(p);
389 if(floor(ik/nDofsPerNode) != floor(jk/nDofsPerNode)) {
390 lWideRangeColPermutations++;
398 Qperm->doExport(*lQperm,*QpermExporter,Xpetra::ABSMAX);
399 ColIdStatus->doExport(*lColIdStatus,*QpermExporter,Xpetra::ABSMAX);
402 if(RowColPairs.size()>0) GetOStream(
Warnings0) <<
"MueLu::PermutationFactory: There are Row/Col pairs left!!!" << std::endl;
407 size_t cntFreeRowIdx = 0;
408 std::queue<GlobalOrdinal> qFreeGRowIdx;
409 for (
size_t lik = 0; lik < RowIdStatus->getLocalLength(); ++lik) {
410 if(RowIdStatusArray[lik] == 0.0) {
412 qFreeGRowIdx.push(RowIdStatus->getMap()->getGlobalElement(lik));
417 for (
size_t lik = 0; lik < RowIdStatus->getLocalLength(); ++lik) {
418 if(RowIdStatusArray[lik] == 0.0) {
419 RowIdStatusArray[lik] = 1.0;
420 Pperm->replaceLocalValue(lik, qFreeGRowIdx.front());
422 if(floor(qFreeGRowIdx.front()/nDofsPerNode) != floor(RowIdStatus->getMap()->getGlobalElement(lik)/nDofsPerNode)) {
423 lWideRangeRowPermutations++;
430 size_t cntFreeColIdx = 0;
431 std::queue<GlobalOrdinal> qFreeGColIdx;
432 for (
size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
433 if(ColIdStatusArray[ljk] == 0.0) {
435 qFreeGColIdx.push(ColIdStatus->getMap()->getGlobalElement(ljk));
439 size_t cntUnusedColIdx = 0;
440 std::queue<GlobalOrdinal> qUnusedGColIdx;
441 for (
size_t ljk = 0; ljk < ColIdUsed->getLocalLength(); ++ljk) {
442 if(ColIdUsedArray[ljk] == 0.0) {
444 qUnusedGColIdx.push(ColIdUsed->getMap()->getGlobalElement(ljk));
449 for (
size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
451 if(cntUnusedColIdx == 0)
break;
453 if(ColIdStatusArray[ljk] == 0.0) {
454 ColIdStatusArray[ljk] = 1.0;
455 Qperm->replaceLocalValue(ljk, qUnusedGColIdx.front());
456 ColIdUsed->replaceGlobalValue(qUnusedGColIdx.front(),1.0);
458 if(floor(qUnusedGColIdx.front()/nDofsPerNode) != floor(ColIdStatus->getMap()->getGlobalElement(ljk)/nDofsPerNode)) {
459 lWideRangeColPermutations++;
461 qUnusedGColIdx.pop();
473 for (
size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
474 if(ColIdStatusArray[ljk] == 0.0) {
479 GlobalOrdinal global_cntFreeColIdx = 0;
480 LocalOrdinal local_cntFreeColIdx = cntFreeColIdx;
481 MueLu_sumAll(comm, Teuchos::as<GlobalOrdinal>(local_cntFreeColIdx), global_cntFreeColIdx);
483 std::cout <<
"global # of empty column idx entries in Qperm: " << global_cntFreeColIdx << std::endl;
487 if(global_cntFreeColIdx > 0) {
490 GlobalOrdinal global_cntUnusedColIdx = 0;
491 LocalOrdinal local_cntUnusedColIdx = cntUnusedColIdx;
492 MueLu_sumAll(comm, Teuchos::as<GlobalOrdinal>(local_cntUnusedColIdx), global_cntUnusedColIdx);
494 std::cout <<
"global # of unused column idx: " << global_cntUnusedColIdx << std::endl;
498 std::vector<LocalOrdinal> local_UnusedColIdxOnProc (numProcs);
499 std::vector<LocalOrdinal> global_UnusedColIdxOnProc(numProcs);
500 local_UnusedColIdxOnProc[myRank] = local_cntUnusedColIdx;
501 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs, &local_UnusedColIdxOnProc[0], &global_UnusedColIdxOnProc[0]);
504 std::cout <<
"PROC " << myRank <<
" global num unused indices per proc: ";
505 for (
size_t ljk = 0; ljk < global_UnusedColIdxOnProc.size(); ++ljk) {
506 std::cout <<
" " << global_UnusedColIdxOnProc[ljk];
508 std::cout << std::endl;
512 std::vector<GlobalOrdinal> local_UnusedColIdxVector(Teuchos::as<size_t>(global_cntUnusedColIdx));
513 std::vector<GlobalOrdinal> global_UnusedColIdxVector(Teuchos::as<size_t>(global_cntUnusedColIdx));
514 GlobalOrdinal global_cntUnusedColIdxStartIter = 0;
515 for(
int proc=0; proc<myRank; proc++) {
516 global_cntUnusedColIdxStartIter += global_UnusedColIdxOnProc[proc];
518 for(GlobalOrdinal k = global_cntUnusedColIdxStartIter; k < global_cntUnusedColIdxStartIter+local_cntUnusedColIdx; k++) {
519 local_UnusedColIdxVector[k] = qUnusedGColIdx.front();
520 qUnusedGColIdx.pop();
522 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(global_cntUnusedColIdx), &local_UnusedColIdxVector[0], &global_UnusedColIdxVector[0]);
524 std::cout <<
"PROC " << myRank <<
" global UnusedGColIdx: ";
525 for (
size_t ljk = 0; ljk < global_UnusedColIdxVector.size(); ++ljk) {
526 std::cout <<
" " << global_UnusedColIdxVector[ljk];
528 std::cout << std::endl;
535 std::vector<LocalOrdinal> local_EmptyColIdxOnProc (numProcs);
536 std::vector<LocalOrdinal> global_EmptyColIdxOnProc(numProcs);
537 local_EmptyColIdxOnProc[myRank] = local_cntFreeColIdx;
538 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs, &local_EmptyColIdxOnProc[0], &global_EmptyColIdxOnProc[0]);
541 std::cout <<
"PROC " << myRank <<
" global num of needed column indices: ";
542 for (
size_t ljk = 0; ljk < global_EmptyColIdxOnProc.size(); ++ljk) {
543 std::cout <<
" " << global_EmptyColIdxOnProc[ljk];
545 std::cout << std::endl;
550 GlobalOrdinal global_UnusedColStartIdx = 0;
551 for(
int proc=0; proc<myRank; proc++) {
552 global_UnusedColStartIdx += global_EmptyColIdxOnProc[proc];
556 GetOStream(
Statistics0) <<
"PROC " << myRank <<
" is allowd to use the following column gids: ";
557 for(GlobalOrdinal k = global_UnusedColStartIdx; k < global_UnusedColStartIdx + Teuchos::as<GlobalOrdinal>(cntFreeColIdx); k++) {
558 GetOStream(
Statistics0) << global_UnusedColIdxVector[k] <<
" ";
564 GlobalOrdinal array_iter = 0;
565 for (
size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
567 if(ColIdStatusArray[ljk] == 0.0) {
568 ColIdStatusArray[ljk] = 1.0;
569 Qperm->replaceLocalValue(ljk, global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter]);
570 ColIdUsed->replaceGlobalValue(global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter],1.0);
572 if(floor(global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter]/nDofsPerNode) != floor(ColIdStatus->getMap()->getGlobalElement(ljk)/nDofsPerNode)) {
573 lWideRangeColPermutations++;
584 Teuchos::RCP<CrsMatrixWrap> permPTmatrix = Teuchos::rcp(
new CrsMatrixWrap(A->getRowMap(),1,Xpetra::StaticProfile));
585 Teuchos::RCP<CrsMatrixWrap> permQTmatrix = Teuchos::rcp(
new CrsMatrixWrap(A->getRowMap(),1,Xpetra::StaticProfile));
587 for(
size_t row=0; row<A->getNodeNumRows(); row++) {
591 Teuchos::ArrayRCP<GlobalOrdinal> indoutP(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(PpermData[row])));
592 Teuchos::ArrayRCP<GlobalOrdinal> indoutQ(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(QpermData[row])));
593 Teuchos::ArrayRCP<Scalar> valout(1,Teuchos::ScalarTraits<Scalar>::one());
594 permPTmatrix->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indoutP.view(0,indoutP.size()), valout.view(0,valout.size()));
595 permQTmatrix->insertGlobalValues (A->getRowMap()->getGlobalElement(row), indoutQ.view(0,indoutQ.size()), valout.view(0,valout.size()));
598 permPTmatrix->fillComplete();
599 permQTmatrix->fillComplete();
603 for(
size_t row=0; row<permPTmatrix->getNodeNumRows(); row++) {
604 if(permPTmatrix->getNumEntriesInLocalRow(row) != 1)
605 GetOStream(
Warnings0) <<
"#entries in row " << row <<
" of permPTmatrix is " << permPTmatrix->getNumEntriesInLocalRow(row) << std::endl;
606 if(permPmatrix->getNumEntriesInLocalRow(row) != 1)
607 GetOStream(
Warnings0) <<
"#entries in row " << row <<
" of permPmatrix is " << permPmatrix->getNumEntriesInLocalRow(row) << std::endl;
608 if(permQTmatrix->getNumEntriesInLocalRow(row) != 1)
609 GetOStream(
Warnings0) <<
"#entries in row " << row <<
" of permQmatrix is " << permQTmatrix->getNumEntriesInLocalRow(row) << std::endl;
613 Teuchos::RCP<Matrix> ApermQt = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *permQTmatrix,
false, GetOStream(
Statistics2));
614 Teuchos::RCP<Matrix> permPApermQt = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*permPmatrix,
false, *ApermQt,
false, GetOStream(
Statistics2));
623 Teuchos::RCP<Vector> diagVec = VectorFactory::Build(permPApermQt->getRowMap(),
true);
624 Teuchos::RCP<Vector> invDiagVec = VectorFactory::Build(permPApermQt->getRowMap(),
true);
625 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
626 Teuchos::ArrayRCP< Scalar > invDiagVecData = invDiagVec->getDataNonConst(0);
628 permPApermQt->getLocalDiagCopy(*diagVec);
629 for(
size_t i = 0; i<diagVec->getMap()->getNodeNumElements(); ++i) {
630 if(diagVecData[i] != 0.0)
631 invDiagVecData[i] = Teuchos::ScalarTraits<Scalar>::one () / diagVecData[i];
633 invDiagVecData[i] = Teuchos::ScalarTraits<Scalar>::one ();
634 GetOStream(
Statistics0) <<
"MueLu::PermutationFactory: found zero on diagonal in row " << i << std::endl;
638 Teuchos::RCP<CrsMatrixWrap> diagScalingOp = Teuchos::rcp(
new CrsMatrixWrap(permPApermQt->getRowMap(),1,Xpetra::StaticProfile));
640 for(
size_t row=0; row<A->getNodeNumRows(); row++) {
641 Teuchos::ArrayRCP<GlobalOrdinal> indout(1,permPApermQt->getRowMap()->getGlobalElement(row));
642 Teuchos::ArrayRCP<Scalar> valout(1,invDiagVecData[row]);
643 diagScalingOp->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
645 diagScalingOp->fillComplete();
647 Teuchos::RCP<Matrix> scaledA = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*diagScalingOp,
false, *permPApermQt,
false, GetOStream(
Statistics2));
648 currentLevel.
Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(scaledA), genFactory);
650 currentLevel.
Set(
"permA", Teuchos::rcp_dynamic_cast<Matrix>(permPApermQt), genFactory);
651 currentLevel.
Set(
"permP", Teuchos::rcp_dynamic_cast<Matrix>(permPmatrix), genFactory);
652 currentLevel.
Set(
"permQT", Teuchos::rcp_dynamic_cast<Matrix>(permQTmatrix), genFactory);
653 currentLevel.
Set(
"permScaling", Teuchos::rcp_dynamic_cast<Matrix>(diagScalingOp), genFactory);
657 Teuchos::RCP<Vector> diagPVec = VectorFactory::Build(permPmatrix->getRowMap(),
true);
658 permPmatrix->getLocalDiagCopy(*diagPVec);
659 Teuchos::ArrayRCP< const Scalar > diagPVecData = diagPVec->getData(0);
660 LocalOrdinal lNumRowPermutations = 0;
661 GlobalOrdinal gNumRowPermutations = 0;
662 for(
size_t i = 0; i<diagPVec->getMap()->getNodeNumElements(); ++i) {
663 if(diagPVecData[i] == 0.0) {
664 lNumRowPermutations++;
669 MueLu_sumAll(diagPVec->getMap()->getComm(), Teuchos::as<GlobalOrdinal>(lNumRowPermutations), gNumRowPermutations);
673 Teuchos::RCP<Vector> diagQTVec = VectorFactory::Build(permQTmatrix->getRowMap(),
true);
674 permQTmatrix->getLocalDiagCopy(*diagQTVec);
675 Teuchos::ArrayRCP< const Scalar > diagQTVecData = diagQTVec->getData(0);
676 LocalOrdinal lNumColPermutations = 0;
677 GlobalOrdinal gNumColPermutations = 0;
678 for(
size_t i = 0; i<diagQTVec->getMap()->getNodeNumElements(); ++i) {
679 if(diagQTVecData[i] == 0.0) {
680 lNumColPermutations++;
685 MueLu_sumAll(diagQTVec->getMap()->getComm(), Teuchos::as<GlobalOrdinal>(lNumColPermutations), gNumColPermutations);
687 currentLevel.
Set(
"#RowPermutations", gNumRowPermutations, genFactory);
688 currentLevel.
Set(
"#ColPermutations", gNumColPermutations, genFactory);
689 currentLevel.
Set(
"#WideRangeRowPermutations", gWideRangeRowPermutations, genFactory);
690 currentLevel.
Set(
"#WideRangeColPermutations", gWideRangeColPermutations, genFactory);
692 GetOStream(
Statistics0) <<
"#Row permutations/max possible permutations: " << gNumRowPermutations <<
"/" << diagPVec->getMap()->getGlobalNumElements() << std::endl;
693 GetOStream(
Statistics0) <<
"#Column permutations/max possible permutations: " << gNumColPermutations <<
"/" << diagQTVec->getMap()->getGlobalNumElements() << std::endl;
694 GetOStream(
Runtime1) <<
"#wide range row permutations: " << gWideRangeRowPermutations <<
" #wide range column permutations: " << gWideRangeColPermutations << std::endl;
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())
void sortingPermutation(const std::vector< Scalar > &values, std::vector< LocalOrdinal > &v)
One-liner description of what is happening.
Namespace for MueLu classes and methods.
Print even more statistics.
Base class for factories (e.g., R, P, and A_coarse).
Print statistics that do not involve significant additional computation.
Class that holds all level-specific information.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void BuildPermutation(const Teuchos::RCP< Matrix > &A, const Teuchos::RCP< const Map > permRowMap, Level ¤tLevel, const FactoryBase *genFactory) const
build permutation operators
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)