MueLu  Version of the Day
MueLu_UtilitiesBase_decl.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_UTILITIESBASE_DECL_HPP
47 #define MUELU_UTILITIESBASE_DECL_HPP
48 
49 #include <unistd.h> //necessary for "sleep" function in debugging methods
50 #include <string>
51 
52 #include "MueLu_ConfigDefs.hpp"
53 
54 #include <Teuchos_DefaultComm.hpp>
55 #include <Teuchos_ScalarTraits.hpp>
56 #include <Teuchos_ParameterList.hpp>
57 
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>
71 
72 #include <Xpetra_Import.hpp>
73 #include <Xpetra_ImportFactory.hpp>
74 #include <Xpetra_MatrixMatrix.hpp>
75 #include <Xpetra_CrsMatrixWrap.hpp>
76 
77 
78 #include "MueLu_Exceptions.hpp"
79 
80 
81 
82 
83 namespace MueLu {
84 
85 // MPI helpers
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))
92 
100  template <class Scalar,
101  class LocalOrdinal = int,
102  class GlobalOrdinal = LocalOrdinal,
103  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
104  class UtilitiesBase {
105  public:
106 #undef MUELU_UTILITIESBASE_SHORT
107 //#include "MueLu_UseShortNames.hpp"
108  private:
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;
115  public:
116  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
117 
118 
119  static RCP<Matrix> Crs2Op(RCP<CrsMatrix> Op) {
120  if (Op.is_null())
121  return Teuchos::null;
122  return rcp(new CrsMatrixWrap(Op));
123  }
124 
131  static Teuchos::ArrayRCP<Scalar> GetMatrixDiagonal(const Matrix& A) {
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);
138  LocalOrdinal j = 0;
139  for (; j < cols.size(); ++j) {
140  if (Teuchos::as<size_t>(cols[j]) == i) {
141  diag[i] = vals[j];
142  break;
143  }
144  }
145  if (j == cols.size()) {
146  // Diagonal entry is absent
147  diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
148  }
149  }
150  return diag;
151  }
152 
159  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100) {
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);
168  LocalOrdinal j = 0;
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];
173  else
174  diagVals[i]=Teuchos::ScalarTraits<Scalar>::zero();
175  break;
176  }
177  }
178  if (j == cols.size()) {
179  // Diagonal entry is absent
180  diagVals[i]=Teuchos::ScalarTraits<Scalar>::zero();
181  }
182  }
183  diagVals=null;
184  return diag;
185  }
186 
187 
188 
195  static Teuchos::ArrayRCP<Scalar> GetLumpedMatrixDiagonal(const Matrix& A) {
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]);
205  }
206  }
207  return diag;
208  }
209 
217  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A) {
218  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
219  RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
220 
221  try {
222  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
223  if (crsOp == NULL) {
224  throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
225  }
226  Teuchos::ArrayRCP<size_t> offsets;
227  crsOp->getLocalDiagOffsets(offsets);
228  crsOp->getLocalDiagCopy(*localDiag,offsets());
229  }
230  catch (...) {
231  ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
232  Teuchos::ArrayRCP<Scalar> diagVals = GetMatrixDiagonal(A);
233  for (LocalOrdinal i = 0; i < localDiagVals.size(); i++)
234  localDiagVals[i] = diagVals[i];
235  localDiagVals = diagVals = null;
236  }
237 
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);
243  }
244  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
245  return diagonal;
246  }
247 
248 
249  // TODO: should NOT return an Array. Definition must be changed to:
250  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
251  // or
252  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
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);
258  RES->norm2(norms);
259  return norms;
260  }
261 
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); // no need to initialize to zero
267  Op.apply(X, *RES, Teuchos::NO_TRANS, one, zero);
268  RES->update(one, RHS, negone);
269  return RES;
270  }
271 
272 #ifndef _WIN32
273 #include <unistd.h>
274  static void PauseForDebugger() {
275  RCP<const Teuchos::Comm<int> > comm = Teuchos::DefaultComm<int>::getComm();
276  int myPID = comm->getRank();
277  int pid = getpid();
278  char hostname[80];
279  for (int i = 0; i <comm->getSize(); i++) {
280  if (i == myPID) {
281  gethostname(hostname, sizeof(hostname));
282  std::cout << "Host: " << hostname << "\tMPI rank: " << myPID << ",\tPID: " << pid << "\n\tattach " << pid << std::endl;
283  sleep(1);
284  }
285  }
286  if (myPID == 0) {
287  std::cout << "** Enter a character to continue > " << std::endl;
288  char go = ' ';
289  int r = scanf("%c", &go);
290  (void)r;
291  assert(r > 0);
292  }
293  comm->barrier();
294  }
295 #else
296  static void PauseForDebugger() {
297  throw(Exceptions::RuntimeError("MueLu Utils: PauseForDebugger not implemented on Windows."));
298  }
299 #endif
300 
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) {
318  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
319  "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
320 
321  // Create three vectors, fill z with random numbers
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());
325 
326  z->setSeed(seed); // seed random number generator
327  z->randomize(true);// use Xpetra implementation: -> same results for Epetra and Tpetra
328 
329  Teuchos::Array<Magnitude> norms(1);
330 
331  typedef Teuchos::ScalarTraits<Scalar> STS;
332 
333  const Scalar zero = STS::zero(), one = STS::one();
334 
335  Scalar lambda = zero;
336  Magnitude residual = STS::magnitude(zero);
337 
338  // power iteration
339  RCP<Vector> diagInvVec;
340  if (scaleByDiag) {
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);
345  }
346 
347  for (int iter = 0; iter < niters; ++iter) {
348  z->norm2(norms); // Compute 2-norm of z
349  q->update(one/norms[0], *z, zero); // Set q = z / normz
350  A.apply(*q, *z); // Compute z = A*q
351  if (scaleByDiag)
352  z->elementWiseMultiply(one, *diagInvVec, *z, zero);
353  lambda = q->dot(*z); // Approximate maximum eigenvalue: lamba = dot(q,z)
354 
355  if (iter % 100 == 0 || iter + 1 == niters) {
356  r->update(1.0, *z, -lambda, *q, zero); // Compute A*q - lambda*q
357  r->norm2(norms);
358  residual = STS::magnitude(norms[0] / lambda);
359  if (verbose) {
360  std::cout << "Iter = " << iter
361  << " Lambda = " << lambda
362  << " Residual of A*q - lambda*q = " << residual
363  << std::endl;
364  }
365  }
366  if (residual < tolerance)
367  break;
368  }
369  return lambda;
370  }
371 
372 
373 
374  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) {
375  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
376  return fancy;
377  }
378 
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();
385 
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]);
390  }
391  return Teuchos::ScalarTraits<Scalar>::magnitude(d);
392  }
393 
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);
410  if (nnz > 1)
411  for (size_t col = 0; col < nnz; col++)
412  if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
413  boundaryNodes[row] = false;
414  break;
415  }
416  }
417  return boundaryNodes;
418  }
419 
424  static Scalar Frobenius(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
425  // We check only row maps. Column may be different. One would hope that they are the same, as we typically
426  // calculate frobenius norm of the specified sparsity pattern with an updated matrix from the previous step,
427  // but matrix addition, even when one is submatrix of the other, changes column map (though change may be as
428  // simple as couple of elements swapped)
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");
431 
432  const Map& AColMap = *A.getColMap();
433  const Map& BColMap = *B.getColMap();
434 
435  Teuchos::ArrayView<const LocalOrdinal> indA, indB;
436  Teuchos::ArrayView<const Scalar> valA, valB;
437  size_t nnzA = 0, nnzB = 0;
438 
439  // We use a simple algorithm
440  // for each row we fill valBAll array with the values in the corresponding row of B
441  // as such, it serves as both sorted array and as storage, so we don't need to do a
442  // tricky problem: "find a value in the row of B corresponding to the specific GID"
443  // Once we do that, we translate LID of entries of row of A to LID of B, and multiply
444  // corresponding entries.
445  // The algorithm should be reasonably cheap, as it does not sort anything, provided
446  // that getLocalElement and getGlobalElement functions are reasonably effective. It
447  // *is* possible that the costs are hidden in those functions, but if maps are close
448  // to linear maps, we should be fine
449  Teuchos::Array<Scalar> valBAll(BColMap.getNodeNumElements());
450 
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);
457  nnzA = indA.size();
458  nnzB = indB.size();
459 
460  // Set up array values
461  for (size_t j = 0; j < nnzB; j++)
462  valBAll[indB[j]] = valB[j];
463 
464  for (size_t j = 0; j < nnzA; j++) {
465  // The cost of the whole Frobenius dot product function depends on the
466  // cost of the getLocalElement and getGlobalElement functions here.
467  LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
468  if (ind != invalid)
469  f += valBAll[ind] * valA[j];
470  }
471 
472  // Clean up array values
473  for (size_t j = 0; j < nnzB; j++)
474  valBAll[indB[j]] = zero;
475  }
476 
477  MueLu_sumAll(AColMap.getComm(), f, gf);
478 
479  return gf;
480  }
481 
491  static void SetRandomSeed(const Teuchos::Comm<int> &comm) {
492  // Distribute the seeds evenly in [1,maxint-1]. This guarantees nothing
493  // about where in random number stream we are, but avoids overflow situations
494  // in parallel when multiplying by a PID. It would be better to use
495  // a good parallel random number generator.
496  double one = 1.0;
497  int maxint = INT_MAX; //= 2^31-1 = 2147483647 for 32-bit integers
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].";
502  throw Exceptions::RuntimeError(errStr.str());
503  }
504  std::srand(mySeed);
505  // For Tpetra, we could use Kokkos' random number generator here.
506  Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
507  // Epetra
508  // MultiVector::Random() -> Epetra_Util::RandomDouble() -> Epetra_Utils::RandomInt()
509  // Its own random number generator, based on Seed_. Seed_ is initialized in Epetra_Util constructor with std::rand()
510  // So our setting std::srand() affects that too
511  }
512 
513 
514  }; // class Utils
515 
516 
517 
519 
520 } //namespace MueLu
521 
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).
#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.