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_CrsGraphFactory_fwd.hpp>
63#include <Xpetra_CrsGraph_fwd.hpp>
64#include <Xpetra_CrsMatrix_fwd.hpp>
65#include <Xpetra_CrsMatrixWrap_fwd.hpp>
66#include <Xpetra_Map_fwd.hpp>
67#include <Xpetra_BlockedMap_fwd.hpp>
68#include <Xpetra_MapFactory_fwd.hpp>
69#include <Xpetra_Matrix_fwd.hpp>
70#include <Xpetra_MatrixFactory_fwd.hpp>
71#include <Xpetra_MultiVector_fwd.hpp>
72#include <Xpetra_MultiVectorFactory_fwd.hpp>
73#include <Xpetra_Operator_fwd.hpp>
74#include <Xpetra_Vector_fwd.hpp>
75#include <Xpetra_BlockedMultiVector.hpp>
76#include <Xpetra_BlockedVector.hpp>
77#include <Xpetra_VectorFactory_fwd.hpp>
78#include <Xpetra_ExportFactory.hpp>
80#include <Xpetra_Import.hpp>
81#include <Xpetra_ImportFactory.hpp>
82#include <Xpetra_MatrixMatrix.hpp>
83#include <Xpetra_CrsGraph.hpp>
84#include <Xpetra_CrsGraphFactory.hpp>
85#include <Xpetra_CrsMatrixWrap.hpp>
86#include <Xpetra_StridedMap.hpp>
94#define MueLu_sumAll(rcpComm, in, out) \
95 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out))
96#define MueLu_minAll(rcpComm, in, out) \
97 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out))
98#define MueLu_maxAll(rcpComm, in, out) \
99 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out))
114#undef MUELU_UTILITIESBASE_SHORT
117 using CrsGraph = Xpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node>;
119 using CrsMatrixWrap = Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
120 using CrsMatrix = Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
121 using Matrix = Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
122 using Vector = Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
123 using BlockedVector = Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
124 using MultiVector = Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
126 using BlockedMap = Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node>;
127 using Map = Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>;
129 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Magnitude;
132 static RCP<Matrix>
Crs2Op(RCP<CrsMatrix> Op) {
134 return Teuchos::null;
146 RCP<const Map> rowmap = Ain->getRowMap();
147 RCP<const Map> colmap = Ain->getColMap();
148 RCP<CrsMatrixWrap> Aout = rcp(
new CrsMatrixWrap(rowmap, expectedNNZperRow <= 0 ? Ain->getGlobalMaxNumRowEntries() : expectedNNZperRow));
150 for(
size_t row=0; row<Ain->getLocalNumRows(); row++)
152 size_t nnz = Ain->getNumEntriesInLocalRow(row);
154 Teuchos::ArrayView<const LocalOrdinal> indices;
155 Teuchos::ArrayView<const Scalar> vals;
156 Ain->getLocalRowView(row, indices, vals);
158 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz,
Exceptions::RuntimeError,
"MueLu::ThresholdAFilterFactory::Build: number of nonzeros not equal to number of indices? Error.");
160 Teuchos::ArrayRCP<GlobalOrdinal> indout(indices.size(),Teuchos::ScalarTraits<GlobalOrdinal>::zero());
161 Teuchos::ArrayRCP<Scalar> valout(indices.size(),Teuchos::ScalarTraits<Scalar>::zero());
162 size_t nNonzeros = 0;
165 LocalOrdinal lclColIdx = colmap->getLocalElement(glbRow);
166 for(
size_t i=0; i<(size_t)indices.size(); i++) {
167 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(threshold) || indices[i]==lclColIdx) {
168 indout[nNonzeros] = colmap->getGlobalElement(indices[i]);
169 valout[nNonzeros] = vals[i];
174 for(
size_t i=0; i<(size_t)indices.size(); i++) {
175 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(threshold)) {
176 indout[nNonzeros] = colmap->getGlobalElement(indices[i]);
177 valout[nNonzeros] = vals[i];
182 indout.resize(nNonzeros);
183 valout.resize(nNonzeros);
185 Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
187 Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());
200 using STS = Teuchos::ScalarTraits<Scalar>;
201 RCP<CrsGraph> sparsityPattern = CrsGraphFactory::Build(A->getRowMap(), expectedNNZperRow <= 0 ? A->getGlobalMaxNumRowEntries() : expectedNNZperRow);
204 ArrayRCP<const Scalar> D = diag->getData(0);
206 for(
size_t row=0; row<A->getLocalNumRows(); row++)
208 ArrayView<const LocalOrdinal> indices;
209 ArrayView<const Scalar> vals;
210 A->getLocalRowView(row, indices, vals);
212 GlobalOrdinal globalRow = A->getRowMap()->getGlobalElement(row);
213 LocalOrdinal col = A->getColMap()->getLocalElement(globalRow);
215 const Scalar Dk = STS::magnitude(D[col]) > 0.0 ? STS::magnitude(D[col]) : 1.0;
216 Array<GlobalOrdinal> indicesNew;
218 for(
size_t i=0; i<size_t(indices.size()); i++)
220 if(col == indices[i] || STS::magnitude(STS::squareroot(Dk)*vals[i]*STS::squareroot(Dk)) > STS::magnitude(threshold))
221 indicesNew.append(A->getColMap()->getGlobalElement(indices[i]));
223 sparsityPattern->insertGlobalIndices(globalRow, ArrayView<const GlobalOrdinal>(indicesNew.data(), indicesNew.length()));
225 sparsityPattern->fillComplete();
227 return sparsityPattern;
237 size_t numRows = A.getRowMap()->getLocalNumElements();
238 Teuchos::ArrayRCP<Scalar> diag(numRows);
239 Teuchos::ArrayView<const LocalOrdinal> cols;
240 Teuchos::ArrayView<const Scalar> vals;
241 for (
size_t i = 0; i < numRows; ++i) {
242 A.getLocalRowView(i, cols, vals);
244 for (; j < cols.size(); ++j) {
245 if (Teuchos::as<size_t>(cols[j]) == i) {
250 if (j == cols.size()) {
252 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
265 Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer(
"UtilitiesBase::GetMatrixDiagonalInverse");
267 RCP<const Map> rowMap = A.getRowMap();
268 RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
270 A.getLocalDiagCopy(*diag);
284 Magnitude tol = Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero()),
285 Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(),
286 const bool replaceSingleEntryRowWithZero =
false,
287 const bool useAverageAbsDiagVal =
false) {
289 typedef Teuchos::ScalarTraits<Scalar> TST;
291 RCP<Vector> diag = Teuchos::null;
292 const Scalar zero = TST::zero();
293 const Scalar one = TST::one();
294 const Scalar two = one + one;
296 Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
298 RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
299 Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
300 if(bA == Teuchos::null) {
301 RCP<const Map> rowMap = rcpA->getRowMap();
302 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
303 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
304 Teuchos::Array<Scalar> regSum(diag->getLocalLength());
305 Teuchos::ArrayView<const LocalOrdinal> cols;
306 Teuchos::ArrayView<const Scalar> vals;
308 std::vector<int> nnzPerRow(rowMap->getLocalNumElements());
313 const Magnitude zeroMagn = TST::magnitude(zero);
314 Magnitude avgAbsDiagVal = TST::magnitude(zero);
315 int numDiagsEqualToOne = 0;
316 for (
size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
318 rcpA->getLocalRowView(i, cols, vals);
321 regSum[i] += vals[j];
322 const Magnitude rowEntryMagn = TST::magnitude(vals[j]);
323 if (rowEntryMagn > zeroMagn)
325 diagVals[i] += rowEntryMagn;
326 if (
static_cast<size_t>(cols[j]) == i)
327 avgAbsDiagVal += rowEntryMagn;
329 if (nnzPerRow[i] == 1 && TST::magnitude(diagVals[i])==1.)
330 numDiagsEqualToOne++;
332 if (useAverageAbsDiagVal)
333 tol = TST::magnitude(100 * Teuchos::ScalarTraits<Scalar>::eps()) * (avgAbsDiagVal-numDiagsEqualToOne) / (rowMap->getLocalNumElements()-numDiagsEqualToOne);
335 for (
size_t i = 0; i < rowMap->getLocalNumElements(); ++i) {
336 if (replaceSingleEntryRowWithZero && nnzPerRow[i] <=
static_cast<int>(1))
338 else if ((diagVals[i] != zero) && (TST::magnitude(diagVals[i]) < TST::magnitude(two*regSum[i])))
339 diagVals[i] = one / TST::magnitude((two*regSum[i]));
341 if(TST::magnitude(diagVals[i]) > tol)
342 diagVals[i] = one / diagVals[i];
344 diagVals[i] = valReplacement;
350 TEUCHOS_TEST_FOR_EXCEPTION(doReciprocal, Xpetra::Exceptions::RuntimeError,
351 "UtilitiesBase::GetLumpedMatrixDiagonal(): extracting reciprocal of diagonal of a blocked matrix is not supported");
352 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),
true);
354 for (
size_t row = 0; row < bA->Rows(); ++row) {
355 for (
size_t col = 0; col < bA->Cols(); ++col) {
356 if (!bA->getMatrix(row,col).is_null()) {
358 bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA->getMatrix(row,col)) == Teuchos::null);
359 RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag,row,bThyraMode);
361 ddtemp->update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
362 bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
379 size_t numRows = A.getRowMap()->getLocalNumElements();
380 Magnitude ZERO = Teuchos::ScalarTraits<Magnitude>::zero();
381 Teuchos::ArrayRCP<Magnitude> maxvec(numRows);
382 Teuchos::ArrayView<const LocalOrdinal> cols;
383 Teuchos::ArrayView<const Scalar> vals;
384 for (
size_t i = 0; i < numRows; ++i) {
385 A.getLocalRowView(i, cols, vals);
388 if (Teuchos::as<size_t>(cols[j]) != i) {
389 mymax = std::max(mymax,-Teuchos::ScalarTraits<Scalar>::real(vals[j]));
397 static Teuchos::ArrayRCP<Magnitude>
GetMatrixMaxMinusOffDiagonal(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> &BlockNumber) {
398 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()),std::runtime_error,
"GetMatrixMaxMinusOffDiagonal: BlockNumber must match's A's column map.");
400 Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
402 size_t numRows = A.getRowMap()->getLocalNumElements();
403 Magnitude ZERO = Teuchos::ScalarTraits<Magnitude>::zero();
404 Teuchos::ArrayRCP<Magnitude> maxvec(numRows);
405 Teuchos::ArrayView<const LocalOrdinal> cols;
406 Teuchos::ArrayView<const Scalar> vals;
407 for (
size_t i = 0; i < numRows; ++i) {
408 A.getLocalRowView(i, cols, vals);
411 if (Teuchos::as<size_t>(cols[j]) != i && block_id[i] == block_id[cols[j]]) {
412 mymax = std::max(mymax,-Teuchos::ScalarTraits<Scalar>::real(vals[j]));
429 static Teuchos::RCP<Vector>
GetInverse(Teuchos::RCP<const Vector> v,
Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100,
Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
431 RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),
true);
434 RCP<const BlockedVector> bv = Teuchos::rcp_dynamic_cast<const BlockedVector>(v);
435 if(bv.is_null() ==
false) {
436 RCP<BlockedVector> bret = Teuchos::rcp_dynamic_cast<BlockedVector>(ret);
437 TEUCHOS_TEST_FOR_EXCEPTION(bret.is_null() ==
true,
MueLu::Exceptions::RuntimeError,
"MueLu::UtilitiesBase::GetInverse: return vector should be of type BlockedVector");
438 RCP<const BlockedMap> bmap = bv->getBlockedMap();
439 for(
size_t r = 0; r < bmap->getNumMaps(); ++r) {
440 RCP<const MultiVector> submvec = bv->getMultiVector(r,bmap->getThyraMode());
441 RCP<const Vector> subvec = submvec->getVector(0);
443 bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
449 ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
450 ArrayRCP<const Scalar> inputVals = v->getData(0);
451 for (
size_t i = 0; i < v->getMap()->getLocalNumElements(); ++i) {
452 if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
453 retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
455 retVals[i] = valReplacement;
468 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
471 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
472 if(!browMap.is_null()) rowMap = browMap->getMap();
474 RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
480 Teuchos::ArrayRCP<size_t> offsets;
481 crsOp->getLocalDiagOffsets(offsets);
482 crsOp->getLocalDiagCopy(*localDiag,offsets());
485 ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
488 localDiagVals[i] = diagVals[i];
489 localDiagVals = diagVals = null;
492 RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
493 RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
494 importer = A.getCrsGraph()->getImporter();
495 if (importer == Teuchos::null) {
496 importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
498 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
514 using STS =
typename Teuchos::ScalarTraits<SC>;
517 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
518 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
519 if(!browMap.is_null()) rowMap = browMap->getMap();
521 RCP<Vector> local = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(rowMap);
522 RCP<Vector> ghosted = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(colMap,
true);
523 ArrayRCP<SC> localVals = local->getDataNonConst(0);
525 for (LO row = 0; row < static_cast<LO>(A.getRowMap()->getLocalNumElements()); ++row) {
526 size_t nnz = A.getNumEntriesInLocalRow(row);
527 ArrayView<const LO> indices;
528 ArrayView<const SC> vals;
529 A.getLocalRowView(row, indices, vals);
533 for (LO colID = 0; colID < static_cast<LO>(nnz); colID++) {
534 if(indices[colID] != row) {
541 RCP< const Xpetra::Import<LO,GO,Node> > importer;
542 importer = A.getCrsGraph()->getImporter();
543 if (importer == Teuchos::null) {
544 importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
546 ghosted->doImport(*local, *(importer), Xpetra::INSERT);
552 static RCP<Xpetra::Vector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> >
554 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
555 using STS =
typename Teuchos::ScalarTraits<Scalar>;
556 using MTS =
typename Teuchos::ScalarTraits<Magnitude>;
561 using RealValuedVector = Xpetra::Vector<MT,LO,GO,Node>;
564 RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
565 if(!browMap.is_null()) rowMap = browMap->getMap();
567 RCP<RealValuedVector> local = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(rowMap);
568 RCP<RealValuedVector> ghosted = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(colMap,
true);
569 ArrayRCP<MT> localVals = local->getDataNonConst(0);
571 for (LO rowIdx = 0; rowIdx < static_cast<LO>(A.getRowMap()->getLocalNumElements()); ++rowIdx) {
572 size_t nnz = A.getNumEntriesInLocalRow(rowIdx);
573 ArrayView<const LO> indices;
574 ArrayView<const SC> vals;
575 A.getLocalRowView(rowIdx, indices, vals);
579 for (LO colID = 0; colID < static_cast<LO>(nnz); ++colID) {
580 if(indices[colID] != rowIdx) {
581 si += STS::magnitude(vals[colID]);
584 localVals[rowIdx] = si;
587 RCP< const Xpetra::Import<LO,GO,Node> > importer;
588 importer = A.getCrsGraph()->getImporter();
589 if (importer == Teuchos::null) {
590 importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
592 ghosted->doImport(*local, *(importer), Xpetra::INSERT);
603 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
604 const size_t numVecs = X.getNumVectors();
605 RCP<MultiVector> RES =
Residual(Op, X, RHS);
606 Teuchos::Array<Magnitude> norms(numVecs);
612 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
613 const size_t numVecs = X.getNumVectors();
615 Teuchos::Array<Magnitude> norms(numVecs);
621 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
622 const size_t numVecs = X.getNumVectors();
624 RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(RHS.getMap(), numVecs,
false);
625 Op.residual(X,RHS,*RES);
631 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides");
632 TEUCHOS_TEST_FOR_EXCEPTION(Resid.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of residual vectors != number of right-hand sides");
633 Op.residual(X,RHS,Resid);
649 LocalOrdinal niters = 10,
Magnitude tolerance = 1e-2,
bool verbose =
false,
unsigned int seed = 123) {
651 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
654 RCP<Vector> diagInvVec;
656 RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
657 A.getLocalDiagCopy(*diagVec);
658 diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
659 diagInvVec->reciprocal(*diagVec);
678 LocalOrdinal niters = 10,
Magnitude tolerance = 1e-2,
bool verbose =
false,
unsigned int seed = 123) {
680 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
683 RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
684 RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
685 RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
690 Teuchos::Array<Magnitude> norms(1);
692 typedef Teuchos::ScalarTraits<Scalar> STS;
694 const Scalar zero = STS::zero(), one = STS::one();
697 Magnitude residual = STS::magnitude(zero);
700 for (
int iter = 0; iter < niters; ++iter) {
702 q->update(one/norms[0], *z, zero);
704 if (diagInvVec != Teuchos::null)
705 z->elementWiseMultiply(one, *diagInvVec, *z, zero);
708 if (iter % 100 == 0 || iter + 1 == niters) {
709 r->update(1.0, *z, -lambda, *q, zero);
711 residual = STS::magnitude(norms[0] / lambda);
713 std::cout <<
"Iter = " << iter
714 <<
" Lambda = " << lambda
715 <<
" Residual of A*q - lambda*q = " << residual
725 static RCP<Teuchos::FancyOStream>
MakeFancy(std::ostream& os) {
726 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
735 const size_t numVectors = v.size();
737 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
738 for (
size_t j = 0; j < numVectors; j++) {
739 d += (v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
741 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
748 static typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Distance2(
const Teuchos::ArrayView<double> & weight,
const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>> &v,
LocalOrdinal i0,
LocalOrdinal i1) {
749 const size_t numVectors = v.size();
750 using MT =
typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
752 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
753 for (
size_t j = 0; j < numVectors; j++) {
754 d += Teuchos::as<MT>(weight[j])*(v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
756 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
772 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) {
774 typedef Teuchos::ScalarTraits<Scalar> STS;
775 ArrayRCP<bool> boundaryNodes(numRows,
true);
776 if (count_twos_as_dirichlet) {
778 ArrayView<const LocalOrdinal> indices;
779 ArrayView<const Scalar> vals;
780 A.getLocalRowView(row, indices, vals);
781 size_t nnz = A.getNumEntriesInLocalRow(row);
784 for (col = 0; col < nnz; col++)
785 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
786 if (!boundaryNodes[row])
788 boundaryNodes[row] =
false;
791 boundaryNodes[row] =
true;
796 ArrayView<const LocalOrdinal> indices;
797 ArrayView<const Scalar> vals;
798 A.getLocalRowView(row, indices, vals);
799 size_t nnz = A.getNumEntriesInLocalRow(row);
801 for (
size_t col = 0; col < nnz; col++)
802 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
803 boundaryNodes[row] =
false;
808 return boundaryNodes;
824 static Teuchos::ArrayRCP<const bool>
DetectDirichletRowsExt(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
bool & bHasZeroDiagonal,
const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
827 bHasZeroDiagonal =
false;
829 Teuchos::RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
830 A.getLocalDiagCopy(*diagVec);
831 Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
834 typedef Teuchos::ScalarTraits<Scalar> STS;
835 ArrayRCP<bool> boundaryNodes(numRows,
false);
837 ArrayView<const LocalOrdinal> indices;
838 ArrayView<const Scalar> vals;
839 A.getLocalRowView(row, indices, vals);
841 bool bHasDiag =
false;
842 for (
decltype(indices.size()) col = 0; col < indices.size(); col++) {
843 if ( indices[col] != row) {
844 if (STS::magnitude(vals[col] / STS::magnitude(sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col]))) ) > tol) {
847 }
else bHasDiag =
true;
849 if (bHasDiag ==
false) bHasZeroDiagonal =
true;
850 else if(nnz == 0) boundaryNodes[row] =
true;
852 return boundaryNodes;
863 Teuchos::ArrayRCP<bool> nonzeros) {
864 TEUCHOS_ASSERT(vals.size() == nonzeros.size());
865 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
866 const magnitudeType eps = 2.0*Teuchos::ScalarTraits<magnitudeType>::eps();
867 for(
size_t i=0; i<static_cast<size_t>(vals.size()); i++) {
868 nonzeros[i] = (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > eps);
881 const Teuchos::ArrayRCP<bool>& dirichletRows,
882 Teuchos::ArrayRCP<bool> dirichletCols,
883 Teuchos::ArrayRCP<bool> dirichletDomain) {
884 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
885 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A .getDomainMap();
886 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
887 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
888 TEUCHOS_ASSERT(
static_cast<size_t>(dirichletRows.size()) == rowMap->getLocalNumElements());
889 TEUCHOS_ASSERT(
static_cast<size_t>(dirichletCols.size()) == colMap->getLocalNumElements());
890 TEUCHOS_ASSERT(
static_cast<size_t>(dirichletDomain.size()) == domMap->getLocalNumElements());
891 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, 1,
true);
893 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
894 if (dirichletRows[i]) {
895 ArrayView<const LocalOrdinal> indices;
896 ArrayView<const Scalar> values;
897 A.getLocalRowView(i,indices,values);
898 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
899 myColsToZero->replaceLocalValue(indices[j],0,one);
903 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
904 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
905 if (!importer.is_null()) {
906 globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, 1,
true);
908 globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
910 myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
913 globalColsToZero = myColsToZero;
915 FindNonZeros(globalColsToZero->getData(0),dirichletDomain);
932 static void ApplyRowSumCriterion(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
933 typedef Teuchos::ScalarTraits<Scalar> STS;
934 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
935 typedef Teuchos::ScalarTraits<MT> MTS;
936 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> rowmap = A.getRowMap();
937 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
938 size_t nnz = A.getNumEntriesInLocalRow(row);
939 ArrayView<const LocalOrdinal> indices;
940 ArrayView<const Scalar> vals;
941 A.getLocalRowView(row, indices, vals);
943 Scalar rowsum = STS::zero();
944 Scalar diagval = STS::zero();
946 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
949 diagval = vals[colID];
950 rowsum += vals[colID];
953 if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
955 dirichletRows[row] =
true;
960 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) {
961 typedef Teuchos::ScalarTraits<Scalar> STS;
962 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
963 typedef Teuchos::ScalarTraits<MT> MTS;
964 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowmap = A.getRowMap();
966 TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()),std::runtime_error,
"ApplyRowSumCriterion: BlockNumber must match's A's column map.");
968 Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
969 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getLocalNumElements()); ++row) {
970 size_t nnz = A.getNumEntriesInLocalRow(row);
971 ArrayView<const LocalOrdinal> indices;
972 ArrayView<const Scalar> vals;
973 A.getLocalRowView(row, indices, vals);
975 Scalar rowsum = STS::zero();
976 Scalar diagval = STS::zero();
977 for (
LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
980 diagval = vals[colID];
981 if(block_id[row] == block_id[col])
982 rowsum += vals[colID];
986 if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
988 dirichletRows[row] =
true;
1005 static Teuchos::ArrayRCP<const bool>
DetectDirichletCols(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
1006 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
1007 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1008 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1009 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
1010 Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
1011 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,1);
1012 myColsToZero->putScalar(zero);
1014 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1015 if (dirichletRows[i]) {
1016 Teuchos::ArrayView<const LocalOrdinal> indices;
1017 Teuchos::ArrayView<const Scalar> values;
1018 A.getLocalRowView(i,indices,values);
1019 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
1020 myColsToZero->replaceLocalValue(indices[j],0,one);
1024 Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,1);
1025 globalColsToZero->putScalar(zero);
1026 Teuchos::RCP<Xpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
1028 globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
1030 myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
1031 Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
1032 Teuchos::ArrayRCP<bool> dirichletCols(colMap->getLocalNumElements(),
true);
1033 Magnitude eps = Teuchos::ScalarTraits<Magnitude>::eps();
1034 for(
size_t i=0; i<colMap->getLocalNumElements(); i++) {
1035 dirichletCols[i] = Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i])>2.0*eps;
1037 return dirichletCols;
1045 static Scalar Frobenius(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
1050 TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()),
Exceptions::Incompatible,
"MueLu::CGSolver::Frobenius: row maps are incompatible");
1051 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(),
Exceptions::RuntimeError,
"Matrices must be fill completed");
1053 const Map& AColMap = *A.getColMap();
1054 const Map& BColMap = *B.getColMap();
1056 Teuchos::ArrayView<const LocalOrdinal> indA, indB;
1057 Teuchos::ArrayView<const Scalar> valA, valB;
1058 size_t nnzA = 0, nnzB = 0;
1070 Teuchos::Array<Scalar> valBAll(BColMap.getLocalNumElements());
1072 LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
1073 Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
1074 size_t numRows = A.getLocalNumRows();
1075 for (
size_t i = 0; i < numRows; i++) {
1076 A.getLocalRowView(i, indA, valA);
1077 B.getLocalRowView(i, indB, valB);
1082 for (
size_t j = 0; j < nnzB; j++)
1083 valBAll[indB[j]] = valB[j];
1085 for (
size_t j = 0; j < nnzA; j++) {
1088 LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
1090 f += valBAll[ind] * valA[j];
1094 for (
size_t j = 0; j < nnzB; j++)
1095 valBAll[indB[j]] = zero;
1118 int maxint = INT_MAX;
1119 int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.getRank()+1)/(comm.getSize()+one)) );
1120 if (mySeed < 1 || mySeed == maxint) {
1121 std::ostringstream errStr;
1122 errStr <<
"Error detected with random seed = " << mySeed <<
". It should be in the interval [1,2^31-2].";
1127 Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
1139 static void FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
1140 std::vector<LocalOrdinal>& dirichletRows,
bool count_twos_as_dirichlet=
false) {
1141 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
1142 dirichletRows.resize(0);
1143 for(
size_t i=0; i<A->getLocalNumRows(); i++) {
1144 Teuchos::ArrayView<const LocalOrdinal> indices;
1145 Teuchos::ArrayView<const Scalar> values;
1146 A->getLocalRowView(i,indices,values);
1148 for (
size_t j=0; j<(size_t)indices.size(); j++) {
1149 if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > Teuchos::ScalarTraits<MT>::eps()) {
1153 if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
1154 dirichletRows.push_back(i);
1162 const std::vector<LocalOrdinal>& dirichletRows) {
1163 RCP<const Map> Rmap = A->getRowMap();
1164 RCP<const Map> Cmap = A->getColMap();
1165 Scalar one =Teuchos::ScalarTraits<Scalar>::one();
1166 Scalar zero =Teuchos::ScalarTraits<Scalar>::zero();
1168 for(
size_t i=0; i<dirichletRows.size(); i++) {
1169 GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
1171 Teuchos::ArrayView<const LocalOrdinal> indices;
1172 Teuchos::ArrayView<const Scalar> values;
1173 A->getLocalRowView(dirichletRows[i],indices,values);
1175 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
1176 for(
size_t j=0; j<(size_t)indices.size(); j++) {
1177 if(Cmap->getGlobalElement(indices[j])==row_gid)
1188 const Teuchos::ArrayRCP<const bool>& dirichletRows) {
1189 TEUCHOS_ASSERT(A->isFillComplete());
1190 RCP<const Map> domMap = A->getDomainMap();
1191 RCP<const Map> ranMap = A->getRangeMap();
1192 RCP<const Map> Rmap = A->getRowMap();
1193 RCP<const Map> Cmap = A->getColMap();
1194 TEUCHOS_ASSERT(
static_cast<size_t>(dirichletRows.size()) == Rmap->getLocalNumElements());
1195 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1196 const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1198 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1199 if (dirichletRows[i]){
1202 Teuchos::ArrayView<const LocalOrdinal> indices;
1203 Teuchos::ArrayView<const Scalar> values;
1204 A->getLocalRowView(i,indices,values);
1206 Teuchos::ArrayRCP<Scalar> valuesNC(values.size());
1207 for(
size_t j=0; j<(size_t)indices.size(); j++) {
1208 if(Cmap->getGlobalElement(indices[j])==row_gid)
1213 A->replaceLocalValues(i,indices,valuesNC());
1216 A->fillComplete(domMap, ranMap);
1221 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1222 const std::vector<LocalOrdinal>& dirichletRows,
1223 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1224 for(
size_t i=0; i<dirichletRows.size(); i++) {
1225 Teuchos::ArrayView<const LocalOrdinal> indices;
1226 Teuchos::ArrayView<const Scalar> values;
1227 A->getLocalRowView(dirichletRows[i],indices,values);
1229 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
1230 for(
size_t j=0; j<(size_t)indices.size(); j++)
1231 valuesNC[j]=replaceWith;
1237 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1238 const Teuchos::ArrayRCP<const bool>& dirichletRows,
1239 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1240 TEUCHOS_ASSERT(
static_cast<size_t>(dirichletRows.size()) == A->getRowMap()->getLocalNumElements());
1241 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1242 if (dirichletRows[i]) {
1243 Teuchos::ArrayView<const LocalOrdinal> indices;
1244 Teuchos::ArrayView<const Scalar> values;
1245 A->getLocalRowView(i,indices,values);
1247 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
1248 for(
size_t j=0; j<(size_t)indices.size(); j++)
1249 valuesNC[j]=replaceWith;
1256 static void ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& X,
1257 const Teuchos::ArrayRCP<const bool>& dirichletRows,
1258 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1259 TEUCHOS_ASSERT(
static_cast<size_t>(dirichletRows.size()) == X->getMap()->getLocalNumElements());
1260 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1261 if (dirichletRows[i]) {
1262 for(
size_t j=0; j<X->getNumVectors(); j++)
1263 X->replaceLocalValue(i,j,replaceWith);
1271 const Teuchos::ArrayRCP<const bool>& dirichletCols,
1272 Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1273 TEUCHOS_ASSERT(
static_cast<size_t>(dirichletCols.size()) == A->getColMap()->getLocalNumElements());
1274 for(
size_t i=0; i<A->getLocalNumRows(); i++) {
1275 Teuchos::ArrayView<const LocalOrdinal> indices;
1276 Teuchos::ArrayView<const Scalar> values;
1277 A->getLocalRowView(i,indices,values);
1279 Scalar* valuesNC =
const_cast<Scalar*
>(values.getRawPtr());
1280 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
1281 if (dirichletCols[indices[j]])
1282 valuesNC[j] = replaceWith;
1288 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletRow,
1289 Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletCol) {
1292 if(!A->getRowMap()->isSameAs(*A->getDomainMap())) {
1293 throw std::runtime_error(
"UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
1295 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A->getCrsGraph()->getImporter();
1296 bool has_import = !importer.is_null();
1299 std::vector<LocalOrdinal> dirichletRows;
1303 printf(
"[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
1304 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++)
1305 printf(
"%d ",dirichletRows[i]);
1310 isDirichletRow = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getRowMap(),
true);
1311 isDirichletCol = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getColMap(),
true);
1314 Teuchos::ArrayRCP<int> dr_rcp = isDirichletRow->getDataNonConst(0);
1315 Teuchos::ArrayView<int> dr = dr_rcp();
1316 Teuchos::ArrayRCP<int> dc_rcp = isDirichletCol->getDataNonConst(0);
1317 Teuchos::ArrayView<int> dc = dc_rcp();
1318 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1319 dr[dirichletRows[i]] = 1;
1320 if(!has_import) dc[dirichletRows[i]] = 1;
1325 isDirichletCol->doImport(*isDirichletRow,*importer,Xpetra::CombineMode::ADD);
1332 static RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> >
GeneratedBlockedTargetMap(
const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> & sourceBlockedMap,
1333 const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> & Importer) {
1334 typedef Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> IntVector;
1335 Xpetra::UnderlyingLib lib = sourceBlockedMap.lib();
1338 RCP<const Map> fullMap = sourceBlockedMap.getMap();
1339 RCP<const Map> stridedMap = Teuchos::rcp_dynamic_cast<const Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >(fullMap);
1340 if(!stridedMap.is_null()) fullMap = stridedMap->getMap();
1343 const size_t numSubMaps = sourceBlockedMap.getNumMaps();
1344 if(!Importer.getSourceMap()->isCompatible(*fullMap))
1345 throw std::runtime_error(
"GenerateBlockedTargetMap(): Map compatibility error");
1348 RCP<IntVector> block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(fullMap);
1350 for(
size_t i=0; i<numSubMaps; i++) {
1351 RCP<const Map> map = sourceBlockedMap.getMap(i);
1353 for(
size_t j=0; j<map->getLocalNumElements(); j++) {
1354 LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
1355 block_ids->replaceLocalValue(jj,(
int)i);
1360 RCP<const Map> targetMap = Importer.getTargetMap();
1361 RCP<IntVector> new_block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(targetMap);
1362 new_block_ids->doImport(*block_ids,Importer,Xpetra::CombineMode::ADD);
1363 Teuchos::ArrayRCP<const int> dataRCP = new_block_ids->getData(0);
1364 Teuchos::ArrayView<const int> data = dataRCP();
1368 Teuchos::Array<Teuchos::Array<GlobalOrdinal> > elementsInSubMap(numSubMaps);
1369 for(
size_t i=0; i<targetMap->getLocalNumElements(); i++) {
1370 elementsInSubMap[data[i]].push_back(targetMap->getGlobalElement(i));
1374 std::vector<RCP<const Map> > subMaps(numSubMaps);
1375 for(
size_t i=0; i<numSubMaps; i++) {
1376 subMaps[i] = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(lib,Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),elementsInSubMap[i](),targetMap->getIndexBase(),targetMap->getComm());
1380 return rcp(
new BlockedMap(targetMap,subMaps));
1386 static bool MapsAreNested(
const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& rowMap,
const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& colMap) {
1387 ArrayView<const GlobalOrdinal> rowElements = rowMap.getLocalElementList();
1388 ArrayView<const GlobalOrdinal> colElements = colMap.getLocalElementList();
1390 const size_t numElements = rowElements.size();
1392 if (
size_t(colElements.size()) < numElements)
1395 bool goodMap =
true;
1396 for (
size_t i = 0; i < numElements; i++)
1397 if (rowElements[i] != colElements[i]) {
1415#define MUELU_UTILITIESBASE_SHORT
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Exception throws to report incompatible objects (like maps).
Exception throws to report errors in the internal logical of the program.
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MultiVector
static RCP< Vector > GetMatrixOverlappedDeletedRowsum(const Matrix &A)
Extract Overlapped Matrix Deleted Rowsum.
Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > Map
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 RCP< CrsMatrixWrap > GetThresholdedMatrix(const RCP< Matrix > &Ain, const Scalar threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a matrix.
static void ZeroDirichletCols(Teuchos::RCP< Matrix > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
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 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.
Teuchos::ScalarTraits< Scalar >::magnitudeType Magnitude
Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrix
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.
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
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)
Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Vector
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
Xpetra::BlockedVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedVector
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 Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i\not=k}(-a_ik), for each for i in the matrix.
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Extract Matrix Diagonal.
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.
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 ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows)
static RCP< Matrix > Crs2Op(RCP< CrsMatrix > Op)
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 FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
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::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
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 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 bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
Xpetra::CrsGraphFactory< LocalOrdinal, GlobalOrdinal, Node > CrsGraphFactory
static void Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS, MultiVector &Resid)
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.
Xpetra::CrsMatrixWrap< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrixWrap
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 RCP< Xpetra::Vector< Magnitude, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedAbsDeletedRowsum(const Matrix &A)
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS, MultiVector &Resid)
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &BlockNumber)
Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedMultiVector
Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > CrsGraph
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > BlockedMap
static RCP< Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > GetThresholdedGraph(const RCP< Matrix > &A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a graph.
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
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.
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 void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
Namespace for MueLu classes and methods.
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode