46 #ifndef MUELU_UTILITIESBASE_DECL_HPP 47 #define MUELU_UTILITIESBASE_DECL_HPP 54 #include <Teuchos_DefaultComm.hpp> 85 #define MueLu_sumAll(rcpComm, in, out) \ 86 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out)) 87 #define MueLu_minAll(rcpComm, in, out) \ 88 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out)) 89 #define MueLu_maxAll(rcpComm, in, out) \ 90 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out)) 99 template <
class Scalar,
100 class LocalOrdinal = int,
101 class GlobalOrdinal = LocalOrdinal,
103 class UtilitiesBase {
105 #undef MUELU_UTILITIESBASE_SHORT 120 return Teuchos::null;
135 for (
size_t i = 0; i < numRows; ++i) {
138 for (; j < cols.
size(); ++j) {
139 if (Teuchos::as<size_t>(cols[j]) == i) {
144 if (j == cols.
size()) {
186 for (
size_t i = 0; i < numRows; ++i) {
189 for (LocalOrdinal j = 0; j < cols.
size(); ++j) {
208 if(bA == Teuchos::null) {
217 for (LocalOrdinal j = 0; j < cols.
size(); ++j) {
228 for (
size_t row = 0; row < bA->Rows(); ++row) {
229 for (
size_t col = 0; col < bA->Cols(); ++col) {
230 if (!bA->getMatrix(row,col).
is_null()) {
233 RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag,row,bThyraMode);
235 ddtemp->
update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
236 bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
259 for (
size_t i = 0; i < v->
getMap()->getNodeNumElements(); ++i) {
263 retVals[i] = tolReplacement;
291 for (LocalOrdinal i = 0; i < localDiagVals.
size(); i++)
292 localDiagVals[i] = diagVals[i];
293 localDiagVals = diagVals = null;
299 if (importer == Teuchos::null) {
326 RES->
update(one, RHS, negone);
334 int myPID = comm->getRank();
337 for (
int i = 0; i <comm->getSize(); i++) {
339 gethostname(hostname,
sizeof(hostname));
340 std::cout <<
"Host: " << hostname <<
"\tMPI rank: " << myPID <<
",\tPID: " << pid <<
"\n\tattach " << pid << std::endl;
345 std::cout <<
"** Enter a character to continue > " << std::endl;
347 int r = scanf(
"%c", &go);
375 LocalOrdinal niters = 10,
Magnitude tolerance = 1e-2,
bool verbose =
false,
unsigned int seed = 123) {
377 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
391 const Scalar zero = STS::zero(), one = STS::one();
393 Scalar lambda = zero;
394 Magnitude residual = STS::magnitude(zero);
405 for (
int iter = 0; iter < niters; ++iter) {
407 q->
update(one/norms[0], *z, zero);
410 z->elementWiseMultiply(one, *diagInvVec, *z, zero);
413 if (iter % 100 == 0 || iter + 1 == niters) {
414 r->update(1.0, *z, -lambda, *q, zero);
416 residual = STS::magnitude(norms[0] / lambda);
418 std::cout <<
"Iter = " << iter
419 <<
" Lambda = " << lambda
420 <<
" Residual of A*q - lambda*q = " << residual
424 if (residual < tolerance)
445 for (
size_t j = 0; j < numVectors; j++) {
447 d += (vv[i0] - vv[i1])*(vv[i0] - vv[i1]);
463 for (LocalOrdinal row = 0; row < numRows; row++) {
469 for (
size_t col = 0; col < nnz; col++)
470 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
471 boundaryNodes[row] =
false;
475 return boundaryNodes;
495 size_t nnzA = 0, nnzB = 0;
512 for (
size_t i = 0; i < numRows; i++) {
519 for (
size_t j = 0; j < nnzB; j++)
520 valBAll[indB[j]] = valB[j];
522 for (
size_t j = 0; j < nnzA; j++) {
525 LocalOrdinal ind = BColMap.
getLocalElement(AColMap.getGlobalElement(indA[j]));
527 f += valBAll[ind] * valA[j];
531 for (
size_t j = 0; j < nnzB; j++)
532 valBAll[indB[j]] = zero;
555 int maxint = INT_MAX;
556 int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.
getRank()+1)/(comm.
getSize()+one)) );
557 if (mySeed < 1 || mySeed == maxint) {
558 std::ostringstream errStr;
559 errStr <<
"Error detected with random seed = " << mySeed <<
". It should be in the interval [1,2^31-2].";
580 #define MUELU_UTILITIESBASE_SHORT 581 #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)
virtual void norm2(const Teuchos::ArrayView< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &norms) const=0
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.
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MultiVector
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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 Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
virtual LocalOrdinal getLocalElement(GlobalOrdinal globalIndex) const=0
Namespace for MueLu classes and methods.
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Exception throws to report incompatible objects (like maps).
virtual void reciprocal(const MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)=0
static void PauseForDebugger()
#define MueLu_sumAll(rcpComm, in, out)
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
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)
virtual bool isFillComplete() const=0
virtual int getRank() const=0
virtual Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > getMap() const =0
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
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)
virtual Scalar dot(const Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &a) const=0
static magnitudeType magnitude(T a)
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
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 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)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows.
virtual Teuchos::RCP< const Map > getRangeMap() 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)
Exception throws to report errors in the internal logical of the program.
Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > Map
virtual size_t getNumVectors() const=0
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
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 void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
virtual Teuchos::RCP< const Map > getDomainMap() const=0
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
virtual int getSize() const=0