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 Matrix> rcpA = Teuchos::rcpFromRef(A);
171 RCP<const Map> rowMap = rcpA->getRowMap();
172 RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
174 rcpA->getLocalDiagCopy(*diag);
191 size_t numRows = A.getRowMap()->getNodeNumElements();
192 Teuchos::ArrayRCP<Scalar> diag(numRows);
193 Teuchos::ArrayView<const LocalOrdinal> cols;
194 Teuchos::ArrayView<const Scalar> vals;
195 for (
size_t i = 0; i < numRows; ++i) {
196 A.getLocalRowView(i, cols, vals);
197 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
199 diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
213 RCP<Vector> diag = Teuchos::null;
215 RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
216 Teuchos::rcp_dynamic_cast<
const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
217 if(bA == Teuchos::null) {
218 RCP<const Map> rowMap = rcpA->getRowMap();
219 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
220 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
221 Teuchos::ArrayView<const LocalOrdinal> cols;
222 Teuchos::ArrayView<const Scalar> vals;
223 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
224 rcpA->getLocalRowView(i, cols, vals);
225 diagVals[i] = Teuchos::ScalarTraits<Scalar>::zero();
227 diagVals[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
235 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),
true);
237 for (
size_t row = 0; row < bA->Rows(); ++row) {
238 for (
size_t col = 0; col < bA->Cols(); ++col) {
239 if (!bA->getMatrix(row,col).is_null()) {
241 bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA->getMatrix(row,col)) == Teuchos::null);
242 RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag,row,bThyraMode);
244 ddtemp->update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
245 bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
263 static Teuchos::RCP<Vector>
GetInverse(Teuchos::RCP<const Vector> v,
Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100,
Scalar tolReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
265 RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),
true);
268 RCP<const BlockedVector> bv = Teuchos::rcp_dynamic_cast<
const BlockedVector>(v);
269 if(bv.is_null() ==
false) {
270 RCP<BlockedVector> bret = Teuchos::rcp_dynamic_cast<
BlockedVector>(ret);
271 TEUCHOS_TEST_FOR_EXCEPTION(bret.is_null() ==
true,
MueLu::Exceptions::RuntimeError,
"MueLu::UtilitiesBase::GetInverse: return vector should be of type BlockedVector");
272 RCP<const BlockedMap> bmap = bv->getBlockedMap();
273 for(
size_t r = 0; r < bmap->getNumMaps(); ++r) {
274 RCP<const MultiVector> submvec = bv->getMultiVector(r,bmap->getThyraMode());
275 RCP<const Vector> subvec = submvec->getVector(0);
277 bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
283 ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
284 ArrayRCP<const Scalar> inputVals = v->getData(0);
285 for (
size_t i = 0; i < v->getMap()->getNodeNumElements(); ++i) {
286 if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
287 retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
289 retVals[i] = tolReplacement;
302 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
305 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(rowMap);
306 if(!browMap.is_null()) rowMap = browMap->getMap();
308 RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
314 Teuchos::ArrayRCP<size_t> offsets;
315 crsOp->getLocalDiagOffsets(offsets);
316 crsOp->getLocalDiagCopy(*localDiag,offsets());
319 ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
322 localDiagVals[i] = diagVals[i];
323 localDiagVals = diagVals = null;
326 RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
327 RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
328 importer = A.getCrsGraph()->getImporter();
329 if (importer == Teuchos::null) {
330 importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
332 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
348 using STS =
typename Teuchos::ScalarTraits<SC>;
351 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
352 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(rowMap);
353 if(!browMap.is_null()) rowMap = browMap->getMap();
355 RCP<Vector> local = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(rowMap);
356 RCP<Vector> ghosted = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(colMap,
true);
357 ArrayRCP<SC> localVals = local->getDataNonConst(0);
359 for (LO row = 0; row < static_cast<LO>(A.getRowMap()->getNodeNumElements()); ++row) {
360 size_t nnz = A.getNumEntriesInLocalRow(row);
361 ArrayView<const LO> indices;
362 ArrayView<const SC> vals;
363 A.getLocalRowView(row, indices, vals);
367 for (LO colID = 0; colID < static_cast<LO>(nnz); colID++) {
368 if(indices[colID] != row) {
375 RCP< const Xpetra::Import<LO,GO,Node> > importer;
376 importer = A.getCrsGraph()->getImporter();
377 if (importer == Teuchos::null) {
378 importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
380 ghosted->doImport(*local, *(importer), Xpetra::INSERT);
386 static RCP<Xpetra::Vector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> >
388 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
389 using STS =
typename Teuchos::ScalarTraits<Scalar>;
390 using MTS =
typename Teuchos::ScalarTraits<Magnitude>;
395 using RealValuedVector = Xpetra::Vector<MT,LO,GO,Node>;
398 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<
const BlockedMap>(rowMap);
399 if(!browMap.is_null()) rowMap = browMap->getMap();
401 RCP<RealValuedVector> local = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(rowMap);
402 RCP<RealValuedVector> ghosted = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(colMap,
true);
403 ArrayRCP<MT> localVals = local->getDataNonConst(0);
405 for (LO rowIdx = 0; rowIdx < static_cast<LO>(A.getRowMap()->getNodeNumElements()); ++rowIdx) {
406 size_t nnz = A.getNumEntriesInLocalRow(rowIdx);
407 ArrayView<const LO> indices;
408 ArrayView<const SC> vals;
409 A.getLocalRowView(rowIdx, indices, vals);
413 for (LO colID = 0; colID < static_cast<LO>(nnz); ++colID) {
414 if(indices[colID] != rowIdx) {
415 si += STS::magnitude(vals[colID]);
418 localVals[rowIdx] = si;
421 RCP< const Xpetra::Import<LO,GO,Node> > importer;
422 importer = A.getCrsGraph()->getImporter();
423 if (importer == Teuchos::null) {
424 importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
426 ghosted->doImport(*local, *(importer), Xpetra::INSERT);
437 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
438 const size_t numVecs = X.getNumVectors();
439 RCP<MultiVector> RES =
Residual(Op, X, RHS);
440 Teuchos::Array<Magnitude> norms(numVecs);
446 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
447 const size_t numVecs = X.getNumVectors();
449 Teuchos::Array<Magnitude> norms(numVecs);
455 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
456 const size_t numVecs = X.getNumVectors();
458 RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(RHS.getMap(), numVecs,
false);
459 Op.residual(X,RHS,*RES);
465 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides");
466 TEUCHOS_TEST_FOR_EXCEPTION(Resid.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of residual vectors != number of right-hand sides");
467 Op.residual(X,RHS,Resid);
473 RCP<const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
474 int myPID = comm->getRank();
477 for (
int i = 0; i <comm->getSize(); i++) {
479 gethostname(hostname,
sizeof(hostname));
480 std::cout <<
"Host: " << hostname <<
"\tMPI rank: " << myPID <<
",\tPID: " << pid <<
"\n\tattach " << pid << std::endl;
485 std::cout <<
"** Enter a character to continue > " << std::endl;
487 int r = scanf(
"%c", &go);
515 LocalOrdinal niters = 10,
Magnitude tolerance = 1e-2,
bool verbose =
false,
unsigned int seed = 123) {
517 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
520 RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
521 RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
522 RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
527 Teuchos::Array<Magnitude> norms(1);
529 typedef Teuchos::ScalarTraits<Scalar> STS;
531 const Scalar zero = STS::zero(), one = STS::one();
534 Magnitude residual = STS::magnitude(zero);
537 RCP<Vector> diagInvVec;
539 RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
540 A.getLocalDiagCopy(*diagVec);
541 diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
542 diagInvVec->reciprocal(*diagVec);
545 for (
int iter = 0; iter < niters; ++iter) {
547 q->update(one/norms[0], *z, zero);
550 z->elementWiseMultiply(one, *diagInvVec, *z, zero);
553 if (iter % 100 == 0 || iter + 1 == niters) {
554 r->update(1.0, *z, -lambda, *q, zero);
556 residual = STS::magnitude(norms[0] / lambda);
558 std::cout <<
"Iter = " << iter
559 <<
" Lambda = " << lambda
560 <<
" Residual of A*q - lambda*q = " << residual
572 static RCP<Teuchos::FancyOStream>
MakeFancy(std::ostream& os) {
573 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
582 const size_t numVectors = v.size();
584 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
585 for (
size_t j = 0; j < numVectors; j++) {
586 d += (v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
588 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
603 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) {
605 typedef Teuchos::ScalarTraits<Scalar> STS;
606 ArrayRCP<bool> boundaryNodes(numRows,
true);
607 if (count_twos_as_dirichlet) {
609 ArrayView<const LocalOrdinal> indices;
610 ArrayView<const Scalar> vals;
611 A.getLocalRowView(row, indices, vals);
612 size_t nnz = A.getNumEntriesInLocalRow(row);
615 for (col = 0; col < nnz; col++)
616 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
617 if (!boundaryNodes[row])
619 boundaryNodes[row] =
false;
622 boundaryNodes[row] =
true;
627 ArrayView<const LocalOrdinal> indices;
628 ArrayView<const Scalar> vals;
629 A.getLocalRowView(row, indices, vals);
630 size_t nnz = A.getNumEntriesInLocalRow(row);
632 for (
size_t col = 0; col < nnz; col++)
633 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
634 boundaryNodes[row] =
false;
639 return boundaryNodes;
655 static Teuchos::ArrayRCP<const bool>
DetectDirichletRowsExt(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
bool & bHasZeroDiagonal,
const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
658 bHasZeroDiagonal =
false;
660 Teuchos::RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
661 A.getLocalDiagCopy(*diagVec);
662 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
665 typedef Teuchos::ScalarTraits<Scalar> STS;
666 ArrayRCP<bool> boundaryNodes(numRows,
false);
668 ArrayView<const LocalOrdinal> indices;
669 ArrayView<const Scalar> vals;
670 A.getLocalRowView(row, indices, vals);
672 bool bHasDiag =
false;
673 for (decltype(indices.size()) col = 0; col < indices.size(); col++) {
674 if ( indices[col] != row) {
675 if (STS::magnitude(vals[col] / STS::magnitude(sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col]))) ) > tol) {
678 }
else bHasDiag =
true;
680 if (bHasDiag ==
false) bHasZeroDiagonal =
true;
681 else if(nnz == 0) boundaryNodes[row] =
true;
683 return boundaryNodes;
697 static Teuchos::ArrayRCP<const bool>
DetectDirichletCols(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
698 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
699 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
700 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
701 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
702 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
703 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,1);
704 myColsToZero->putScalar(zero);
706 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
707 if (dirichletRows[i]) {
708 Teuchos::ArrayView<const LocalOrdinal> indices;
709 Teuchos::ArrayView<const Scalar> values;
710 A.getLocalRowView(i,indices,values);
711 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
712 myColsToZero->replaceLocalValue(indices[j],0,one);
716 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,1);
717 globalColsToZero->putScalar(zero);
718 Teuchos::RCP<Xpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
720 globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
722 myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
723 Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
724 Teuchos::ArrayRCP<bool> dirichletCols(colMap->getNodeNumElements(),
true);
725 Magnitude eps = Teuchos::ScalarTraits<Magnitude>::eps();
726 for(
size_t i=0; i<colMap->getNodeNumElements(); i++) {
727 dirichletCols[i] = Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i])>2.0*eps;
729 return dirichletCols;
737 static Scalar Frobenius(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
742 TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()),
Exceptions::Incompatible,
"MueLu::CGSolver::Frobenius: row maps are incompatible");
743 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(),
Exceptions::RuntimeError,
"Matrices must be fill completed");
745 const Map& AColMap = *A.getColMap();
746 const Map& BColMap = *B.getColMap();
748 Teuchos::ArrayView<const LocalOrdinal> indA, indB;
749 Teuchos::ArrayView<const Scalar> valA, valB;
750 size_t nnzA = 0, nnzB = 0;
762 Teuchos::Array<Scalar> valBAll(BColMap.getNodeNumElements());
764 LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
765 Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
766 size_t numRows = A.getNodeNumRows();
767 for (
size_t i = 0; i < numRows; i++) {
768 A.getLocalRowView(i, indA, valA);
769 B.getLocalRowView(i, indB, valB);
774 for (
size_t j = 0; j < nnzB; j++)
775 valBAll[indB[j]] = valB[j];
777 for (
size_t j = 0; j < nnzA; j++) {
780 LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
782 f += valBAll[ind] * valA[j];
786 for (
size_t j = 0; j < nnzB; j++)
787 valBAll[indB[j]] = zero;
810 int maxint = INT_MAX;
811 int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.getRank()+1)/(comm.getSize()+one)) );
812 if (mySeed < 1 || mySeed == maxint) {
813 std::ostringstream errStr;
814 errStr <<
"Error detected with random seed = " << mySeed <<
". It should be in the interval [1,2^31-2].";
819 Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
831 static void FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
832 std::vector<LocalOrdinal>& dirichletRows,
bool count_twos_as_dirichlet=
false) {
833 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
834 dirichletRows.resize(0);
835 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
836 Teuchos::ArrayView<const LocalOrdinal> indices;
837 Teuchos::ArrayView<const Scalar> values;
838 A->getLocalRowView(i,indices,values);
840 for (
size_t j=0; j<(size_t)indices.size(); j++) {
841 if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > Teuchos::ScalarTraits<MT>::eps()) {
845 if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
846 dirichletRows.push_back(i);
854 const std::vector<LocalOrdinal>& dirichletRows) {
855 RCP<const Map> Rmap = A->getRowMap();
856 RCP<const Map> Cmap = A->getColMap();
857 Scalar one =Teuchos::ScalarTraits<Scalar>::one();
858 Scalar zero =Teuchos::ScalarTraits<Scalar>::zero();
860 for(
size_t i=0; i<dirichletRows.size(); i++) {
861 GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
863 Teuchos::ArrayView<const LocalOrdinal> indices;
864 Teuchos::ArrayView<const Scalar> values;
865 A->getLocalRowView(dirichletRows[i],indices,values);
867 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
868 for(
size_t j=0; j<(size_t)indices.size(); j++) {
869 if(Cmap->getGlobalElement(indices[j])==row_gid)
880 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
881 TEUCHOS_ASSERT(A->isFillComplete());
882 RCP<const Map> domMap = A->getDomainMap();
883 RCP<const Map> ranMap = A->getRangeMap();
884 RCP<const Map> Rmap = A->getRowMap();
885 RCP<const Map> Cmap = A->getColMap();
886 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getNodeNumElements());
887 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
888 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
890 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
891 if (dirichletRows[i]){
894 Teuchos::ArrayView<const LocalOrdinal> indices;
895 Teuchos::ArrayView<const Scalar> values;
896 A->getLocalRowView(i,indices,values);
898 Teuchos::ArrayRCP<Scalar> valuesNC(values.size());
899 for(
size_t j=0; j<(size_t)indices.size(); j++) {
900 if(Cmap->getGlobalElement(indices[j])==row_gid)
905 A->replaceLocalValues(i,indices,valuesNC());
908 A->fillComplete(domMap, ranMap);
913 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
914 const std::vector<LocalOrdinal>& dirichletRows,
915 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
916 for(
size_t i=0; i<dirichletRows.size(); i++) {
917 Teuchos::ArrayView<const LocalOrdinal> indices;
918 Teuchos::ArrayView<const Scalar> values;
919 A->getLocalRowView(dirichletRows[i],indices,values);
921 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
922 for(
size_t j=0; j<(size_t)indices.size(); j++)
923 valuesNC[j]=replaceWith;
929 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
930 const Teuchos::ArrayRCP<const bool>& dirichletRows,
931 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
932 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == A->getRowMap()->getNodeNumElements());
933 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
934 if (dirichletRows[i]) {
935 Teuchos::ArrayView<const LocalOrdinal> indices;
936 Teuchos::ArrayView<const Scalar> values;
937 A->getLocalRowView(i,indices,values);
939 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
940 for(
size_t j=0; j<(size_t)indices.size(); j++)
941 valuesNC[j]=replaceWith;
948 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& X,
949 const Teuchos::ArrayRCP<const bool>& dirichletRows,
950 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
951 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == X->getMap()->getNodeNumElements());
952 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
953 if (dirichletRows[i]) {
954 for(
size_t j=0; j<X->getNumVectors(); j++)
955 X->replaceLocalValue(i,j,replaceWith);
963 const Teuchos::ArrayRCP<const bool>& dirichletCols,
964 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
965 TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == A->getColMap()->getNodeNumElements());
966 for(
size_t i=0; i<A->getNodeNumRows(); i++) {
967 Teuchos::ArrayView<const LocalOrdinal> indices;
968 Teuchos::ArrayView<const Scalar> values;
969 A->getLocalRowView(i,indices,values);
971 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
972 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
973 if (dirichletCols[indices[j]])
974 valuesNC[j] = replaceWith;
980 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletRow,
981 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletCol) {
984 if(!A->getRowMap()->isSameAs(*A->getDomainMap())) {
985 throw std::runtime_error(
"UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
987 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A->getCrsGraph()->getImporter();
988 bool has_import = !importer.is_null();
991 std::vector<LocalOrdinal> dirichletRows;
995 printf(
"[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
996 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++)
997 printf(
"%d ",dirichletRows[i]);
1002 isDirichletRow = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getRowMap(),
true);
1003 isDirichletCol = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getColMap(),
true);
1005 Teuchos::ArrayRCP<int> dr_rcp = isDirichletRow->getDataNonConst(0);
1006 Teuchos::ArrayView<int> dr = dr_rcp();
1007 Teuchos::ArrayRCP<int> dc_rcp = isDirichletCol->getDataNonConst(0);
1008 Teuchos::ArrayView<int> dc = dc_rcp();
1009 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1010 dr[dirichletRows[i]] = 1;
1011 if(!has_import) dc[dirichletRows[i]] = 1;
1015 isDirichletCol->doImport(*isDirichletRow,*importer,Xpetra::CombineMode::ADD);
1022 static RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> >
GeneratedBlockedTargetMap(
const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> & sourceBlockedMap,
1023 const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> & Importer) {
1024 typedef Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> IntVector;
1025 Xpetra::UnderlyingLib lib = sourceBlockedMap.lib();
1028 RCP<const Map> fullMap = sourceBlockedMap.getMap();
1029 RCP<const Map> stridedMap = Teuchos::rcp_dynamic_cast<
const Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >(fullMap);
1030 if(!stridedMap.is_null()) fullMap = stridedMap->getMap();
1033 const size_t numSubMaps = sourceBlockedMap.getNumMaps();
1034 if(!Importer.getSourceMap()->isCompatible(*fullMap))
1035 throw std::runtime_error(
"GenerateBlockedTargetMap(): Map compatibility error");
1038 RCP<IntVector> block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(fullMap);
1040 for(
size_t i=0; i<numSubMaps; i++) {
1041 RCP<const Map> map = sourceBlockedMap.getMap(i);
1043 for(
size_t j=0; j<map->getNodeNumElements(); j++) {
1044 LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
1045 block_ids->replaceLocalValue(jj,(
int)i);
1050 RCP<const Map> targetMap = Importer.getTargetMap();
1051 RCP<IntVector> new_block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(targetMap);
1052 new_block_ids->doImport(*block_ids,Importer,Xpetra::CombineMode::ADD);
1053 Teuchos::ArrayRCP<const int> dataRCP = new_block_ids->getData(0);
1054 Teuchos::ArrayView<const int> data = dataRCP();
1058 Teuchos::Array<Teuchos::Array<GlobalOrdinal> > elementsInSubMap(numSubMaps);
1059 for(
size_t i=0; i<targetMap->getNodeNumElements(); i++) {
1060 elementsInSubMap[data[i]].push_back(targetMap->getGlobalElement(i));
1064 std::vector<RCP<const Map> > subMaps(numSubMaps);
1065 for(
size_t i=0; i<numSubMaps; i++) {
1066 subMaps[i] = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(lib,Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),elementsInSubMap[i](),targetMap->getIndexBase(),targetMap->getComm());
1070 return rcp(
new BlockedMap(targetMap,subMaps));
1081 #define MUELU_UTILITIESBASE_SHORT 1082 #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)
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)
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)
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 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)
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.