46 #ifndef MUELU_UTILITIESBASE_DECL_HPP 47 #define MUELU_UTILITIESBASE_DECL_HPP 57 #include <Teuchos_DefaultComm.hpp> 58 #include <Teuchos_ScalarTraits.hpp> 59 #include <Teuchos_ParameterList.hpp> 61 #include <Xpetra_BlockedCrsMatrix_fwd.hpp> 62 #include <Xpetra_CrsMatrix_fwd.hpp> 63 #include <Xpetra_CrsMatrixWrap_fwd.hpp> 64 #include <Xpetra_Map_fwd.hpp> 65 #include <Xpetra_BlockedMap_fwd.hpp> 66 #include <Xpetra_MapFactory_fwd.hpp> 67 #include <Xpetra_Matrix_fwd.hpp> 68 #include <Xpetra_MatrixFactory_fwd.hpp> 69 #include <Xpetra_MultiVector_fwd.hpp> 70 #include <Xpetra_MultiVectorFactory_fwd.hpp> 71 #include <Xpetra_Operator_fwd.hpp> 72 #include <Xpetra_Vector_fwd.hpp> 73 #include <Xpetra_BlockedMultiVector.hpp> 74 #include <Xpetra_BlockedVector.hpp> 75 #include <Xpetra_VectorFactory_fwd.hpp> 76 #include <Xpetra_ExportFactory.hpp> 78 #include <Xpetra_Import.hpp> 79 #include <Xpetra_ImportFactory.hpp> 80 #include <Xpetra_MatrixMatrix.hpp> 81 #include <Xpetra_CrsMatrixWrap.hpp> 82 #include <Xpetra_StridedMap.hpp> 90 #define MueLu_sumAll(rcpComm, in, out) \ 91 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out)) 92 #define MueLu_minAll(rcpComm, in, out) \ 93 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out)) 94 #define MueLu_maxAll(rcpComm, in, out) \ 95 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out)) 108 class UtilitiesBase {
110 #undef MUELU_UTILITIESBASE_SHORT 113 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>
CrsMatrixWrap;
114 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
CrsMatrix;
115 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
Matrix;
116 typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
Vector;
117 typedef Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
BlockedVector;
118 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
MultiVector;
120 typedef Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>
BlockedMap;
121 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>
Map;
123 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Magnitude;
126 static RCP<Matrix>
Crs2Op(RCP<CrsMatrix> Op) {
128 return Teuchos::null;
139 size_t numRows = A.getRowMap()->getNodeNumElements();
140 Teuchos::ArrayRCP<Scalar> diag(numRows);
141 Teuchos::ArrayView<const LocalOrdinal> cols;
142 Teuchos::ArrayView<const Scalar> vals;
143 for (
size_t i = 0; i < numRows; ++i) {
144 A.getLocalRowView(i, cols, vals);
146 for (; j < cols.size(); ++j) {
147 if (Teuchos::as<size_t>(cols[j]) == i) {
152 if (j == cols.size()) {
154 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
167 Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer(
"UtilitiesBase::GetMatrixDiagonalInverse");
169 RCP<const Map> rowMap = A.getRowMap();
170 RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
172 A.getLocalDiagCopy(*diag);
186 Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100,
187 Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(),
188 const bool replaceSingleEntryRowWithZero =
false) {
190 typedef Teuchos::ScalarTraits<Scalar> TST;
192 RCP<Vector> diag = Teuchos::null;
193 const Scalar zero = TST::zero();
194 const Scalar one = TST::one();
195 const Scalar two = one + one;
199 Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
201 RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
202 Teuchos::rcp_dynamic_cast<
const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
203 if(bA == Teuchos::null) {
204 RCP<const Map> rowMap = rcpA->getRowMap();
205 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
206 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
207 Teuchos::Array<Scalar> regSum(diag->getLocalLength());
208 Teuchos::ArrayView<const LocalOrdinal> cols;
209 Teuchos::ArrayView<const Scalar> vals;
211 std::vector<int> nnzPerRow(rowMap->getNodeNumElements());
213 const Magnitude zeroMagn = TST::magnitude(zero);
214 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
216 rcpA->getLocalRowView(i, cols, vals);
219 regSum[i] += vals[j];
220 const Magnitude rowEntryMagn = TST::magnitude(vals[j]);
221 if (rowEntryMagn > zeroMagn)
223 diagVals[i] += rowEntryMagn;
227 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
228 if (replaceSingleEntryRowWithZero && nnzPerRow[i] <= static_cast<int>(1))
230 else if (replaceSingleEntryRowWithZero && diagVals[i] != zero && TST::magnitude(diagVals[i]) < TST::magnitude(two*regSum[i]))
231 diagVals[i] = one / (two*regSum[i]);
233 if(TST::magnitude(diagVals[i]) > tol)
234 diagVals[i] = one / diagVals[i];
236 diagVals[i] = valReplacement;
242 TEUCHOS_TEST_FOR_EXCEPTION(doReciprocal, Xpetra::Exceptions::RuntimeError,
243 "UtilitiesBase::GetLumpedMatrixDiagonal(): extracting reciprocal of diagonal of a blocked matrix is not supported");
244 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),
true);
246 for (
size_t row = 0; row < bA->Rows(); ++row) {
247 for (
size_t col = 0; col < bA->Cols(); ++col) {
248 if (!bA->getMatrix(row,col).is_null()) {
250 bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA->getMatrix(row,col)) == Teuchos::null);
251 RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag,row,bThyraMode);
253 ddtemp->update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
254 bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
271 size_t numRows = A.getRowMap()->getNodeNumElements();
272 Magnitude ZERO = Teuchos::ScalarTraits<Magnitude>::zero();
273 Teuchos::ArrayRCP<Magnitude> maxvec(numRows);
274 Teuchos::ArrayView<const LocalOrdinal> cols;
275 Teuchos::ArrayView<const Scalar> vals;
276 for (
size_t i = 0; i < numRows; ++i) {
277 A.getLocalRowView(i, cols, vals);
280 if (Teuchos::as<size_t>(cols[j]) != i) {
281 mymax = std::max(mymax,-Teuchos::ScalarTraits<Scalar>::real(vals[j]));
289 static Teuchos::ArrayRCP<Magnitude>
GetMatrixMaxMinusOffDiagonal(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> &BlockNumber) {
290 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()),std::runtime_error,
"GetMatrixMaxMinusOffDiagonal: BlockNumber must match's A's column map.");
292 Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
294 size_t numRows = A.getRowMap()->getNodeNumElements();
295 Magnitude ZERO = Teuchos::ScalarTraits<Magnitude>::zero();
296 Teuchos::ArrayRCP<Magnitude> maxvec(numRows);
297 Teuchos::ArrayView<const LocalOrdinal> cols;
298 Teuchos::ArrayView<const Scalar> vals;
299 for (
size_t i = 0; i < numRows; ++i) {
300 A.getLocalRowView(i, cols, vals);
303 if (Teuchos::as<size_t>(cols[j]) != i && block_id[i] == block_id[cols[j]]) {
304 mymax = std::max(mymax,-Teuchos::ScalarTraits<Scalar>::real(vals[j]));
321 static Teuchos::RCP<Vector>
GetInverse(Teuchos::RCP<const Vector> v,
Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100,
Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
323 RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),
true);
326 RCP<const BlockedVector> bv = Teuchos::rcp_dynamic_cast<
const BlockedVector>(v);
327 if(bv.is_null() ==
false) {
328 RCP<BlockedVector> bret = Teuchos::rcp_dynamic_cast<
BlockedVector>(ret);
329 TEUCHOS_TEST_FOR_EXCEPTION(bret.is_null() ==
true,
MueLu::Exceptions::RuntimeError,
"MueLu::UtilitiesBase::GetInverse: return vector should be of type BlockedVector");
330 RCP<const BlockedMap> bmap = bv->getBlockedMap();
331 for(
size_t r = 0; r < bmap->getNumMaps(); ++r) {
332 RCP<const MultiVector> submvec = bv->getMultiVector(r,bmap->getThyraMode());
333 RCP<const Vector> subvec = submvec->getVector(0);
335 bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
341 ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
342 ArrayRCP<const Scalar> inputVals = v->getData(0);
343 for (
size_t i = 0; i < v->getMap()->getNodeNumElements(); ++i) {
344 if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
345 retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
347 retVals[i] = valReplacement;
360 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
363 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(rowMap);
364 if(!browMap.is_null()) rowMap = browMap->getMap();
366 RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
372 Teuchos::ArrayRCP<size_t> offsets;
373 crsOp->getLocalDiagOffsets(offsets);
374 crsOp->getLocalDiagCopy(*localDiag,offsets());
377 ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
380 localDiagVals[i] = diagVals[i];
381 localDiagVals = diagVals = null;
384 RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
385 RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
386 importer = A.getCrsGraph()->getImporter();
387 if (importer == Teuchos::null) {
388 importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
390 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
406 using STS =
typename Teuchos::ScalarTraits<SC>;
409 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
410 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(rowMap);
411 if(!browMap.is_null()) rowMap = browMap->getMap();
413 RCP<Vector> local = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(rowMap);
414 RCP<Vector> ghosted = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(colMap,
true);
415 ArrayRCP<SC> localVals = local->getDataNonConst(0);
417 for (LO row = 0; row < static_cast<LO>(A.getRowMap()->getNodeNumElements()); ++row) {
418 size_t nnz = A.getNumEntriesInLocalRow(row);
419 ArrayView<const LO> indices;
420 ArrayView<const SC> vals;
421 A.getLocalRowView(row, indices, vals);
425 for (LO colID = 0; colID < static_cast<LO>(nnz); colID++) {
426 if(indices[colID] != row) {
433 RCP< const Xpetra::Import<LO,GO,Node> > importer;
434 importer = A.getCrsGraph()->getImporter();
435 if (importer == Teuchos::null) {
436 importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
438 ghosted->doImport(*local, *(importer), Xpetra::INSERT);
444 static RCP<Xpetra::Vector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> >
446 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
447 using STS =
typename Teuchos::ScalarTraits<Scalar>;
448 using MTS =
typename Teuchos::ScalarTraits<Magnitude>;
453 using RealValuedVector = Xpetra::Vector<MT,LO,GO,Node>;
456 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(rowMap);
457 if(!browMap.is_null()) rowMap = browMap->getMap();
459 RCP<RealValuedVector> local = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(rowMap);
460 RCP<RealValuedVector> ghosted = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(colMap,
true);
461 ArrayRCP<MT> localVals = local->getDataNonConst(0);
463 for (LO rowIdx = 0; rowIdx < static_cast<LO>(A.getRowMap()->getNodeNumElements()); ++rowIdx) {
464 size_t nnz = A.getNumEntriesInLocalRow(rowIdx);
465 ArrayView<const LO> indices;
466 ArrayView<const SC> vals;
467 A.getLocalRowView(rowIdx, indices, vals);
471 for (LO colID = 0; colID < static_cast<LO>(nnz); ++colID) {
472 if(indices[colID] != rowIdx) {
473 si += STS::magnitude(vals[colID]);
476 localVals[rowIdx] = si;
479 RCP< const Xpetra::Import<LO,GO,Node> > importer;
480 importer = A.getCrsGraph()->getImporter();
481 if (importer == Teuchos::null) {
482 importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
484 ghosted->doImport(*local, *(importer), Xpetra::INSERT);
495 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
496 const size_t numVecs = X.getNumVectors();
497 RCP<MultiVector> RES =
Residual(Op, X, RHS);
498 Teuchos::Array<Magnitude> norms(numVecs);
504 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
505 const size_t numVecs = X.getNumVectors();
507 Teuchos::Array<Magnitude> norms(numVecs);
513 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
514 const size_t numVecs = X.getNumVectors();
516 RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(RHS.getMap(), numVecs,
false);
517 Op.residual(X,RHS,*RES);
523 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides");
524 TEUCHOS_TEST_FOR_EXCEPTION(Resid.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of residual vectors != number of right-hand sides");
525 Op.residual(X,RHS,Resid);
541 LocalOrdinal niters = 10,
Magnitude tolerance = 1e-2,
bool verbose =
false,
unsigned int seed = 123) {
543 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
546 RCP<Vector> diagInvVec;
548 RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
549 A.getLocalDiagCopy(*diagVec);
550 diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
551 diagInvVec->reciprocal(*diagVec);
570 LocalOrdinal niters = 10,
Magnitude tolerance = 1e-2,
bool verbose =
false,
unsigned int seed = 123) {
572 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
575 RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
576 RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
577 RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
582 Teuchos::Array<Magnitude> norms(1);
584 typedef Teuchos::ScalarTraits<Scalar> STS;
586 const Scalar zero = STS::zero(), one = STS::one();
589 Magnitude residual = STS::magnitude(zero);
592 for (
int iter = 0; iter < niters; ++iter) {
594 q->update(one/norms[0], *z, zero);
596 if (diagInvVec != Teuchos::null)
597 z->elementWiseMultiply(one, *diagInvVec, *z, zero);
600 if (iter % 100 == 0 || iter + 1 == niters) {
601 r->update(1.0, *z, -lambda, *q, zero);
603 residual = STS::magnitude(norms[0] / lambda);
605 std::cout <<
"Iter = " << iter
606 <<
" Lambda = " << lambda
607 <<
" Residual of A*q - lambda*q = " << residual
617 static RCP<Teuchos::FancyOStream>
MakeFancy(std::ostream& os) {
618 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
627 const size_t numVectors = v.size();
629 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
630 for (
size_t j = 0; j < numVectors; j++) {
631 d += (v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
633 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
640 static typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Distance2(
const Teuchos::ArrayView<double> & weight,
const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>> &v,
LocalOrdinal i0,
LocalOrdinal i1) {
641 const size_t numVectors = v.size();
643 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
644 for (
size_t j = 0; j < numVectors; j++) {
645 d += weight[j]*(v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
647 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
663 static Teuchos::ArrayRCP<const bool>
DetectDirichletRows(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero(),
bool count_twos_as_dirichlet=
false) {
665 typedef Teuchos::ScalarTraits<Scalar> STS;
666 ArrayRCP<bool> boundaryNodes(numRows,
true);
667 if (count_twos_as_dirichlet) {
669 ArrayView<const LocalOrdinal> indices;
670 ArrayView<const Scalar> vals;
671 A.getLocalRowView(row, indices, vals);
672 size_t nnz = A.getNumEntriesInLocalRow(row);
675 for (col = 0; col < nnz; col++)
676 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
677 if (!boundaryNodes[row])
679 boundaryNodes[row] =
false;
682 boundaryNodes[row] =
true;
687 ArrayView<const LocalOrdinal> indices;
688 ArrayView<const Scalar> vals;
689 A.getLocalRowView(row, indices, vals);
690 size_t nnz = A.getNumEntriesInLocalRow(row);
692 for (
size_t col = 0; col < nnz; col++)
693 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
694 boundaryNodes[row] =
false;
699 return boundaryNodes;
715 static Teuchos::ArrayRCP<const bool>
DetectDirichletRowsExt(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
bool & bHasZeroDiagonal,
const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
718 bHasZeroDiagonal =
false;
720 Teuchos::RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
721 A.getLocalDiagCopy(*diagVec);
722 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
725 typedef Teuchos::ScalarTraits<Scalar> STS;
726 ArrayRCP<bool> boundaryNodes(numRows,
false);
728 ArrayView<const LocalOrdinal> indices;
729 ArrayView<const Scalar> vals;
730 A.getLocalRowView(row, indices, vals);
732 bool bHasDiag =
false;
733 for (decltype(indices.size()) col = 0; col < indices.size(); col++) {
734 if ( indices[col] != row) {
735 if (STS::magnitude(vals[col] / STS::magnitude(sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col]))) ) > tol) {
738 }
else bHasDiag =
true;
740 if (bHasDiag ==
false) bHasZeroDiagonal =
true;
741 else if(nnz == 0) boundaryNodes[row] =
true;
743 return boundaryNodes;
754 Teuchos::ArrayRCP<bool> nonzeros) {
755 TEUCHOS_ASSERT(vals.size() == nonzeros.size());
756 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
757 const magnitudeType eps = 2.0*Teuchos::ScalarTraits<magnitudeType>::eps();
758 for(
size_t i=0; i<static_cast<size_t>(vals.size()); i++) {
759 nonzeros[i] = (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > eps);
772 const Teuchos::ArrayRCP<bool>& dirichletRows,
773 Teuchos::ArrayRCP<bool> dirichletCols,
774 Teuchos::ArrayRCP<bool> dirichletDomain) {
775 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
776 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A .getDomainMap();
777 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
778 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
779 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == rowMap->getNodeNumElements());
780 TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == colMap->getNodeNumElements());
781 TEUCHOS_ASSERT(static_cast<size_t>(dirichletDomain.size()) == domMap->getNodeNumElements());
782 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, 1,
true);
784 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
785 if (dirichletRows[i]) {
786 ArrayView<const LocalOrdinal> indices;
787 ArrayView<const Scalar> values;
788 A.getLocalRowView(i,indices,values);
789 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
790 myColsToZero->replaceLocalValue(indices[j],0,one);
794 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
795 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
796 if (!importer.is_null()) {
797 globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, 1,
true);
799 globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
801 myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
804 globalColsToZero = myColsToZero;
806 FindNonZeros(globalColsToZero->getData(0),dirichletDomain);
823 static void ApplyRowSumCriterion(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
824 typedef Teuchos::ScalarTraits<Scalar> STS;
825 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
826 typedef Teuchos::ScalarTraits<MT> MTS;
827 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> rowmap = A.getRowMap();
828 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getNodeNumElements()); ++row) {
829 size_t nnz = A.getNumEntriesInLocalRow(row);
830 ArrayView<const LocalOrdinal> indices;
831 ArrayView<const Scalar> vals;
832 A.getLocalRowView(row, indices, vals);
834 Scalar rowsum = STS::zero();
835 Scalar diagval = STS::zero();
837 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
840 diagval = vals[colID];
841 rowsum += vals[colID];
844 if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
846 dirichletRows[row] =
true;
851 static void ApplyRowSumCriterion(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> &BlockNumber,
const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
852 typedef Teuchos::ScalarTraits<Scalar> STS;
853 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
854 typedef Teuchos::ScalarTraits<MT> MTS;
855 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowmap = A.getRowMap();
857 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()),std::runtime_error,
"ApplyRowSumCriterion: BlockNumber must match's A's column map.");
859 Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
860 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getNodeNumElements()); ++row) {
861 size_t nnz = A.getNumEntriesInLocalRow(row);
862 ArrayView<const LocalOrdinal> indices;
863 ArrayView<const Scalar> vals;
864 A.getLocalRowView(row, indices, vals);
866 Scalar rowsum = STS::zero();
867 Scalar diagval = STS::zero();
868 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
871 diagval = vals[colID];
872 if(block_id[row] == block_id[col])
873 rowsum += vals[colID];
877 if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
879 dirichletRows[row] =
true;
896 static Teuchos::ArrayRCP<const bool>
DetectDirichletCols(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
897 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
898 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
899 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
900 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
901 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
902 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,1);
903 myColsToZero->putScalar(zero);
905 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
906 if (dirichletRows[i]) {
907 Teuchos::ArrayView<const LocalOrdinal> indices;
908 Teuchos::ArrayView<const Scalar> values;
909 A.getLocalRowView(i,indices,values);
910 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
911 myColsToZero->replaceLocalValue(indices[j],0,one);
915 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,1);
916 globalColsToZero->putScalar(zero);
917 Teuchos::RCP<Xpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
919 globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
921 myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
922 Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
923 Teuchos::ArrayRCP<bool> dirichletCols(colMap->getNodeNumElements(),
true);
924 Magnitude eps = Teuchos::ScalarTraits<Magnitude>::eps();
925 for(
size_t i=0; i<colMap->getNodeNumElements(); i++) {
926 dirichletCols[i] = Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i])>2.0*eps;
928 return dirichletCols;
936 static Scalar Frobenius(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
941 TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()),
Exceptions::Incompatible,
"MueLu::CGSolver::Frobenius: row maps are incompatible");
942 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(),
Exceptions::RuntimeError,
"Matrices must be fill completed");
944 const Map& AColMap = *A.getColMap();
945 const Map& BColMap = *B.getColMap();
947 Teuchos::ArrayView<const LocalOrdinal> indA, indB;
948 Teuchos::ArrayView<const Scalar> valA, valB;
949 size_t nnzA = 0, nnzB = 0;
961 Teuchos::Array<Scalar> valBAll(BColMap.getNodeNumElements());
963 LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
964 Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
965 size_t numRows = A.getNodeNumRows();
966 for (
size_t i = 0; i < numRows; i++) {
967 A.getLocalRowView(i, indA, valA);
968 B.getLocalRowView(i, indB, valB);
973 for (
size_t j = 0; j < nnzB; j++)
974 valBAll[indB[j]] = valB[j];
976 for (
size_t j = 0; j < nnzA; j++) {
979 LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
981 f += valBAll[ind] * valA[j];
985 for (
size_t j = 0; j < nnzB; j++)
986 valBAll[indB[j]] = zero;
1009 int maxint = INT_MAX;
1010 int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.getRank()+1)/(comm.getSize()+one)) );
1011 if (mySeed < 1 || mySeed == maxint) {
1012 std::ostringstream errStr;
1013 errStr <<
"Error detected with random seed = " << mySeed <<
". It should be in the interval [1,2^31-2].";
1018 Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
1030 static void FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
1031 std::vector<LocalOrdinal>& dirichletRows,
bool count_twos_as_dirichlet=
false) {
1032 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
1033 dirichletRows.resize(0);
1034 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
1035 Teuchos::ArrayView<const LocalOrdinal> indices;
1036 Teuchos::ArrayView<const Scalar> values;
1037 A->getLocalRowView(i,indices,values);
1039 for (
size_t j=0; j<(size_t)indices.size(); j++) {
1040 if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > Teuchos::ScalarTraits<MT>::eps()) {
1044 if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
1045 dirichletRows.push_back(i);
1053 const std::vector<LocalOrdinal>& dirichletRows) {
1054 RCP<const Map> Rmap = A->getRowMap();
1055 RCP<const Map> Cmap = A->getColMap();
1056 Scalar one =Teuchos::ScalarTraits<Scalar>::one();
1057 Scalar zero =Teuchos::ScalarTraits<Scalar>::zero();
1059 for(
size_t i=0; i<dirichletRows.size(); i++) {
1060 GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
1062 Teuchos::ArrayView<const LocalOrdinal> indices;
1063 Teuchos::ArrayView<const Scalar> values;
1064 A->getLocalRowView(dirichletRows[i],indices,values);
1066 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
1067 for(
size_t j=0; j<(size_t)indices.size(); j++) {
1068 if(Cmap->getGlobalElement(indices[j])==row_gid)
1079 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
1080 TEUCHOS_ASSERT(A->isFillComplete());
1081 RCP<const Map> domMap = A->getDomainMap();
1082 RCP<const Map> ranMap = A->getRangeMap();
1083 RCP<const Map> Rmap = A->getRowMap();
1084 RCP<const Map> Cmap = A->getColMap();
1085 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getNodeNumElements());
1086 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1087 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1089 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1090 if (dirichletRows[i]){
1093 Teuchos::ArrayView<const LocalOrdinal> indices;
1094 Teuchos::ArrayView<const Scalar> values;
1095 A->getLocalRowView(i,indices,values);
1097 Teuchos::ArrayRCP<Scalar> valuesNC(values.size());
1098 for(
size_t j=0; j<(size_t)indices.size(); j++) {
1099 if(Cmap->getGlobalElement(indices[j])==row_gid)
1104 A->replaceLocalValues(i,indices,valuesNC());
1107 A->fillComplete(domMap, ranMap);
1112 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1113 const std::vector<LocalOrdinal>& dirichletRows,
1114 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1115 for(
size_t i=0; i<dirichletRows.size(); i++) {
1116 Teuchos::ArrayView<const LocalOrdinal> indices;
1117 Teuchos::ArrayView<const Scalar> values;
1118 A->getLocalRowView(dirichletRows[i],indices,values);
1120 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
1121 for(
size_t j=0; j<(size_t)indices.size(); j++)
1122 valuesNC[j]=replaceWith;
1128 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1129 const Teuchos::ArrayRCP<const bool>& dirichletRows,
1130 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1131 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == A->getRowMap()->getNodeNumElements());
1132 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1133 if (dirichletRows[i]) {
1134 Teuchos::ArrayView<const LocalOrdinal> indices;
1135 Teuchos::ArrayView<const Scalar> values;
1136 A->getLocalRowView(i,indices,values);
1138 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
1139 for(
size_t j=0; j<(size_t)indices.size(); j++)
1140 valuesNC[j]=replaceWith;
1147 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& X,
1148 const Teuchos::ArrayRCP<const bool>& dirichletRows,
1149 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1150 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == X->getMap()->getNodeNumElements());
1151 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1152 if (dirichletRows[i]) {
1153 for(
size_t j=0; j<X->getNumVectors(); j++)
1154 X->replaceLocalValue(i,j,replaceWith);
1162 const Teuchos::ArrayRCP<const bool>& dirichletCols,
1163 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1164 TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == A->getColMap()->getNodeNumElements());
1165 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
1166 Teuchos::ArrayView<const LocalOrdinal> indices;
1167 Teuchos::ArrayView<const Scalar> values;
1168 A->getLocalRowView(i,indices,values);
1170 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
1171 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
1172 if (dirichletCols[indices[j]])
1173 valuesNC[j] = replaceWith;
1179 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletRow,
1180 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletCol) {
1183 if(!A->getRowMap()->isSameAs(*A->getDomainMap())) {
1184 throw std::runtime_error(
"UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
1186 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A->getCrsGraph()->getImporter();
1187 bool has_import = !importer.is_null();
1190 std::vector<LocalOrdinal> dirichletRows;
1194 printf(
"[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
1195 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++)
1196 printf(
"%d ",dirichletRows[i]);
1201 isDirichletRow = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getRowMap(),
true);
1202 isDirichletCol = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getColMap(),
true);
1204 Teuchos::ArrayRCP<int> dr_rcp = isDirichletRow->getDataNonConst(0);
1205 Teuchos::ArrayView<int> dr = dr_rcp();
1206 Teuchos::ArrayRCP<int> dc_rcp = isDirichletCol->getDataNonConst(0);
1207 Teuchos::ArrayView<int> dc = dc_rcp();
1208 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1209 dr[dirichletRows[i]] = 1;
1210 if(!has_import) dc[dirichletRows[i]] = 1;
1214 isDirichletCol->doImport(*isDirichletRow,*importer,Xpetra::CombineMode::ADD);
1221 static RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> >
GeneratedBlockedTargetMap(
const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> & sourceBlockedMap,
1222 const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> & Importer) {
1223 typedef Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> IntVector;
1224 Xpetra::UnderlyingLib lib = sourceBlockedMap.lib();
1227 RCP<const Map> fullMap = sourceBlockedMap.getMap();
1228 RCP<const Map> stridedMap = Teuchos::rcp_dynamic_cast<
const Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >(fullMap);
1229 if(!stridedMap.is_null()) fullMap = stridedMap->getMap();
1232 const size_t numSubMaps = sourceBlockedMap.getNumMaps();
1233 if(!Importer.getSourceMap()->isCompatible(*fullMap))
1234 throw std::runtime_error(
"GenerateBlockedTargetMap(): Map compatibility error");
1237 RCP<IntVector> block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(fullMap);
1239 for(
size_t i=0; i<numSubMaps; i++) {
1240 RCP<const Map> map = sourceBlockedMap.getMap(i);
1242 for(
size_t j=0; j<map->getNodeNumElements(); j++) {
1243 LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
1244 block_ids->replaceLocalValue(jj,(
int)i);
1249 RCP<const Map> targetMap = Importer.getTargetMap();
1250 RCP<IntVector> new_block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(targetMap);
1251 new_block_ids->doImport(*block_ids,Importer,Xpetra::CombineMode::ADD);
1252 Teuchos::ArrayRCP<const int> dataRCP = new_block_ids->getData(0);
1253 Teuchos::ArrayView<const int> data = dataRCP();
1257 Teuchos::Array<Teuchos::Array<GlobalOrdinal> > elementsInSubMap(numSubMaps);
1258 for(
size_t i=0; i<targetMap->getNodeNumElements(); i++) {
1259 elementsInSubMap[data[i]].push_back(targetMap->getGlobalElement(i));
1263 std::vector<RCP<const Map> > subMaps(numSubMaps);
1264 for(
size_t i=0; i<numSubMaps; i++) {
1265 subMaps[i] = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(lib,Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),elementsInSubMap[i](),targetMap->getIndexBase(),targetMap->getComm());
1269 return rcp(
new BlockedMap(targetMap,subMaps));
1275 static bool MapsAreNested(
const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& rowMap,
const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& colMap) {
1276 ArrayView<const GlobalOrdinal> rowElements = rowMap.getNodeElementList();
1277 ArrayView<const GlobalOrdinal> colElements = colMap.getNodeElementList();
1279 const size_t numElements = rowElements.size();
1281 if (
size_t(colElements.size()) < numElements)
1284 bool goodMap =
true;
1285 for (
size_t i = 0; i < numElements; i++)
1286 if (rowElements[i] != colElements[i]) {
1304 #define MUELU_UTILITIESBASE_SHORT 1305 #endif // MUELU_UTILITIESBASE_DECL_HPP Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedMultiVector
MueLu::DefaultLocalOrdinal LocalOrdinal
static Teuchos::ArrayRCP< const bool > DetectDirichletCols(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
static void Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS, MultiVector &Resid)
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MultiVector
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Apply Rowsum Criterion.
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &BlockNumber)
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows)
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
static RCP< Vector > GetMatrixOverlappedDeletedRowsum(const Matrix &A)
Extract Overlapped Matrix Deleted Rowsum.
Namespace for MueLu classes and methods.
Exception throws to report incompatible objects (like maps).
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows (extended version)
#define MueLu_sumAll(rcpComm, in, out)
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i=k}(-a_ik), for each for i in the matrix.
static bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
static RCP< Xpetra::Vector< Magnitude, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedAbsDeletedRowsum(const Matrix &A)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > > &isDirichletRow, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > > &isDirichletCol)
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::ArrayView< double > &weight, const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Weighted squared distance between two rows in a multivector.
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false)
Extract Matrix Diagonal of lumped matrix.
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
Xpetra::BlockedVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedVector
Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrix
static Scalar Frobenius(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &B)
Frobenius inner product of two matrices.
Teuchos::ScalarTraits< Scalar >::magnitudeType Magnitude
Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Vector
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Extract Matrix Diagonal.
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &BlockNumber, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > BlockedMap
static void ZeroDirichletCols(Teuchos::RCP< Matrix > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
static RCP< Matrix > Crs2Op(RCP< CrsMatrix > Op)
Xpetra::CrsMatrixWrap< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrixWrap
static void FindNonZeros(const Teuchos::ArrayRCP< const Scalar > vals, Teuchos::ArrayRCP< bool > nonzeros)
Find non-zero values in an ArrayRCP Compares the value to 2 * machine epsilon.
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static Scalar PowerMethod(const Matrix &A, const RCP< Vector > &diagInvVec, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
Exception throws to report errors in the internal logical of the program.
Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > Map
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS, MultiVector &Resid)
static void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomain)
Detects Dirichlet columns & domains from a list of Dirichlet rows.
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &X, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.