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> 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> 86 #define MueLu_sumAll(rcpComm, in, out) \ 87 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out)) 88 #define MueLu_minAll(rcpComm, in, out) \ 89 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out)) 90 #define MueLu_maxAll(rcpComm, in, out) \ 91 Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out)) 100 template <
class Scalar,
101 class LocalOrdinal = int,
102 class GlobalOrdinal = LocalOrdinal,
103 class Node = KokkosClassic::DefaultNode::DefaultNodeType>
104 class UtilitiesBase {
106 #undef MUELU_UTILITIESBASE_SHORT 109 typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>
CrsMatrixWrap;
110 typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
CrsMatrix;
111 typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>
Matrix;
112 typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
Vector;
113 typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>
MultiVector;
114 typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>
Map;
116 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Magnitude;
119 static RCP<Matrix>
Crs2Op(RCP<CrsMatrix> Op) {
121 return Teuchos::null;
132 size_t numRows = A.getRowMap()->getNodeNumElements();
133 Teuchos::ArrayRCP<Scalar> diag(numRows);
134 Teuchos::ArrayView<const LocalOrdinal> cols;
135 Teuchos::ArrayView<const Scalar> vals;
136 for (
size_t i = 0; i < numRows; ++i) {
137 A.getLocalRowView(i, cols, vals);
139 for (; j < cols.size(); ++j) {
140 if (Teuchos::as<size_t>(cols[j]) == i) {
145 if (j == cols.size()) {
147 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
160 RCP<const Map> rowMap = A.getRowMap();
161 RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
162 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
163 size_t numRows = rowMap->getNodeNumElements();
164 Teuchos::ArrayView<const LocalOrdinal> cols;
165 Teuchos::ArrayView<const Scalar> vals;
166 for (
size_t i = 0; i < numRows; ++i) {
167 A.getLocalRowView(i, cols, vals);
169 for (; j < cols.size(); ++j) {
170 if (Teuchos::as<size_t>(cols[j]) == i) {
171 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]) > tol)
172 diagVals[i] = Teuchos::ScalarTraits<Scalar>::one() / vals[j];
174 diagVals[i]=Teuchos::ScalarTraits<Scalar>::zero();
178 if (j == cols.size()) {
180 diagVals[i]=Teuchos::ScalarTraits<Scalar>::zero();
196 size_t numRows = A.getRowMap()->getNodeNumElements();
197 Teuchos::ArrayRCP<Scalar> diag(numRows);
198 Teuchos::ArrayView<const LocalOrdinal> cols;
199 Teuchos::ArrayView<const Scalar> vals;
200 for (
size_t i = 0; i < numRows; ++i) {
201 A.getLocalRowView(i, cols, vals);
202 diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
203 for (LocalOrdinal j = 0; j < cols.size(); ++j) {
204 diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
218 RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
219 RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
222 const CrsMatrixWrap* crsOp =
dynamic_cast<const CrsMatrixWrap*
>(&A);
226 Teuchos::ArrayRCP<size_t> offsets;
227 crsOp->getLocalDiagOffsets(offsets);
228 crsOp->getLocalDiagCopy(*localDiag,offsets());
231 ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
233 for (LocalOrdinal i = 0; i < localDiagVals.size(); i++)
234 localDiagVals[i] = diagVals[i];
235 localDiagVals = diagVals = null;
238 RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
239 RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
240 importer = A.getCrsGraph()->getImporter();
241 if (importer == Teuchos::null) {
242 importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
244 diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
253 static Teuchos::Array<Magnitude>
ResidualNorm(
const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op,
const MultiVector& X,
const MultiVector& RHS) {
254 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
255 const size_t numVecs = X.getNumVectors();
256 RCP<MultiVector> RES =
Residual(Op, X, RHS);
257 Teuchos::Array<Magnitude> norms(numVecs);
262 static RCP<MultiVector>
Residual(
const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op,
const MultiVector& X,
const MultiVector& RHS) {
263 TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(),
Exceptions::RuntimeError,
"Number of solution vectors != number of right-hand sides")
264 const size_t numVecs = X.getNumVectors();
265 Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
266 RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRangeMap(), numVecs,
false);
267 Op.apply(X, *RES, Teuchos::NO_TRANS, one, zero);
268 RES->update(one, RHS, negone);
275 RCP<const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
276 int myPID = comm->getRank();
279 for (
int i = 0; i <comm->getSize(); i++) {
281 gethostname(hostname,
sizeof(hostname));
282 std::cout <<
"Host: " << hostname <<
"\tMPI rank: " << myPID <<
",\tPID: " << pid <<
"\n\tattach " << pid << std::endl;
287 std::cout <<
"** Enter a character to continue > " << std::endl;
289 int r = scanf(
"%c", &go);
316 static Scalar
PowerMethod(
const Matrix& A,
bool scaleByDiag =
true,
317 LocalOrdinal niters = 10, Magnitude tolerance = 1e-2,
bool verbose =
false,
unsigned int seed = 123) {
319 "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
322 RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
323 RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
324 RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
329 Teuchos::Array<Magnitude> norms(1);
331 typedef Teuchos::ScalarTraits<Scalar> STS;
333 const Scalar zero = STS::zero(), one = STS::one();
335 Scalar lambda = zero;
336 Magnitude residual = STS::magnitude(zero);
339 RCP<Vector> diagInvVec;
341 RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
342 A.getLocalDiagCopy(*diagVec);
343 diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
344 diagInvVec->reciprocal(*diagVec);
347 for (
int iter = 0; iter < niters; ++iter) {
349 q->update(one/norms[0], *z, zero);
352 z->elementWiseMultiply(one, *diagInvVec, *z, zero);
355 if (iter % 100 == 0 || iter + 1 == niters) {
356 r->update(1.0, *z, -lambda, *q, zero);
358 residual = STS::magnitude(norms[0] / lambda);
360 std::cout <<
"Iter = " << iter
361 <<
" Lambda = " << lambda
362 <<
" Residual of A*q - lambda*q = " << residual
366 if (residual < tolerance)
374 static RCP<Teuchos::FancyOStream>
MakeFancy(std::ostream& os) {
375 RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
383 static typename Teuchos::ScalarTraits<Scalar>::magnitudeType
Distance2(
const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& v, LocalOrdinal i0, LocalOrdinal i1) {
384 size_t numVectors = v.getNumVectors();
386 Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
387 for (
size_t j = 0; j < numVectors; j++) {
388 Teuchos::ArrayRCP<const Scalar> vv = v.getData(j);
389 d += (vv[i0] - vv[i1])*(vv[i0] - vv[i1]);
391 return Teuchos::ScalarTraits<Scalar>::magnitude(d);
401 static Teuchos::ArrayRCP<const bool>
DetectDirichletRows(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
402 LocalOrdinal numRows = A.getNodeNumRows();
403 typedef Teuchos::ScalarTraits<Scalar> STS;
404 ArrayRCP<bool> boundaryNodes(numRows,
true);
405 for (LocalOrdinal row = 0; row < numRows; row++) {
406 ArrayView<const LocalOrdinal> indices;
407 ArrayView<const Scalar> vals;
408 A.getLocalRowView(row, indices, vals);
409 size_t nnz = A.getNumEntriesInLocalRow(row);
411 for (
size_t col = 0; col < nnz; col++)
412 if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
413 boundaryNodes[row] =
false;
417 return boundaryNodes;
424 static Scalar
Frobenius(
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
429 TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()),
Exceptions::Incompatible,
"MueLu::CGSolver::Frobenius: row maps are incompatible");
430 TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(),
Exceptions::RuntimeError,
"Matrices must be fill completed");
432 const Map& AColMap = *A.getColMap();
433 const Map& BColMap = *B.getColMap();
435 Teuchos::ArrayView<const LocalOrdinal> indA, indB;
436 Teuchos::ArrayView<const Scalar> valA, valB;
437 size_t nnzA = 0, nnzB = 0;
449 Teuchos::Array<Scalar> valBAll(BColMap.getNodeNumElements());
451 LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
452 Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
453 size_t numRows = A.getNodeNumRows();
454 for (
size_t i = 0; i < numRows; i++) {
455 A.getLocalRowView(i, indA, valA);
456 B.getLocalRowView(i, indB, valB);
461 for (
size_t j = 0; j < nnzB; j++)
462 valBAll[indB[j]] = valB[j];
464 for (
size_t j = 0; j < nnzA; j++) {
467 LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
469 f += valBAll[ind] * valA[j];
473 for (
size_t j = 0; j < nnzB; j++)
474 valBAll[indB[j]] = zero;
497 int maxint = INT_MAX;
498 int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.getRank()+1)/(comm.getSize()+one)) );
499 if (mySeed < 1 || mySeed == maxint) {
500 std::ostringstream errStr;
501 errStr <<
"Error detected with random seed = " << mySeed <<
". It should be in the interval [1,2^31-2].";
506 Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
522 #define MUELU_UTILITIESBASE_SHORT 523 #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
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.
Exception throws to report incompatible objects (like maps).
static void PauseForDebugger()
#define MueLu_sumAll(rcpComm, in, out)
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
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)
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< 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 RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100)
Extract Matrix Diagonal.
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.