46 #ifndef MUELU_UTILITIESBASE_DECL_HPP 47 #define MUELU_UTILITIESBASE_DECL_HPP 54 #include <Teuchos_DefaultComm.hpp> 58 #include <Xpetra_BlockedCrsMatrix_fwd.hpp> 59 #include <Xpetra_CrsMatrix_fwd.hpp> 60 #include <Xpetra_CrsMatrixWrap_fwd.hpp> 61 #include <Xpetra_Map_fwd.hpp> 62 #include <Xpetra_MapFactory_fwd.hpp> 63 #include <Xpetra_Matrix_fwd.hpp> 64 #include <Xpetra_MatrixFactory_fwd.hpp> 65 #include <Xpetra_MultiVector_fwd.hpp> 66 #include <Xpetra_MultiVectorFactory_fwd.hpp> 67 #include <Xpetra_Operator_fwd.hpp> 68 #include <Xpetra_Vector_fwd.hpp> 69 #include <Xpetra_VectorFactory_fwd.hpp> 70 #include <Xpetra_ExportFactory.hpp> 72 #include <Xpetra_Import.hpp> 73 #include <Xpetra_ImportFactory.hpp> 74 #include <Xpetra_MatrixMatrix.hpp> 75 #include <Xpetra_CrsMatrixWrap.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,
102 class Node = KokkosClassic::DefaultNode::DefaultNodeType>
103 class UtilitiesBase {
105 #undef MUELU_UTILITIESBASE_SHORT 108 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>
CrsMatrixWrap;
109 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
CrsMatrix;
110 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
Matrix;
111 typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
Vector;
112 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
MultiVector;
113 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>
Map;
120 return Teuchos::null;
131 size_t numRows = A.getRowMap()->getNodeNumElements();
135 for (
size_t i = 0; i < numRows; ++i) {
136 A.getLocalRowView(i, cols, vals);
138 for (; j < cols.
size(); ++j) {
139 if (Teuchos::as<size_t>(cols[j]) == i) {
144 if (j == cols.
size()) {
163 RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
165 rcpA->getLocalDiagCopy(*diag);
182 size_t numRows = A.getRowMap()->getNodeNumElements();
186 for (
size_t i = 0; i < numRows; ++i) {
187 A.getLocalRowView(i, cols, vals);
189 for (LocalOrdinal j = 0; j < cols.
size(); ++j) {
207 Teuchos::rcp_dynamic_cast<
const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
208 if(bA == Teuchos::null) {
210 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,
true);
214 for (
size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
215 rcpA->getLocalRowView(i, cols, vals);
217 for (LocalOrdinal j = 0; j < cols.
size(); ++j) {
226 diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),
true);
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()) {
232 bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA->getMatrix(row,col)) == Teuchos::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);
256 RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),
true);
259 for (
size_t i = 0; i < v->getMap()->getNodeNumElements(); ++i) {
263 retVals[i] = tolReplacement;
277 RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
285 crsOp->getLocalDiagOffsets(offsets);
286 crsOp->getLocalDiagCopy(*localDiag,offsets());
291 for (LocalOrdinal i = 0; i < localDiagVals.
size(); i++)
292 localDiagVals[i] = diagVals[i];
293 localDiagVals = diagVals = null;
296 RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
298 importer = A.getCrsGraph()->getImporter();
299 if (importer == Teuchos::null) {
300 importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
302 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
313 const size_t numVecs = X.getNumVectors();
322 const size_t numVecs = X.getNumVectors();
324 RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRangeMap(), numVecs,
false);
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.");
380 RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
381 RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
382 RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
391 const Scalar zero = STS::zero(), one = STS::one();
393 Scalar lambda = zero;
394 Magnitude residual = STS::magnitude(zero);
399 RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
400 A.getLocalDiagCopy(*diagVec);
401 diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
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)
442 size_t numVectors = v.getNumVectors();
445 for (
size_t j = 0; j < numVectors; j++) {
447 d += (vv[i0] - vv[i1])*(vv[i0] - vv[i1]);
460 LocalOrdinal numRows = A.getNodeNumRows();
463 for (LocalOrdinal row = 0; row < numRows; row++) {
466 A.getLocalRowView(row, indices, vals);
467 size_t nnz = A.getNumEntriesInLocalRow(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;
482 static Scalar
Frobenius(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
490 const Map& AColMap = *A.getColMap();
491 const Map& BColMap = *B.getColMap();
495 size_t nnzA = 0, nnzB = 0;
511 size_t numRows = A.getNodeNumRows();
512 for (
size_t i = 0; i < numRows; i++) {
513 A.getLocalRowView(i, indA, valA);
514 B.getLocalRowView(i, indB, valB);
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 Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
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)
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.
Namespace for MueLu classes and methods.
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Exception throws to report incompatible objects (like maps).
static void PauseForDebugger()
#define MueLu_sumAll(rcpComm, in, out)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual int getRank() 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
static void seedrandom(unsigned int s)
Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Vector
static magnitudeType magnitude(T a)
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 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)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows.
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::RCP< Vector > GetLumpedMatrixDiagonal(Teuchos::RCP< const Matrix > rcpA)
Extract Matrix Diagonal of lumped matrix.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
virtual int getSize() const=0