46 #ifndef MUELU_UTILITIESBASE_DECL_HPP 47 #define MUELU_UTILITIESBASE_DECL_HPP 54 #include <Teuchos_DefaultComm.hpp> 55 #include <Teuchos_ScalarTraits.hpp> 56 #include <Teuchos_ParameterList.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 115 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Magnitude;
118 static RCP<Matrix>
Crs2Op(RCP<CrsMatrix> Op) {
120 return Teuchos::null;
131 size_t numRows = A.
getRowMap()->getNodeNumElements();
132 Teuchos::ArrayRCP<Scalar> diag(numRows);
133 Teuchos::ArrayView<const LocalOrdinal> cols;
134 Teuchos::ArrayView<const Scalar> vals;
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()) {
146 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
160 RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
162 RCP<const Map> rowMap = rcpA->getRowMap();
165 rcpA->getLocalDiagCopy(*diag);
182 size_t numRows = A.
getRowMap()->getNodeNumElements();
183 Teuchos::ArrayRCP<Scalar> diag(numRows);
184 Teuchos::ArrayView<const LocalOrdinal> cols;
185 Teuchos::ArrayView<const Scalar> vals;
186 for (
size_t i = 0; i < numRows; ++i) {
188 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
189 for (LocalOrdinal j = 0; j < cols.size(); ++j) {
190 diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
204 RCP<Vector> diag = Teuchos::null;
206 RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
208 if(bA == Teuchos::null) {
209 RCP<const Map> rowMap = rcpA->
getRowMap();
211 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
212 Teuchos::ArrayView<const LocalOrdinal> cols;
213 Teuchos::ArrayView<const Scalar> vals;
214 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
215 rcpA->getLocalRowView(i, cols, vals);
216 diagVals[i] = Teuchos::ScalarTraits<Scalar>::zero();
217 for (LocalOrdinal j = 0; j < cols.size(); ++j) {
218 diagVals[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[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);
254 static Teuchos::RCP<Vector>
GetInverse(Teuchos::RCP<const Vector> v,
Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar tolReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
257 ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
258 ArrayRCP<const Scalar> inputVals = v->getData(0);
259 for (
size_t i = 0; i < v->getMap()->getNodeNumElements(); ++i) {
260 if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
261 retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
263 retVals[i] = tolReplacement;
284 Teuchos::ArrayRCP<size_t> offsets;
289 ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
291 for (LocalOrdinal i = 0; i < localDiagVals.size(); i++)
292 localDiagVals[i] = diagVals[i];
293 localDiagVals = diagVals = null;
297 RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
299 if (importer == Teuchos::null) {
314 RCP<MultiVector> RES =
Residual(Op, X, RHS);
315 Teuchos::Array<Magnitude> norms(numVecs);
323 Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
325 Op.
apply(X, *RES, Teuchos::NO_TRANS, one, zero);
326 RES->update(one, RHS, negone);
333 RCP<const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
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.");
387 Teuchos::Array<Magnitude> norms(1);
389 typedef Teuchos::ScalarTraits<Scalar> STS;
391 const Scalar zero = STS::zero(), one = STS::one();
393 Scalar lambda = zero;
394 Magnitude residual = STS::magnitude(zero);
397 RCP<Vector> diagInvVec;
402 diagInvVec->reciprocal(*diagVec);
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)
432 static RCP<Teuchos::FancyOStream>
MakeFancy(std::ostream& os) {
433 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
444 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
445 for (
size_t j = 0; j < numVectors; j++) {
446 Teuchos::ArrayRCP<const Scalar> vv = v.
getData(j);
447 d += (vv[i0] - vv[i1])*(vv[i0] - vv[i1]);
449 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
461 typedef Teuchos::ScalarTraits<Scalar> STS;
462 ArrayRCP<bool> boundaryNodes(numRows,
true);
463 for (LocalOrdinal row = 0; row < numRows; row++) {
464 ArrayView<const LocalOrdinal> indices;
465 ArrayView<const Scalar> vals;
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;
493 Teuchos::ArrayView<const LocalOrdinal> indA, indB;
494 Teuchos::ArrayView<const Scalar> valA, valB;
495 size_t nnzA = 0, nnzB = 0;
509 LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
510 Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
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].";
564 Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
580 #define MUELU_UTILITIESBASE_SHORT 581 #endif // MUELU_UTILITIESBASE_DECL_HPP 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 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
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.
Exception throws to report incompatible objects (like maps).
static void PauseForDebugger()
#define MueLu_sumAll(rcpComm, in, out)
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
virtual bool isFillComplete() 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
Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Vector
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
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)
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.