35 BuildPermutation(
const Teuchos::RCP<Matrix> & A,
const Teuchos::RCP<const Map>& permRowMap,
38 const Scalar SC_ZERO = Teuchos::ScalarTraits<Scalar>::zero();
39 const Scalar SC_ONE = Teuchos::ScalarTraits<Scalar>::one();
41 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
42 const MT MT_ZERO = Teuchos::ScalarTraits<MT>::zero();
43 const MT MT_ONE = Teuchos::ScalarTraits<MT>::one();
45 const Teuchos::RCP< const Teuchos::Comm< int > > comm = A->getRowMap()->getComm();
46 int numProcs = comm->getSize();
47 int myRank = comm->getRank();
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();
55 std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > permutedDiagCandidates;
56 std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > keepDiagonalEntries;
57 std::vector<MT> Weights;
61 for (
size_t row = 0; row < A->getRowMap()->getLocalNumElements(); row++) {
64 if(permRowMap->isNodeGlobalElement(grow) ==
true)
continue;
66 size_t nnz = A->getNumEntriesInLocalRow(row);
69 Teuchos::ArrayView<const LocalOrdinal> indices;
70 Teuchos::ArrayView<const Scalar> vals;
71 A->getLocalRowView(row, indices, vals);
74 "MueLu::PermutationFactory::Build: number of nonzeros not equal to number of indices? Error.");
80 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
81 norm1 += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
82 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
83 maxVal = Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
84 gMaxValIdx = A->getColMap()->getGlobalElement(indices[j]);
88 if(grow == gMaxValIdx)
89 keepDiagonalEntries.push_back(std::make_pair(grow,grow));
94 for (
size_t row = 0; row < permRowMap->getLocalNumElements(); row++) {
96 LocalOrdinal lArow = A->getRowMap()->getLocalElement(grow);
97 size_t nnz = A->getNumEntriesInLocalRow(lArow);
100 Teuchos::ArrayView<const LocalOrdinal> indices;
101 Teuchos::ArrayView<const Scalar> vals;
102 A->getLocalRowView(lArow, indices, vals);
105 "MueLu::PermutationFactory::Build: number of nonzeros not equal to number of indices? Error.");
111 for (
size_t j = 0; j < Teuchos::as<size_t>(indices.size()); j++) {
112 norm1 += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
113 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > maxVal) {
114 maxVal = Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
115 gMaxValIdx = A->getColMap()->getGlobalElement(indices[j]);
119 if(maxVal > Teuchos::ScalarTraits<MT>::zero ()) {
120 permutedDiagCandidates.push_back(std::make_pair(grow,gMaxValIdx));
121 Weights.push_back(maxVal/(norm1*Teuchos::as<MT>(nnz)));
123 std::cout <<
"ATTENTION: row " << grow <<
" has only zero entries -> singular matrix!" << std::endl;
129 std::vector<int> permutation;
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);
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);
148 Teuchos::RCP<Export> exporter = ExportFactory::Build(gColVec->getMap(), gDomVec->getMap());
149 gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD);
150 gColVec->doImport(*gDomVec,*exporter,Xpetra::INSERT);
152 std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > permutedDiagCandidatesFiltered;
153 std::map<GlobalOrdinal, MT> gColId2Weight;
156 Teuchos::ArrayRCP< Scalar > ddata = gColVec->getDataNonConst(0);
157 for(
size_t i = 0; i < permutedDiagCandidates.size(); ++i) {
159 std::pair<GlobalOrdinal, GlobalOrdinal> pp = permutedDiagCandidates[permutation[i]];
163 LocalOrdinal lcol = A->getColMap()->getLocalElement(gcol);
164 if(Teuchos::ScalarTraits<Scalar>::real (ddata[lcol]) > MT_ZERO){
169 ddata[lcol] += SC_ONE;
171 permutedDiagCandidatesFiltered.push_back(std::make_pair(grow,gcol));
172 gColId2Weight[gcol] = Weights[permutation[i]];
177 gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD);
178 gColVec->doImport(*gDomVec,*exporter,Xpetra::INSERT);
185 std::vector<GlobalOrdinal> multipleColRequests;
189 std::queue<GlobalOrdinal> unusedColIdx;
191 for(
size_t sz = 0; sz < gDomVec->getLocalLength(); ++sz) {
192 Teuchos::ArrayRCP< const Scalar > arrDomVec = gDomVec->getData(0);
198 if(Teuchos::ScalarTraits<Scalar>::real(arrDomVec[sz]) > MT_ONE) {
199 multipleColRequests.push_back(gDomVec->getMap()->getGlobalElement(sz));
200 }
else if(Teuchos::ScalarTraits<Scalar>::real (arrDomVec[sz]) == MT_ZERO) {
201 unusedColIdx.push(gDomVec->getMap()->getGlobalElement(sz));
206 LocalOrdinal localMultColRequests = Teuchos::as<LocalOrdinal>(multipleColRequests.size());
216 if(globalMultColRequests > 0) {
222 std::vector<GlobalOrdinal> numMyMultColRequests(numProcs, 0);
223 std::vector<GlobalOrdinal> numGlobalMultColRequests(numProcs, 0);
224 numMyMultColRequests[myRank] = localMultColRequests;
225 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
226 numMyMultColRequests.data(),
227 numGlobalMultColRequests.data());
231 for (
int i=0; i<myRank-1; i++)
232 nMyOffset += numGlobalMultColRequests[i];
235 std::vector<GlobalOrdinal> procMultRequestedColIds(globalMultColRequests, zero);
236 std::vector<GlobalOrdinal> global_procMultRequestedColIds(globalMultColRequests, zero);
239 for(
size_t i = 0; i < multipleColRequests.size(); i++) {
240 procMultRequestedColIds[nMyOffset + i] = multipleColRequests[i];
244 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(globalMultColRequests),
245 procMultRequestedColIds.data(),
246 global_procMultRequestedColIds.data());
249 for (
size_t k = 0; k<global_procMultRequestedColIds.size(); k++) {
252 std::vector<MT> MyWeightForColId(numProcs, MT_ZERO);
253 std::vector<MT> GlobalWeightForColId(numProcs, MT_ZERO);
255 if(gColVec->getMap()->isNodeGlobalElement(globColId)) {
256 MyWeightForColId[myRank] = gColId2Weight[globColId];
258 MyWeightForColId[myRank] = MT_ZERO;
261 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
262 MyWeightForColId.data(),
263 GlobalWeightForColId.data());
265 if(gColVec->getMap()->isNodeGlobalElement(globColId)) {
268 MT winnerValue = MT_ZERO;
269 int winnerProcRank = 0;
270 for (
int proc = 0; proc < numProcs; proc++) {
271 if(GlobalWeightForColId[proc] > winnerValue) {
272 winnerValue = GlobalWeightForColId[proc];
273 winnerProcRank = proc;
280 if(myRank != winnerProcRank) {
282 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = permutedDiagCandidatesFiltered.begin();
283 while(p != permutedDiagCandidatesFiltered.end() )
285 if((*p).second == globColId)
286 p = permutedDiagCandidatesFiltered.erase(p);
298 std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> > RowColPairs;
299 RowColPairs.insert( RowColPairs.end(), keepDiagonalEntries.begin(), keepDiagonalEntries.end());
300 RowColPairs.insert( RowColPairs.end(), permutedDiagCandidatesFiltered.begin(), permutedDiagCandidatesFiltered.end());
305 gColVec->putScalar(SC_ZERO);
306 gDomVec->putScalar(SC_ZERO);
307 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator pl = RowColPairs.begin();
308 while(pl != RowColPairs.end() )
313 gColVec->sumIntoGlobalValue(jk,1.0);
316 gDomVec->doExport(*gColVec,*exporter,Xpetra::ADD);
317 for(
size_t sz = 0; sz < gDomVec->getLocalLength(); ++sz) {
318 Teuchos::ArrayRCP< const Scalar > arrDomVec = gDomVec->getData(0);
319 if(arrDomVec[sz] > 1.0) {
320 GetOStream(
Runtime0) <<
"RowColPairs has multiple column [" << sz <<
"]=" << arrDomVec[sz] << std::endl;
321 }
else if(arrDomVec[sz] == SC_ZERO) {
322 GetOStream(
Runtime0) <<
"RowColPairs has empty column [" << sz <<
"]=" << arrDomVec[sz] << std::endl;
355 Teuchos::RCP<Vector> Pperm = VectorFactory::Build(A->getRowMap());
356 Teuchos::RCP<Vector> Qperm = VectorFactory::Build(A->getDomainMap());
357 Teuchos::RCP<Vector> lQperm = VectorFactory::Build(A->getColMap());
359 Pperm->putScalar(SC_ZERO);
360 Qperm->putScalar(SC_ZERO);
361 lQperm->putScalar(SC_ZERO);
364 Teuchos::RCP<Export> QpermExporter = ExportFactory::Build(lQperm->getMap(), Qperm->getMap());
366 Teuchos::RCP<Vector> RowIdStatus = VectorFactory::Build(A->getRowMap());
367 Teuchos::RCP<Vector> ColIdStatus = VectorFactory::Build(A->getDomainMap());
368 Teuchos::RCP<Vector> lColIdStatus = VectorFactory::Build(A->getColMap());
369 Teuchos::RCP<Vector> ColIdUsed = VectorFactory::Build(A->getDomainMap());
370 RowIdStatus->putScalar(SC_ZERO);
371 ColIdStatus->putScalar(SC_ZERO);
372 lColIdStatus->putScalar(SC_ZERO);
373 ColIdUsed->putScalar(SC_ZERO);
386 Teuchos::ArrayRCP< Scalar > RowIdStatusArray = RowIdStatus->getDataNonConst(0);
387 Teuchos::ArrayRCP< Scalar > lColIdStatusArray = lColIdStatus->getDataNonConst(0);
389 typename std::vector<std::pair<GlobalOrdinal, GlobalOrdinal> >::iterator p = RowColPairs.begin();
390 while(p != RowColPairs.end() )
398 if(RowIdStatusArray[lik] == SC_ZERO) {
399 RowIdStatusArray[lik] = SC_ONE;
400 lColIdStatusArray[ljk] = SC_ONE;
401 Pperm->replaceLocalValue(lik, ik);
402 lQperm->replaceLocalValue(ljk, ik);
403 ColIdUsed->replaceGlobalValue(ik,SC_ONE);
404 p = RowColPairs.erase(p);
407 if(floor(ik/nDofsPerNode) != floor(jk/nDofsPerNode)) {
408 lWideRangeColPermutations++;
417 Qperm->doExport(*lQperm, *QpermExporter, Xpetra::ABSMAX);
418 ColIdStatus->doExport(*lColIdStatus, *QpermExporter, Xpetra::ABSMAX);
421 if(RowColPairs.size()>0) GetOStream(
Warnings0) <<
"MueLu::PermutationFactory: There are Row/Col pairs left!!!" << std::endl;
426 size_t cntFreeRowIdx = 0;
427 std::queue<GlobalOrdinal> qFreeGRowIdx;
428 for (
size_t lik = 0; lik < RowIdStatus->getLocalLength(); ++lik) {
429 if(RowIdStatusArray[lik] == SC_ZERO) {
431 qFreeGRowIdx.push(RowIdStatus->getMap()->getGlobalElement(lik));
436 for (
size_t lik = 0; lik < RowIdStatus->getLocalLength(); ++lik) {
437 if(RowIdStatusArray[lik] == SC_ZERO) {
438 RowIdStatusArray[lik] = SC_ONE;
439 Pperm->replaceLocalValue(lik, qFreeGRowIdx.front());
441 if(floor(qFreeGRowIdx.front()/nDofsPerNode) != floor(RowIdStatus->getMap()->getGlobalElement(lik)/nDofsPerNode)) {
442 lWideRangeRowPermutations++;
449 size_t cntUnusedColIdx = 0;
450 std::queue<GlobalOrdinal> qUnusedGColIdx;
451 size_t cntFreeColIdx = 0;
452 std::queue<GlobalOrdinal> qFreeGColIdx;
454 Teuchos::ArrayRCP< Scalar > ColIdStatusArray = ColIdStatus->getDataNonConst(0);
455 Teuchos::ArrayRCP< Scalar > ColIdUsedArray = ColIdUsed->getDataNonConst(0);
458 for (
size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
459 if(ColIdStatusArray[ljk] == SC_ZERO) {
461 qFreeGColIdx.push(ColIdStatus->getMap()->getGlobalElement(ljk));
465 for (
size_t ljk = 0; ljk < ColIdUsed->getLocalLength(); ++ljk) {
466 if(ColIdUsedArray[ljk] == SC_ZERO) {
468 qUnusedGColIdx.push(ColIdUsed->getMap()->getGlobalElement(ljk));
474 for (
size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
476 if(cntUnusedColIdx == 0)
break;
478 if(ColIdStatusArray[ljk] == SC_ZERO) {
479 ColIdStatusArray[ljk] = SC_ONE;
480 Qperm->replaceLocalValue(ljk, qUnusedGColIdx.front());
481 ColIdUsed->replaceGlobalValue(qUnusedGColIdx.front(),SC_ONE);
483 if(floor(qUnusedGColIdx.front()/nDofsPerNode) != floor(ColIdStatus->getMap()->getGlobalElement(ljk)/nDofsPerNode)) {
484 lWideRangeColPermutations++;
486 qUnusedGColIdx.pop();
499 for (
size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
500 if(ColIdStatusArray[ljk] == SC_ZERO) {
508 MueLu_sumAll(comm, Teuchos::as<GlobalOrdinal>(local_cntFreeColIdx), global_cntFreeColIdx);
510 std::cout <<
"global # of empty column idx entries in Qperm: " << global_cntFreeColIdx << std::endl;
514 if(global_cntFreeColIdx > 0) {
519 MueLu_sumAll(comm, Teuchos::as<GlobalOrdinal>(local_cntUnusedColIdx), global_cntUnusedColIdx);
521 std::cout <<
"global # of unused column idx: " << global_cntUnusedColIdx << std::endl;
525 std::vector<LocalOrdinal> local_UnusedColIdxOnProc (numProcs);
526 std::vector<LocalOrdinal> global_UnusedColIdxOnProc(numProcs);
527 local_UnusedColIdxOnProc[myRank] = local_cntUnusedColIdx;
528 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
529 local_UnusedColIdxOnProc.data(),
530 global_UnusedColIdxOnProc.data());
533 std::cout <<
"PROC " << myRank <<
" global num unused indices per proc: ";
534 for (
size_t ljk = 0; ljk < global_UnusedColIdxOnProc.size(); ++ljk) {
535 std::cout <<
" " << global_UnusedColIdxOnProc[ljk];
537 std::cout << std::endl;
541 std::vector<GlobalOrdinal> local_UnusedColIdxVector(Teuchos::as<size_t>(global_cntUnusedColIdx));
542 std::vector<GlobalOrdinal> global_UnusedColIdxVector(Teuchos::as<size_t>(global_cntUnusedColIdx));
544 for(
int proc=0; proc<myRank; proc++) {
545 global_cntUnusedColIdxStartIter += global_UnusedColIdxOnProc[proc];
547 for(
GlobalOrdinal k = global_cntUnusedColIdxStartIter; k < global_cntUnusedColIdxStartIter+local_cntUnusedColIdx; k++) {
548 local_UnusedColIdxVector[k] = qUnusedGColIdx.front();
549 qUnusedGColIdx.pop();
551 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, Teuchos::as<int>(global_cntUnusedColIdx),
552 local_UnusedColIdxVector.data(),
553 global_UnusedColIdxVector.data());
555 std::cout <<
"PROC " << myRank <<
" global UnusedGColIdx: ";
556 for (
size_t ljk = 0; ljk < global_UnusedColIdxVector.size(); ++ljk) {
557 std::cout <<
" " << global_UnusedColIdxVector[ljk];
559 std::cout << std::endl;
566 std::vector<LocalOrdinal> local_EmptyColIdxOnProc (numProcs);
567 std::vector<LocalOrdinal> global_EmptyColIdxOnProc(numProcs);
568 local_EmptyColIdxOnProc[myRank] = local_cntFreeColIdx;
569 Teuchos::reduceAll(*comm, Teuchos::REDUCE_MAX, numProcs,
570 local_EmptyColIdxOnProc.data(),
571 global_EmptyColIdxOnProc.data());
574 std::cout <<
"PROC " << myRank <<
" global num of needed column indices: ";
575 for (
size_t ljk = 0; ljk < global_EmptyColIdxOnProc.size(); ++ljk) {
576 std::cout <<
" " << global_EmptyColIdxOnProc[ljk];
578 std::cout << std::endl;
584 for(
int proc=0; proc<myRank; proc++) {
585 global_UnusedColStartIdx += global_EmptyColIdxOnProc[proc];
589 GetOStream(
Statistics0) <<
"PROC " << myRank <<
" is allowd to use the following column gids: ";
590 for(
GlobalOrdinal k = global_UnusedColStartIdx; k < global_UnusedColStartIdx + Teuchos::as<GlobalOrdinal>(cntFreeColIdx); k++) {
591 GetOStream(
Statistics0) << global_UnusedColIdxVector[k] <<
" ";
598 Teuchos::ArrayRCP< Scalar > ColIdStatusArray = ColIdStatus->getDataNonConst(0);
600 for (
size_t ljk = 0; ljk < ColIdStatus->getLocalLength(); ++ljk) {
602 if(ColIdStatusArray[ljk] == SC_ZERO) {
603 ColIdStatusArray[ljk] = SC_ONE;
604 Qperm->replaceLocalValue(ljk, global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter]);
605 ColIdUsed->replaceGlobalValue(global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter],SC_ONE);
607 if(floor(global_UnusedColIdxVector[global_UnusedColStartIdx + array_iter]/nDofsPerNode) != floor(ColIdStatus->getMap()->getGlobalElement(ljk)/nDofsPerNode)) {
608 lWideRangeColPermutations++;
620 Teuchos::RCP<CrsMatrixWrap> permPTmatrix = Teuchos::rcp(
new CrsMatrixWrap(A->getRowMap(),1));
621 Teuchos::RCP<CrsMatrixWrap> permQTmatrix = Teuchos::rcp(
new CrsMatrixWrap(A->getRowMap(),1));
624 Teuchos::ArrayRCP< Scalar > PpermData = Pperm->getDataNonConst(0);
625 Teuchos::ArrayRCP< Scalar > QpermData = Qperm->getDataNonConst(0);
627 for(
size_t row=0; row<A->getLocalNumRows(); row++) {
631 Teuchos::ArrayRCP<GlobalOrdinal> indoutP(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(PpermData[row])));
632 Teuchos::ArrayRCP<GlobalOrdinal> indoutQ(1,Teuchos::as<GO>(Teuchos::ScalarTraits<Scalar>::real(QpermData[row])));
633 Teuchos::ArrayRCP<Scalar> valout(1,SC_ONE);
634 permPTmatrix->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indoutP.view(0,indoutP.size()), valout.view(0,valout.size()));
635 permQTmatrix->insertGlobalValues (A->getRowMap()->getGlobalElement(row), indoutQ.view(0,indoutQ.size()), valout.view(0,valout.size()));
639 permPTmatrix->fillComplete();
640 permQTmatrix->fillComplete();
644 for(
size_t row=0; row<permPTmatrix->getLocalNumRows(); row++) {
645 if(permPTmatrix->getNumEntriesInLocalRow(row) != 1)
646 GetOStream(
Warnings0) <<
"#entries in row " << row <<
" of permPTmatrix is " << permPTmatrix->getNumEntriesInLocalRow(row) << std::endl;
647 if(permPmatrix->getNumEntriesInLocalRow(row) != 1)
648 GetOStream(
Warnings0) <<
"#entries in row " << row <<
" of permPmatrix is " << permPmatrix->getNumEntriesInLocalRow(row) << std::endl;
649 if(permQTmatrix->getNumEntriesInLocalRow(row) != 1)
650 GetOStream(
Warnings0) <<
"#entries in row " << row <<
" of permQmatrix is " << permQTmatrix->getNumEntriesInLocalRow(row) << std::endl;
654 Teuchos::RCP<Matrix> ApermQt = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *permQTmatrix,
false, GetOStream(
Statistics2));
655 Teuchos::RCP<Matrix> permPApermQt = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*permPmatrix,
false, *ApermQt,
false, GetOStream(
Statistics2));
664 Teuchos::RCP<Vector> diagVec = VectorFactory::Build(permPApermQt->getRowMap(),
true);
665 Teuchos::RCP<Vector> invDiagVec = VectorFactory::Build(permPApermQt->getRowMap(),
true);
666 Teuchos::RCP<CrsMatrixWrap> diagScalingOp = Teuchos::rcp(
new CrsMatrixWrap(permPApermQt->getRowMap(),1));
668 permPApermQt->getLocalDiagCopy(*diagVec);
669 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
670 Teuchos::ArrayRCP< Scalar > invDiagVecData = invDiagVec->getDataNonConst(0);
671 for(
size_t i = 0; i<diagVec->getMap()->getLocalNumElements(); ++i) {
672 if(diagVecData[i] != SC_ZERO)
673 invDiagVecData[i] = SC_ONE / diagVecData[i];
675 invDiagVecData[i] = SC_ONE;
676 GetOStream(
Statistics0) <<
"MueLu::PermutationFactory: found zero on diagonal in row " << i << std::endl;
680 for(
size_t row=0; row<A->getLocalNumRows(); row++) {
681 Teuchos::ArrayRCP<GlobalOrdinal> indout(1,permPApermQt->getRowMap()->getGlobalElement(row));
682 Teuchos::ArrayRCP<Scalar> valout(1,invDiagVecData[row]);
683 diagScalingOp->insertGlobalValues(A->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
686 diagScalingOp->fillComplete();
688 Teuchos::RCP<Matrix> scaledA = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*diagScalingOp,
false, *permPApermQt,
false, GetOStream(
Statistics2));
689 currentLevel.
Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(scaledA), genFactory);
691 currentLevel.
Set(
"permA", Teuchos::rcp_dynamic_cast<Matrix>(permPApermQt), genFactory);
692 currentLevel.
Set(
"permP", Teuchos::rcp_dynamic_cast<Matrix>(permPmatrix), genFactory);
693 currentLevel.
Set(
"permQT", Teuchos::rcp_dynamic_cast<Matrix>(permQTmatrix), genFactory);
694 currentLevel.
Set(
"permScaling", Teuchos::rcp_dynamic_cast<Matrix>(diagScalingOp), genFactory);
698 Teuchos::RCP<Vector> diagPVec = VectorFactory::Build(permPmatrix->getRowMap(),
true);
699 permPmatrix->getLocalDiagCopy(*diagPVec);
703 Teuchos::ArrayRCP< const Scalar > diagPVecData = diagPVec->getData(0);
704 for(
size_t i = 0; i<diagPVec->getMap()->getLocalNumElements(); ++i) {
705 if(diagPVecData[i] == SC_ZERO) {
706 lNumRowPermutations++;
712 MueLu_sumAll(diagPVec->getMap()->getComm(), Teuchos::as<GlobalOrdinal>(lNumRowPermutations), gNumRowPermutations);
716 Teuchos::RCP<Vector> diagQTVec = VectorFactory::Build(permQTmatrix->getRowMap(),
true);
717 permQTmatrix->getLocalDiagCopy(*diagQTVec);
721 Teuchos::ArrayRCP< const Scalar > diagQTVecData = diagQTVec->getData(0);
722 for(
size_t i = 0; i<diagQTVec->getMap()->getLocalNumElements(); ++i) {
723 if(diagQTVecData[i] == SC_ZERO) {
724 lNumColPermutations++;
730 MueLu_sumAll(diagQTVec->getMap()->getComm(), Teuchos::as<GlobalOrdinal>(lNumColPermutations), gNumColPermutations);
732 currentLevel.
Set(
"#RowPermutations", gNumRowPermutations, genFactory);
733 currentLevel.
Set(
"#ColPermutations", gNumColPermutations, genFactory);
734 currentLevel.
Set(
"#WideRangeRowPermutations", gWideRangeRowPermutations, genFactory);
735 currentLevel.
Set(
"#WideRangeColPermutations", gWideRangeColPermutations, genFactory);
737 GetOStream(
Statistics0) <<
"#Row permutations/max possible permutations: " << gNumRowPermutations <<
"/" << diagPVec->getMap()->getGlobalNumElements() << std::endl;
738 GetOStream(
Statistics0) <<
"#Column permutations/max possible permutations: " << gNumColPermutations <<
"/" << diagQTVec->getMap()->getGlobalNumElements() << std::endl;
739 GetOStream(
Runtime1) <<
"#wide range row permutations: " << gWideRangeRowPermutations <<
" #wide range column permutations: " << gWideRangeColPermutations << std::endl;