46 #ifndef MUELU_UTILITIESBASE_DECL_HPP 47 #define MUELU_UTILITIESBASE_DECL_HPP 57 #include <Teuchos_DefaultComm.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,
110 class UtilitiesBase {
112 #undef MUELU_UTILITIESBASE_SHORT 130 return Teuchos::null;
145 for (
size_t i = 0; i < numRows; ++i) {
148 for (; j < cols.
size(); ++j) {
149 if (Teuchos::as<size_t>(cols[j]) == i) {
154 if (j == cols.
size()) {
196 for (
size_t i = 0; i < numRows; ++i) {
199 for (LocalOrdinal j = 0; j < cols.
size(); ++j) {
218 if(bA == Teuchos::null) {
227 for (LocalOrdinal j = 0; j < cols.
size(); ++j) {
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()) {
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);
270 if(bv.is_null() ==
false) {
274 for(
size_t r = 0; r < bmap->getNumMaps(); ++r) {
286 for (
size_t i = 0; i < v->
getMap()->getNodeNumElements(); ++i) {
290 retVals[i] = tolReplacement;
322 for (LocalOrdinal i = 0; i < localDiagVals.
size(); i++)
323 localDiagVals[i] = diagVals[i];
324 localDiagVals = diagVals = null;
330 if (importer == Teuchos::null) {
367 RES->
update(one, RHS, negone);
378 Resid.
update(one, RHS, negone);
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.");
442 const Scalar zero = STS::zero(), one = STS::one();
444 Scalar lambda = zero;
445 Magnitude residual = STS::magnitude(zero);
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)
493 const size_t numVectors = v.size();
496 for (
size_t j = 0; j < numVectors; j++) {
497 d += (v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
518 if (count_twos_as_dirichlet) {
519 for (LocalOrdinal row = 0; row < numRows; 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++) {
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;
569 bHasZeroDiagonal =
false;
578 for (LocalOrdinal row = 0; row < numRows; row++) {
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;
615 myColsToZero->putScalar(zero);
617 for(
size_t i=0; i<(size_t) dirichletRows.
size(); i++) {
618 if (dirichletRows[i]) {
622 for(
size_t j=0; j<static_cast<size_t>(indices.
size()); j++)
623 myColsToZero->replaceLocalValue(indices[j],0,one);
628 globalColsToZero->putScalar(zero);
631 globalColsToZero->doExport(*myColsToZero,*exporter,
Xpetra::ADD);
633 myColsToZero->doImport(*globalColsToZero,*exporter,
Xpetra::INSERT);
637 for(
size_t i=0; i<colMap->getNodeNumElements(); i++) {
640 return dirichletCols;
661 size_t nnzA = 0, nnzB = 0;
678 for (
size_t i = 0; i < numRows; i++) {
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].";
743 std::vector<LocalOrdinal>& dirichletRows,
bool count_twos_as_dirichlet=
false) {
745 dirichletRows.resize(0);
751 for (
size_t j=0; j<(size_t)indices.
size(); j++) {
756 if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
757 dirichletRows.push_back(i);
765 const std::vector<LocalOrdinal>& dirichletRows) {
771 for(
size_t i=0; i<dirichletRows.size(); i++) {
778 Scalar* valuesNC =
const_cast<Scalar*
>(values.
getRawPtr());
779 for(
size_t j=0; j<(size_t)indices.
size(); j++) {
797 for(
size_t i=0; i<(size_t) dirichletRows.
size(); i++) {
798 if (dirichletRows[i]){
805 Scalar* valuesNC =
const_cast<Scalar*
>(values.
getRawPtr());
806 for(
size_t j=0; j<(size_t)indices.
size(); j++) {
819 const std::vector<LocalOrdinal>& dirichletRows,
821 for(
size_t i=0; i<dirichletRows.size(); i++) {
826 Scalar* valuesNC =
const_cast<Scalar*
>(values.
getRawPtr());
827 for(
size_t j=0; j<(size_t)indices.
size(); j++)
828 valuesNC[j]=replaceWith;
837 for(
size_t i=0; i<(size_t) dirichletRows.
size(); i++) {
838 if (dirichletRows[i]) {
843 Scalar* valuesNC =
const_cast<Scalar*
>(values.
getRawPtr());
844 for(
size_t j=0; j<(size_t)indices.
size(); j++)
845 valuesNC[j]=replaceWith;
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);
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;
887 throw std::runtime_error(
"UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
890 bool has_import = !importer.
is_null();
893 std::vector<LocalOrdinal> dirichletRows;
898 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++)
899 printf(
"%d ",dirichletRows[i]);
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);
935 const size_t numSubMaps = sourceBlockedMap.
getNumMaps();
937 throw std::runtime_error(
"GenerateBlockedTargetMap(): Map compatibility error");
942 for(
size_t i=0; i<numSubMaps; i++) {
947 block_ids->replaceLocalValue(jj,(
int)i);
954 new_block_ids->doImport(*block_ids,Importer,Xpetra::CombineMode::ADD);
966 std::vector<RCP<const Map> > subMaps(numSubMaps);
967 for(
size_t i=0; i<numSubMaps; i++) {
984 #define MUELU_UTILITIESBASE_SHORT 985 #endif // MUELU_UTILITIESBASE_DECL_HPP virtual Teuchos::ArrayRCP< Scalar > getDataNonConst(size_t j)=0
Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
virtual void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const=0
virtual size_t getNodeNumRows() const=0
static RCP< Import< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
static RCP< Export< LocalOrdinal, GlobalOrdinal, Node > > Build(const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &source, const RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &target)
virtual void norm2(const Teuchos::ArrayView< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &norms) const=0
Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedMultiVector
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm() const=0
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.
size_t getNumMaps() const
static magnitudeType eps()
virtual const RCP< const Map > & getRowMap() const
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.
void setMultiVector(size_t r, Teuchos::RCP< const Vector > v, bool bThyraMode)
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
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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)
virtual size_t getNodeNumElements() const=0
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)
virtual LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const=0
Namespace for MueLu classes and methods.
virtual GlobalOrdinal getIndexBase() const=0
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
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)
virtual void reciprocal(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)=0
static void PauseForDebugger()
#define MueLu_sumAll(rcpComm, in, out)
const RCP< const Map > getMap(size_t i, bool bThyraMode=false) const
virtual void update(const Scalar &alpha, const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Scalar &beta)=0
virtual void doImport(const DistObject< Scalar, LocalOrdinal, GlobalOrdinal, Node > &source, const Import< LocalOrdinal, GlobalOrdinal, Node > &importer, CombineMode CM)=0
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.
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getTargetMap() const=0
virtual const RCP< const Map > & getColMap() const
virtual void getLocalDiagCopy(Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &diag) const=0
virtual void getLocalRowView(LocalOrdinal LocalRow, ArrayView< const LocalOrdinal > &indices, ArrayView< const Scalar > &values) const=0
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)
virtual bool isFillComplete() const=0
virtual int getRank() const=0
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
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.
virtual GlobalOrdinal getGlobalElement(LocalOrdinal localIndex) const=0
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
virtual size_t getNumEntriesInLocalRow(LocalOrdinal localRow) const=0
static void seedrandom(unsigned int s)
Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Vector
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
virtual Scalar dot(const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &a) const=0
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
static magnitudeType magnitude(T a)
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)
void push_back(const value_type &x)
virtual RCP< const CrsGraph > getCrsGraph() const=0
void getLocalDiagOffsets(Teuchos::ArrayRCP< size_t > &offsets) const
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
virtual bool isSameAs(const Map< LocalOrdinal, GlobalOrdinal, Node > &map) const=0
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)
void getLocalDiagCopy(Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &diag) const
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)
virtual Teuchos::RCP< const Map > getRangeMap() const=0
virtual Teuchos::RCP< Node > getNode() const=0
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
virtual RCP< const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const
Exception throws to report errors in the internal logical of the program.
Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > Map
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getSourceMap() const=0
virtual size_t getNumVectors() const=0
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
virtual UnderlyingLib lib() const
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.
virtual Teuchos::ArrayRCP< const Scalar > getData(size_t j) const=0
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed, const Teuchos::RCP< Node > &=Teuchos::null)
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.
virtual Teuchos::RCP< const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getVector(size_t j) const=0
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &X, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
virtual Teuchos::RCP< const Map > getDomainMap() const=0
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
virtual int getSize() const=0