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>
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 #include "MueLu_Exceptions.hpp"
78 
79 
80 
81 
82 namespace MueLu {
83 
84 // MPI helpers
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))
91 
99  template <class Scalar,
100  class LocalOrdinal = int,
101  class GlobalOrdinal = LocalOrdinal,
102  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
103  class UtilitiesBase {
104  public:
105 #undef MUELU_UTILITIESBASE_SHORT
106 //#include "MueLu_UseShortNames.hpp"
107  private:
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;
114  public:
116 
117 
119  if (Op.is_null())
120  return Teuchos::null;
121  return rcp(new CrsMatrixWrap(Op));
122  }
123 
131  size_t numRows = A.getRowMap()->getNodeNumElements();
132  Teuchos::ArrayRCP<Scalar> diag(numRows);
135  for (size_t i = 0; i < numRows; ++i) {
136  A.getLocalRowView(i, cols, vals);
137  LocalOrdinal j = 0;
138  for (; j < cols.size(); ++j) {
139  if (Teuchos::as<size_t>(cols[j]) == i) {
140  diag[i] = vals[j];
141  break;
142  }
143  }
144  if (j == cols.size()) {
145  // Diagonal entry is absent
147  }
148  }
149  return diag;
150  }
151 
159 
160  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
161 
162  RCP<const Map> rowMap = rcpA->getRowMap();
163  RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
164 
165  rcpA->getLocalDiagCopy(*diag);
166 
168 
169  return inv;
170  }
171 
172 
173 
181 
182  size_t numRows = A.getRowMap()->getNodeNumElements();
183  Teuchos::ArrayRCP<Scalar> diag(numRows);
186  for (size_t i = 0; i < numRows; ++i) {
187  A.getLocalRowView(i, cols, vals);
189  for (LocalOrdinal j = 0; j < cols.size(); ++j) {
190  diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
191  }
192  }
193  return diag;
194  }
195 
203 
204  RCP<Vector> diag = Teuchos::null;
205 
207  Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
208  if(bA == Teuchos::null) {
209  RCP<const Map> rowMap = rcpA->getRowMap();
210  diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
211  ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
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]);
219  }
220  }
221 
222  } else {
223  //TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != bA->Cols(), Xpetra::Exceptions::RuntimeError,
224  // "UtilitiesBase::GetLumpedMatrixDiagonal(): you cannot extract the diagonal of a "<< bA->Rows() << "x"<< bA->Cols() << " blocked matrix." );
225 
226  diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),true);
227 
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()) {
231  // if we are in Thyra mode, but the block (row,row) is again a blocked operator, we have to use (pseudo) Xpetra-style GIDs with offset!
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);
234  RCP<const Vector> dd = GetLumpedMatrixDiagonal(bA->getMatrix(row,col));
235  ddtemp->update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
236  bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
237  }
238  }
239  }
240 
241  }
242 
243  // we should never get here...
244  return diag;
245  }
246 
255 
256  RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),true);
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];
262  else
263  retVals[i] = tolReplacement;
264  }
265  return ret;
266  }
267 
276  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
277  RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
278 
279  try {
280  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
281  if (crsOp == NULL) {
282  throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
283  }
285  crsOp->getLocalDiagOffsets(offsets);
286  crsOp->getLocalDiagCopy(*localDiag,offsets());
287  }
288  catch (...) {
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;
294  }
295 
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);
301  }
302  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
303  return diagonal;
304  }
305 
306 
307  // TODO: should NOT return an Array. Definition must be changed to:
308  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
309  // or
310  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
311  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS) {
312  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
313  const size_t numVecs = X.getNumVectors();
314  RCP<MultiVector> RES = Residual(Op, X, RHS);
315  Teuchos::Array<Magnitude> norms(numVecs);
316  RES->norm2(norms);
317  return norms;
318  }
319 
320  static RCP<MultiVector> Residual(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS) {
321  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
322  const size_t numVecs = X.getNumVectors();
323  Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
324  RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(Op.getRangeMap(), numVecs, false); // no need to initialize to zero
325  Op.apply(X, *RES, Teuchos::NO_TRANS, one, zero);
326  RES->update(one, RHS, negone);
327  return RES;
328  }
329 
330 #ifndef _WIN32
331 #include <unistd.h>
332  static void PauseForDebugger() {
334  int myPID = comm->getRank();
335  int pid = getpid();
336  char hostname[80];
337  for (int i = 0; i <comm->getSize(); i++) {
338  if (i == myPID) {
339  gethostname(hostname, sizeof(hostname));
340  std::cout << "Host: " << hostname << "\tMPI rank: " << myPID << ",\tPID: " << pid << "\n\tattach " << pid << std::endl;
341  sleep(1);
342  }
343  }
344  if (myPID == 0) {
345  std::cout << "** Enter a character to continue > " << std::endl;
346  char go = ' ';
347  int r = scanf("%c", &go);
348  (void)r;
349  assert(r > 0);
350  }
351  comm->barrier();
352  }
353 #else
354  static void PauseForDebugger() {
355  throw(Exceptions::RuntimeError("MueLu Utils: PauseForDebugger not implemented on Windows."));
356  }
357 #endif
358 
374  static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true,
375  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
376  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
377  "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
378 
379  // Create three vectors, fill z with random numbers
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());
383 
384  z->setSeed(seed); // seed random number generator
385  z->randomize(true);// use Xpetra implementation: -> same results for Epetra and Tpetra
386 
387  Teuchos::Array<Magnitude> norms(1);
388 
389  typedef Teuchos::ScalarTraits<Scalar> STS;
390 
391  const Scalar zero = STS::zero(), one = STS::one();
392 
393  Scalar lambda = zero;
394  Magnitude residual = STS::magnitude(zero);
395 
396  // power iteration
397  RCP<Vector> diagInvVec;
398  if (scaleByDiag) {
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);
403  }
404 
405  for (int iter = 0; iter < niters; ++iter) {
406  z->norm2(norms); // Compute 2-norm of z
407  q->update(one/norms[0], *z, zero); // Set q = z / normz
408  A.apply(*q, *z); // Compute z = A*q
409  if (scaleByDiag)
410  z->elementWiseMultiply(one, *diagInvVec, *z, zero);
411  lambda = q->dot(*z); // Approximate maximum eigenvalue: lamba = dot(q,z)
412 
413  if (iter % 100 == 0 || iter + 1 == niters) {
414  r->update(1.0, *z, -lambda, *q, zero); // Compute A*q - lambda*q
415  r->norm2(norms);
416  residual = STS::magnitude(norms[0] / lambda);
417  if (verbose) {
418  std::cout << "Iter = " << iter
419  << " Lambda = " << lambda
420  << " Residual of A*q - lambda*q = " << residual
421  << std::endl;
422  }
423  }
424  if (residual < tolerance)
425  break;
426  }
427  return lambda;
428  }
429 
430 
431 
432  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) {
433  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
434  return fancy;
435  }
436 
441  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& v, LocalOrdinal i0, LocalOrdinal i1) {
442  size_t numVectors = v.getNumVectors();
443 
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]);
448  }
450  }
451 
459  static Teuchos::ArrayRCP<const bool> DetectDirichletRows(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
460  LocalOrdinal numRows = A.getNodeNumRows();
461  typedef Teuchos::ScalarTraits<Scalar> STS;
462  ArrayRCP<bool> boundaryNodes(numRows, true);
463  for (LocalOrdinal row = 0; row < numRows; row++) {
466  A.getLocalRowView(row, indices, vals);
467  size_t nnz = A.getNumEntriesInLocalRow(row);
468  if (nnz > 1)
469  for (size_t col = 0; col < nnz; col++)
470  if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
471  boundaryNodes[row] = false;
472  break;
473  }
474  }
475  return boundaryNodes;
476  }
477 
482  static Scalar Frobenius(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
483  // We check only row maps. Column may be different. One would hope that they are the same, as we typically
484  // calculate frobenius norm of the specified sparsity pattern with an updated matrix from the previous step,
485  // but matrix addition, even when one is submatrix of the other, changes column map (though change may be as
486  // simple as couple of elements swapped)
487  TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()), Exceptions::Incompatible, "MueLu::CGSolver::Frobenius: row maps are incompatible");
488  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(), Exceptions::RuntimeError, "Matrices must be fill completed");
489 
490  const Map& AColMap = *A.getColMap();
491  const Map& BColMap = *B.getColMap();
492 
495  size_t nnzA = 0, nnzB = 0;
496 
497  // We use a simple algorithm
498  // for each row we fill valBAll array with the values in the corresponding row of B
499  // as such, it serves as both sorted array and as storage, so we don't need to do a
500  // tricky problem: "find a value in the row of B corresponding to the specific GID"
501  // Once we do that, we translate LID of entries of row of A to LID of B, and multiply
502  // corresponding entries.
503  // The algorithm should be reasonably cheap, as it does not sort anything, provided
504  // that getLocalElement and getGlobalElement functions are reasonably effective. It
505  // *is* possible that the costs are hidden in those functions, but if maps are close
506  // to linear maps, we should be fine
507  Teuchos::Array<Scalar> valBAll(BColMap.getNodeNumElements());
508 
509  LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
510  Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
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);
515  nnzA = indA.size();
516  nnzB = indB.size();
517 
518  // Set up array values
519  for (size_t j = 0; j < nnzB; j++)
520  valBAll[indB[j]] = valB[j];
521 
522  for (size_t j = 0; j < nnzA; j++) {
523  // The cost of the whole Frobenius dot product function depends on the
524  // cost of the getLocalElement and getGlobalElement functions here.
525  LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
526  if (ind != invalid)
527  f += valBAll[ind] * valA[j];
528  }
529 
530  // Clean up array values
531  for (size_t j = 0; j < nnzB; j++)
532  valBAll[indB[j]] = zero;
533  }
534 
535  MueLu_sumAll(AColMap.getComm(), f, gf);
536 
537  return gf;
538  }
539 
549  static void SetRandomSeed(const Teuchos::Comm<int> &comm) {
550  // Distribute the seeds evenly in [1,maxint-1]. This guarantees nothing
551  // about where in random number stream we are, but avoids overflow situations
552  // in parallel when multiplying by a PID. It would be better to use
553  // a good parallel random number generator.
554  double one = 1.0;
555  int maxint = INT_MAX; //= 2^31-1 = 2147483647 for 32-bit integers
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].";
560  throw Exceptions::RuntimeError(errStr.str());
561  }
562  std::srand(mySeed);
563  // For Tpetra, we could use Kokkos' random number generator here.
565  // Epetra
566  // MultiVector::Random() -> Epetra_Util::RandomDouble() -> Epetra_Utils::RandomInt()
567  // Its own random number generator, based on Seed_. Seed_ is initialized in Epetra_Util constructor with std::rand()
568  // So our setting std::srand() affects that too
569  }
570 
571 
572  }; // class Utils
573 
574 
575 
577 
578 } //namespace MueLu
579 
580 #define MUELU_UTILITIESBASE_SHORT
581 #endif // MUELU_UTILITIESBASE_DECL_HPP
Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
size_type size() 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)
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.
size_type size() const
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).
#define MueLu_sumAll(rcpComm, in, out)
bool is_null() const
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