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> 92 #define MueLu_sumAll(rcpComm, in, out) \ 93 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out)) 94 #define MueLu_minAll(rcpComm, in, out) \ 95 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out)) 96 #define MueLu_maxAll(rcpComm, in, out) \ 97 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out)) 106 template <
class Scalar,
107 class LocalOrdinal = int,
108 class GlobalOrdinal = LocalOrdinal,
109 class Node = KokkosClassic::DefaultNode::DefaultNodeType>
110 class UtilitiesBase {
112 #undef MUELU_UTILITIESBASE_SHORT 115 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>
CrsMatrixWrap;
116 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
CrsMatrix;
117 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
Matrix;
118 typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
Vector;
119 typedef Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
BlockedVector;
120 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
MultiVector;
122 typedef Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>
BlockedMap;
123 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>
Map;
125 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Magnitude;
128 static RCP<Matrix>
Crs2Op(RCP<CrsMatrix> Op) {
130 return Teuchos::null;
141 size_t numRows = A.getRowMap()->getNodeNumElements();
142 Teuchos::ArrayRCP<Scalar> diag(numRows);
143 Teuchos::ArrayView<const LocalOrdinal> cols;
144 Teuchos::ArrayView<const Scalar> vals;
145 for (
size_t i = 0; i < numRows; ++i) {
146 A.getLocalRowView(i, cols, vals);
148 for (; j < cols.size(); ++j) {
149 if (Teuchos::as<size_t>(cols[j]) == i) {
154 if (j == cols.size()) {
156 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
170 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
172 RCP<const Map> rowMap = rcpA->getRowMap();
173 RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
175 rcpA->getLocalDiagCopy(*diag);
192 size_t numRows = A.getRowMap()->getNodeNumElements();
193 Teuchos::ArrayRCP<Scalar> diag(numRows);
194 Teuchos::ArrayView<const LocalOrdinal> cols;
195 Teuchos::ArrayView<const Scalar> vals;
196 for (
size_t i = 0; i < numRows; ++i) {
197 A.getLocalRowView(i, cols, vals);
198 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
199 for (LocalOrdinal j = 0; j < cols.size(); ++j) {
200 diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
214 RCP<Vector> diag = Teuchos::null;
216 RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
217 Teuchos::rcp_dynamic_cast<
const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
218 if(bA == Teuchos::null) {
219 RCP<const Map> rowMap = rcpA->getRowMap();
220 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
221 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
222 Teuchos::ArrayView<const LocalOrdinal> cols;
223 Teuchos::ArrayView<const Scalar> vals;
224 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
225 rcpA->getLocalRowView(i, cols, vals);
226 diagVals[i] = Teuchos::ScalarTraits<Scalar>::zero();
227 for (LocalOrdinal j = 0; j < cols.size(); ++j) {
228 diagVals[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
236 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),
true);
238 for (
size_t row = 0; row < bA->Rows(); ++row) {
239 for (
size_t col = 0; col < bA->Cols(); ++col) {
240 if (!bA->getMatrix(row,col).is_null()) {
242 bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA->getMatrix(row,col)) == Teuchos::null);
243 RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag,row,bThyraMode);
245 ddtemp->update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
246 bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
264 static Teuchos::RCP<Vector>
GetInverse(Teuchos::RCP<const Vector> v,
Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar tolReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
266 RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),
true);
269 RCP<const BlockedVector> bv = Teuchos::rcp_dynamic_cast<
const BlockedVector>(v);
270 if(bv.is_null() ==
false) {
271 RCP<BlockedVector> bret = Teuchos::rcp_dynamic_cast<
BlockedVector>(ret);
272 TEUCHOS_TEST_FOR_EXCEPTION(bret.is_null() ==
true,
MueLu::Exceptions::RuntimeError,
"MueLu::UtilitiesBase::GetInverse: return vector should be of type BlockedVector");
273 RCP<const BlockedMap> bmap = bv->getBlockedMap();
274 for(
size_t r = 0; r < bmap->getNumMaps(); ++r) {
275 RCP<const MultiVector> submvec = bv->getMultiVector(r,bmap->getThyraMode());
276 RCP<const Vector> subvec = submvec->getVector(0);
278 bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
284 ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
285 ArrayRCP<const Scalar> inputVals = v->getData(0);
286 for (
size_t i = 0; i < v->getMap()->getNodeNumElements(); ++i) {
287 if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
288 retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
290 retVals[i] = tolReplacement;
303 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
306 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(rowMap);
307 if(!browMap.is_null()) rowMap = browMap->getMap();
309 RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
315 Teuchos::ArrayRCP<size_t> offsets;
316 crsOp->getLocalDiagOffsets(offsets);
317 crsOp->getLocalDiagCopy(*localDiag,offsets());
320 ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
322 for (LocalOrdinal i = 0; i < localDiagVals.size(); i++)
323 localDiagVals[i] = diagVals[i];
324 localDiagVals = diagVals = null;
327 RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
328 RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
329 importer = A.getCrsGraph()->getImporter();
330 if (importer == Teuchos::null) {
331 importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
333 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
343 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
344 const size_t numVecs = X.getNumVectors();
345 RCP<MultiVector> RES =
Residual(Op, X, RHS);
346 Teuchos::Array<Magnitude> norms(numVecs);
352 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
353 const size_t numVecs = X.getNumVectors();
355 Teuchos::Array<Magnitude> norms(numVecs);
361 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
362 const size_t numVecs = X.getNumVectors();
363 Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
365 RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(RHS.getMap(), numVecs,
false);
366 Op.apply(X, *RES, Teuchos::NO_TRANS, one, zero);
367 RES->update(one, RHS, negone);
373 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides");
374 TEUCHOS_TEST_FOR_EXCEPTION(Resid.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of residual vectors != number of right-hand sides");
375 Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
377 Op.apply(X, Resid, Teuchos::NO_TRANS, one, zero);
378 Resid.update(one, RHS, negone);
384 RCP<const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
385 int myPID = comm->getRank();
388 for (
int i = 0; i <comm->getSize(); i++) {
390 gethostname(hostname,
sizeof(hostname));
391 std::cout <<
"Host: " << hostname <<
"\tMPI rank: " << myPID <<
",\tPID: " << pid <<
"\n\tattach " << pid << std::endl;
396 std::cout <<
"** Enter a character to continue > " << std::endl;
398 int r = scanf(
"%c", &go);
426 LocalOrdinal niters = 10,
Magnitude tolerance = 1e-2,
bool verbose =
false,
unsigned int seed = 123) {
428 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
431 RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
432 RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
433 RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
438 Teuchos::Array<Magnitude> norms(1);
440 typedef Teuchos::ScalarTraits<Scalar> STS;
442 const Scalar zero = STS::zero(), one = STS::one();
444 Scalar lambda = zero;
445 Magnitude residual = STS::magnitude(zero);
448 RCP<Vector> diagInvVec;
450 RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
451 A.getLocalDiagCopy(*diagVec);
452 diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
453 diagInvVec->reciprocal(*diagVec);
456 for (
int iter = 0; iter < niters; ++iter) {
458 q->update(one/norms[0], *z, zero);
461 z->elementWiseMultiply(one, *diagInvVec, *z, zero);
464 if (iter % 100 == 0 || iter + 1 == niters) {
465 r->update(1.0, *z, -lambda, *q, zero);
467 residual = STS::magnitude(norms[0] / lambda);
469 std::cout <<
"Iter = " << iter
470 <<
" Lambda = " << lambda
471 <<
" Residual of A*q - lambda*q = " << residual
475 if (residual < tolerance)
483 static RCP<Teuchos::FancyOStream>
MakeFancy(std::ostream& os) {
484 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
492 static typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Distance2(
const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>> &v, LocalOrdinal i0, LocalOrdinal i1) {
493 const size_t numVectors = v.size();
495 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
496 for (
size_t j = 0; j < numVectors; j++) {
497 d += (v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
499 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
514 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) {
515 LocalOrdinal numRows = A.getNodeNumRows();
516 typedef Teuchos::ScalarTraits<Scalar> STS;
517 ArrayRCP<bool> boundaryNodes(numRows,
true);
518 if (count_twos_as_dirichlet) {
519 for (LocalOrdinal row = 0; row < numRows; row++) {
520 ArrayView<const LocalOrdinal> indices;
521 ArrayView<const Scalar> vals;
522 A.getLocalRowView(row, indices, vals);
523 size_t nnz = A.getNumEntriesInLocalRow(row);
526 for (col = 0; col < nnz; col++)
527 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
528 if (!boundaryNodes[row])
530 boundaryNodes[row] =
false;
533 boundaryNodes[row] =
true;
537 for (LocalOrdinal row = 0; row < numRows; row++) {
538 ArrayView<const LocalOrdinal> indices;
539 ArrayView<const Scalar> vals;
540 A.getLocalRowView(row, indices, vals);
541 size_t nnz = A.getNumEntriesInLocalRow(row);
543 for (
size_t col = 0; col < nnz; col++)
544 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
545 boundaryNodes[row] =
false;
550 return boundaryNodes;
566 static Teuchos::ArrayRCP<const bool>
DetectDirichletRowsExt(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
bool & bHasZeroDiagonal,
const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
569 bHasZeroDiagonal =
false;
571 Teuchos::RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
572 A.getLocalDiagCopy(*diagVec);
573 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
575 LocalOrdinal numRows = A.getNodeNumRows();
576 typedef Teuchos::ScalarTraits<Scalar> STS;
577 ArrayRCP<bool> boundaryNodes(numRows,
false);
578 for (LocalOrdinal row = 0; row < numRows; row++) {
579 ArrayView<const LocalOrdinal> indices;
580 ArrayView<const Scalar> vals;
581 A.getLocalRowView(row, indices, vals);
583 bool bHasDiag =
false;
584 for (decltype(indices.size()) col = 0; col < indices.size(); col++) {
585 if ( indices[col] != row) {
586 if (STS::magnitude(vals[col] / sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col])) ) > tol) {
589 }
else bHasDiag =
true;
591 if (bHasDiag ==
false) bHasZeroDiagonal =
true;
592 else if(nnz == 0) boundaryNodes[row] =
true;
594 return boundaryNodes;
608 static Teuchos::ArrayRCP<const bool>
DetectDirichletCols(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
609 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
610 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
611 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
612 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
613 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
614 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,1);
615 myColsToZero->putScalar(zero);
617 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
618 if (dirichletRows[i]) {
619 Teuchos::ArrayView<const LocalOrdinal> indices;
620 Teuchos::ArrayView<const Scalar> values;
621 A.getLocalRowView(i,indices,values);
622 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
623 myColsToZero->replaceLocalValue(indices[j],0,one);
627 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,1);
628 globalColsToZero->putScalar(zero);
629 Teuchos::RCP<Xpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
631 globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
633 myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
634 Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
635 Teuchos::ArrayRCP<bool> dirichletCols(colMap->getNodeNumElements(),
true);
636 Magnitude eps = Teuchos::ScalarTraits<Magnitude>::eps();
637 for(
size_t i=0; i<colMap->getNodeNumElements(); i++) {
638 dirichletCols[i] = Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i])>2.0*eps;
640 return dirichletCols;
648 static Scalar
Frobenius(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
653 TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()),
Exceptions::Incompatible,
"MueLu::CGSolver::Frobenius: row maps are incompatible");
654 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(),
Exceptions::RuntimeError,
"Matrices must be fill completed");
656 const Map& AColMap = *A.getColMap();
657 const Map& BColMap = *B.getColMap();
659 Teuchos::ArrayView<const LocalOrdinal> indA, indB;
660 Teuchos::ArrayView<const Scalar> valA, valB;
661 size_t nnzA = 0, nnzB = 0;
673 Teuchos::Array<Scalar> valBAll(BColMap.getNodeNumElements());
675 LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
676 Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
677 size_t numRows = A.getNodeNumRows();
678 for (
size_t i = 0; i < numRows; i++) {
679 A.getLocalRowView(i, indA, valA);
680 B.getLocalRowView(i, indB, valB);
685 for (
size_t j = 0; j < nnzB; j++)
686 valBAll[indB[j]] = valB[j];
688 for (
size_t j = 0; j < nnzA; j++) {
691 LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
693 f += valBAll[ind] * valA[j];
697 for (
size_t j = 0; j < nnzB; j++)
698 valBAll[indB[j]] = zero;
721 int maxint = INT_MAX;
722 int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.getRank()+1)/(comm.getSize()+one)) );
723 if (mySeed < 1 || mySeed == maxint) {
724 std::ostringstream errStr;
725 errStr <<
"Error detected with random seed = " << mySeed <<
". It should be in the interval [1,2^31-2].";
730 Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
742 static void FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
743 std::vector<LocalOrdinal>& dirichletRows,
bool count_twos_as_dirichlet=
false) {
744 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
745 dirichletRows.resize(0);
746 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
747 Teuchos::ArrayView<const LocalOrdinal> indices;
748 Teuchos::ArrayView<const Scalar> values;
749 A->getLocalRowView(i,indices,values);
751 for (
size_t j=0; j<(size_t)indices.size(); j++) {
752 if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > Teuchos::ScalarTraits<MT>::eps()) {
756 if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
757 dirichletRows.push_back(i);
765 const std::vector<LocalOrdinal>& dirichletRows) {
766 RCP<const Map> Rmap = A->getColMap();
767 RCP<const Map> Cmap = A->getColMap();
768 Scalar one =Teuchos::ScalarTraits<Scalar>::one();
769 Scalar zero =Teuchos::ScalarTraits<Scalar>::zero();
771 for(
size_t i=0; i<dirichletRows.size(); i++) {
772 GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
774 Teuchos::ArrayView<const LocalOrdinal> indices;
775 Teuchos::ArrayView<const Scalar> values;
776 A->getLocalRowView(dirichletRows[i],indices,values);
778 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
779 for(
size_t j=0; j<(size_t)indices.size(); j++) {
780 if(Cmap->getGlobalElement(indices[j])==row_gid)
791 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
792 RCP<const Map> Rmap = A->getColMap();
793 RCP<const Map> Cmap = A->getColMap();
794 Scalar one =Teuchos::ScalarTraits<Scalar>::one();
795 Scalar zero =Teuchos::ScalarTraits<Scalar>::zero();
797 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
798 if (dirichletRows[i]){
799 GlobalOrdinal row_gid = Rmap->getGlobalElement(i);
801 Teuchos::ArrayView<const LocalOrdinal> indices;
802 Teuchos::ArrayView<const Scalar> values;
803 A->getLocalRowView(i,indices,values);
805 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
806 for(
size_t j=0; j<(size_t)indices.size(); j++) {
807 if(Cmap->getGlobalElement(indices[j])==row_gid)
818 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
819 const std::vector<LocalOrdinal>& dirichletRows,
820 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
821 for(
size_t i=0; i<dirichletRows.size(); i++) {
822 Teuchos::ArrayView<const LocalOrdinal> indices;
823 Teuchos::ArrayView<const Scalar> values;
824 A->getLocalRowView(dirichletRows[i],indices,values);
826 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
827 for(
size_t j=0; j<(size_t)indices.size(); j++)
828 valuesNC[j]=replaceWith;
834 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
835 const Teuchos::ArrayRCP<const bool>& dirichletRows,
836 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
837 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
838 if (dirichletRows[i]) {
839 Teuchos::ArrayView<const LocalOrdinal> indices;
840 Teuchos::ArrayView<const Scalar> values;
841 A->getLocalRowView(i,indices,values);
843 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
844 for(
size_t j=0; j<(size_t)indices.size(); j++)
845 valuesNC[j]=replaceWith;
852 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& X,
853 const Teuchos::ArrayRCP<const bool>& dirichletRows,
854 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
855 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
856 if (dirichletRows[i]) {
857 for(
size_t j=0; j<X->getNumVectors(); j++)
858 X->replaceLocalValue(i,j,replaceWith);
866 const Teuchos::ArrayRCP<const bool>& dirichletCols,
867 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
868 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
869 Teuchos::ArrayView<const LocalOrdinal> indices;
870 Teuchos::ArrayView<const Scalar> values;
871 A->getLocalRowView(i,indices,values);
873 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
874 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
875 if (dirichletCols[indices[j]])
876 valuesNC[j] = replaceWith;
882 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletRow,
883 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletCol) {
886 if(!A->getRowMap()->isSameAs(*A->getDomainMap())) {
887 throw std::runtime_error(
"UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
889 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A->getCrsGraph()->getImporter();
890 bool has_import = !importer.is_null();
893 std::vector<LocalOrdinal> dirichletRows;
897 printf(
"[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
898 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++)
899 printf(
"%d ",dirichletRows[i]);
904 isDirichletRow = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getRowMap(),
true);
905 isDirichletCol = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getColMap(),
true);
907 Teuchos::ArrayRCP<int> dr_rcp = isDirichletRow->getDataNonConst(0);
908 Teuchos::ArrayView<int> dr = dr_rcp();
909 Teuchos::ArrayRCP<int> dc_rcp = isDirichletCol->getDataNonConst(0);
910 Teuchos::ArrayView<int> dc = dc_rcp();
911 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
912 dr[dirichletRows[i]] = 1;
913 if(!has_import) dc[dirichletRows[i]] = 1;
917 isDirichletCol->doImport(*isDirichletRow,*importer,Xpetra::CombineMode::ADD);
924 static RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> >
GeneratedBlockedTargetMap(
const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> & sourceBlockedMap,
925 const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> & Importer) {
926 typedef Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> IntVector;
927 Xpetra::UnderlyingLib lib = sourceBlockedMap.lib();
930 RCP<const Map> fullMap = sourceBlockedMap.getMap();
931 RCP<const Map> stridedMap = Teuchos::rcp_dynamic_cast<
const Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >(fullMap);
932 if(!stridedMap.is_null()) fullMap = stridedMap->getMap();
935 const size_t numSubMaps = sourceBlockedMap.getNumMaps();
936 if(!Importer.getSourceMap()->isCompatible(*fullMap))
937 throw std::runtime_error(
"GenerateBlockedTargetMap(): Map compatibility error");
940 RCP<IntVector> block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(fullMap);
942 for(
size_t i=0; i<numSubMaps; i++) {
943 RCP<const Map> map = sourceBlockedMap.getMap(i);
945 for(
size_t j=0; j<map->getNodeNumElements(); j++) {
946 LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
947 block_ids->replaceLocalValue(jj,(
int)i);
952 RCP<const Map> targetMap = Importer.getTargetMap();
953 RCP<IntVector> new_block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(targetMap);
954 new_block_ids->doImport(*block_ids,Importer,Xpetra::CombineMode::ADD);
955 Teuchos::ArrayRCP<const int> dataRCP = new_block_ids->getData(0);
956 Teuchos::ArrayView<const int> data = dataRCP();
960 Teuchos::Array<Teuchos::Array<GlobalOrdinal> > elementsInSubMap(numSubMaps);
961 for(
size_t i=0; i<targetMap->getNodeNumElements(); i++) {
962 elementsInSubMap[data[i]].push_back(targetMap->getGlobalElement(i));
966 std::vector<RCP<const Map> > subMaps(numSubMaps);
967 for(
size_t i=0; i<numSubMaps; i++) {
968 subMaps[i] = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(lib,Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),elementsInSubMap[i](),targetMap->getIndexBase(),targetMap->getComm(),targetMap->getNode());
972 return rcp(
new BlockedMap(targetMap,subMaps));
984 #define MUELU_UTILITIESBASE_SHORT 985 #endif // MUELU_UTILITIESBASE_DECL_HPP Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedMultiVector
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.
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Simple transpose for Tpetra::CrsMatrix types.
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 RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows)
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
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)
static void PauseForDebugger()
#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 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)
Xpetra::BlockedVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedVector
Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrix
static Teuchos::ArrayRCP< Scalar > GetLumpedMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal of lumped matrix.
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 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 RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100)
Extract Matrix Diagonal.
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
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 Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Teuchos::RCP< const Matrix > rcpA)
Extract Matrix Diagonal of lumped matrix.
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.