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 #ifndef _WIN32
50 #include <unistd.h> //necessary for "sleep" function in debugging methods (PauseForDebugging)
51 #endif
52 
53 #include <string>
54 
55 #include "MueLu_ConfigDefs.hpp"
56 
57 #include <Teuchos_DefaultComm.hpp>
58 #include <Teuchos_ScalarTraits.hpp>
60 
61 #include <Xpetra_BlockedCrsMatrix_fwd.hpp>
62 #include <Xpetra_CrsMatrix_fwd.hpp>
63 #include <Xpetra_CrsMatrixWrap_fwd.hpp>
64 #include <Xpetra_Map_fwd.hpp>
65 #include <Xpetra_BlockedMap_fwd.hpp>
66 #include <Xpetra_MapFactory_fwd.hpp>
67 #include <Xpetra_Matrix_fwd.hpp>
68 #include <Xpetra_MatrixFactory_fwd.hpp>
69 #include <Xpetra_MultiVector_fwd.hpp>
70 #include <Xpetra_MultiVectorFactory_fwd.hpp>
71 #include <Xpetra_Operator_fwd.hpp>
72 #include <Xpetra_Vector_fwd.hpp>
73 #include <Xpetra_BlockedMultiVector.hpp>
74 #include <Xpetra_BlockedVector.hpp>
75 #include <Xpetra_VectorFactory_fwd.hpp>
76 #include <Xpetra_ExportFactory.hpp>
77 
78 #include <Xpetra_Import.hpp>
79 #include <Xpetra_ImportFactory.hpp>
80 #include <Xpetra_MatrixMatrix.hpp>
81 #include <Xpetra_CrsMatrixWrap.hpp>
82 #include <Xpetra_StridedMap.hpp>
83 
84 #include "MueLu_Exceptions.hpp"
85 
86 
87 
88 
89 namespace MueLu {
90 
91 // MPI helpers
92 #define MueLu_sumAll(rcpComm, in, out) \
93  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out))
94 #define MueLu_minAll(rcpComm, in, out) \
95  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out))
96 #define MueLu_maxAll(rcpComm, in, out) \
97  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out))
98 
106  template <class Scalar,
107  class LocalOrdinal = int,
108  class GlobalOrdinal = LocalOrdinal,
109  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
110  class UtilitiesBase {
111  public:
112 #undef MUELU_UTILITIESBASE_SHORT
113 //#include "MueLu_UseShortNames.hpp"
114  private:
115  typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> CrsMatrixWrap;
116  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> CrsMatrix;
117  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Matrix;
118  typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Vector;
119  typedef Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> BlockedVector;
120  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MultiVector;
121  typedef Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> BlockedMultiVector;
122  typedef Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> BlockedMap;
123  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
124  public:
126 
127 
129  if (Op.is_null())
130  return Teuchos::null;
131  return rcp(new CrsMatrixWrap(Op));
132  }
133 
141  size_t numRows = A.getRowMap()->getNodeNumElements();
142  Teuchos::ArrayRCP<Scalar> diag(numRows);
145  for (size_t i = 0; i < numRows; ++i) {
146  A.getLocalRowView(i, cols, vals);
147  LocalOrdinal j = 0;
148  for (; j < cols.size(); ++j) {
149  if (Teuchos::as<size_t>(cols[j]) == i) {
150  diag[i] = vals[j];
151  break;
152  }
153  }
154  if (j == cols.size()) {
155  // Diagonal entry is absent
157  }
158  }
159  return diag;
160  }
161 
169 
170  RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
171 
172  RCP<const Map> rowMap = rcpA->getRowMap();
173  RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
174 
175  rcpA->getLocalDiagCopy(*diag);
176 
178 
179  return inv;
180  }
181 
182 
183 
191 
192  size_t numRows = A.getRowMap()->getNodeNumElements();
193  Teuchos::ArrayRCP<Scalar> diag(numRows);
196  for (size_t i = 0; i < numRows; ++i) {
197  A.getLocalRowView(i, cols, vals);
199  for (LocalOrdinal j = 0; j < cols.size(); ++j) {
200  diag[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
201  }
202  }
203  return diag;
204  }
205 
213 
214  RCP<Vector> diag = Teuchos::null;
215 
217  Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
218  if(bA == Teuchos::null) {
219  RCP<const Map> rowMap = rcpA->getRowMap();
220  diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
221  ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
224  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
225  rcpA->getLocalRowView(i, cols, vals);
226  diagVals[i] = Teuchos::ScalarTraits<Scalar>::zero();
227  for (LocalOrdinal j = 0; j < cols.size(); ++j) {
228  diagVals[i] += Teuchos::ScalarTraits<Scalar>::magnitude(vals[j]);
229  }
230  }
231 
232  } else {
233  //TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != bA->Cols(), Xpetra::Exceptions::RuntimeError,
234  // "UtilitiesBase::GetLumpedMatrixDiagonal(): you cannot extract the diagonal of a "<< bA->Rows() << "x"<< bA->Cols() << " blocked matrix." );
235 
236  diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),true);
237 
238  for (size_t row = 0; row < bA->Rows(); ++row) {
239  for (size_t col = 0; col < bA->Cols(); ++col) {
240  if (!bA->getMatrix(row,col).is_null()) {
241  // 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!
242  bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA->getMatrix(row,col)) == Teuchos::null);
243  RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag,row,bThyraMode);
244  RCP<const Vector> dd = GetLumpedMatrixDiagonal(bA->getMatrix(row,col));
245  ddtemp->update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
246  bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
247  }
248  }
249  }
250 
251  }
252 
253  // we should never get here...
254  return diag;
255  }
256 
265 
266  RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),true);
267 
268  // check whether input vector "v" is a BlockedVector
269  RCP<const BlockedVector> bv = Teuchos::rcp_dynamic_cast<const BlockedVector>(v);
270  if(bv.is_null() == false) {
271  RCP<BlockedVector> bret = Teuchos::rcp_dynamic_cast<BlockedVector>(ret);
272  TEUCHOS_TEST_FOR_EXCEPTION(bret.is_null() == true, MueLu::Exceptions::RuntimeError,"MueLu::UtilitiesBase::GetInverse: return vector should be of type BlockedVector");
273  RCP<const BlockedMap> bmap = bv->getBlockedMap();
274  for(size_t r = 0; r < bmap->getNumMaps(); ++r) {
275  RCP<const MultiVector> submvec = bv->getMultiVector(r,bmap->getThyraMode());
276  RCP<const Vector> subvec = submvec->getVector(0);
278  bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
279  }
280  return ret;
281  }
282 
283  // v is an {Epetra,Tpetra}Vector: work with the underlying raw data
284  ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
285  ArrayRCP<const Scalar> inputVals = v->getData(0);
286  for (size_t i = 0; i < v->getMap()->getNodeNumElements(); ++i) {
287  if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
288  retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
289  else
290  retVals[i] = tolReplacement;
291  }
292  return ret;
293  }
294 
303  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
304 
305  // Undo block map (if we have one)
306  RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
307  if(!browMap.is_null()) rowMap = browMap->getMap();
308 
309  RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
310  try {
311  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
312  if (crsOp == NULL) {
313  throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
314  }
316  crsOp->getLocalDiagOffsets(offsets);
317  crsOp->getLocalDiagCopy(*localDiag,offsets());
318  }
319  catch (...) {
320  ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
322  for (LocalOrdinal i = 0; i < localDiagVals.size(); i++)
323  localDiagVals[i] = diagVals[i];
324  localDiagVals = diagVals = null;
325  }
326 
327  RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
329  importer = A.getCrsGraph()->getImporter();
330  if (importer == Teuchos::null) {
331  importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
332  }
333  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
334  return diagonal;
335  }
336 
337 
338  // TODO: should NOT return an Array. Definition must be changed to:
339  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
340  // or
341  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
342  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS) {
343  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
344  const size_t numVecs = X.getNumVectors();
345  RCP<MultiVector> RES = Residual(Op, X, RHS);
346  Teuchos::Array<Magnitude> norms(numVecs);
347  RES->norm2(norms);
348  return norms;
349  }
350 
351  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector & Resid) {
352  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
353  const size_t numVecs = X.getNumVectors();
354  Residual(Op,X,RHS,Resid);
355  Teuchos::Array<Magnitude> norms(numVecs);
356  Resid.norm2(norms);
357  return norms;
358  }
359 
360  static RCP<MultiVector> Residual(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS) {
361  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
362  const size_t numVecs = X.getNumVectors();
363  Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
364  // TODO Op.getRangeMap should return a BlockedMap if it is a BlockedCrsOperator
365  RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(RHS.getMap(), numVecs, false); // no need to initialize to zero
366  Op.apply(X, *RES, Teuchos::NO_TRANS, one, zero);
367  RES->update(one, RHS, negone);
368  return RES;
369  }
370 
371 
372  static void Residual(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector & Resid) {
373  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides");
374  TEUCHOS_TEST_FOR_EXCEPTION(Resid.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of residual vectors != number of right-hand sides");
375  Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
376  // TODO Op.getRangeMap should return a BlockedMap if it is a BlockedCrsOperator
377  Op.apply(X, Resid, Teuchos::NO_TRANS, one, zero);
378  Resid.update(one, RHS, negone);
379  }
380 
381 
382 #ifndef _WIN32
383  static void PauseForDebugger() {
385  int myPID = comm->getRank();
386  int pid = getpid();
387  char hostname[80];
388  for (int i = 0; i <comm->getSize(); i++) {
389  if (i == myPID) {
390  gethostname(hostname, sizeof(hostname));
391  std::cout << "Host: " << hostname << "\tMPI rank: " << myPID << ",\tPID: " << pid << "\n\tattach " << pid << std::endl;
392  sleep(1);
393  }
394  }
395  if (myPID == 0) {
396  std::cout << "** Enter a character to continue > " << std::endl;
397  char go = ' ';
398  int r = scanf("%c", &go);
399  (void)r;
400  assert(r > 0);
401  }
402  comm->barrier();
403  }
404 #else
405  static void PauseForDebugger() {
406  throw(Exceptions::RuntimeError("MueLu Utils: PauseForDebugger not implemented on Windows."));
407  }
408 #endif
409 
425  static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true,
426  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
427  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
428  "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
429 
430  // Create three vectors, fill z with random numbers
431  RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
432  RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
433  RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
434 
435  z->setSeed(seed); // seed random number generator
436  z->randomize(true);// use Xpetra implementation: -> same results for Epetra and Tpetra
437 
438  Teuchos::Array<Magnitude> norms(1);
439 
440  typedef Teuchos::ScalarTraits<Scalar> STS;
441 
442  const Scalar zero = STS::zero(), one = STS::one();
443 
444  Scalar lambda = zero;
445  Magnitude residual = STS::magnitude(zero);
446 
447  // power iteration
448  RCP<Vector> diagInvVec;
449  if (scaleByDiag) {
450  RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
451  A.getLocalDiagCopy(*diagVec);
452  diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
453  diagInvVec->reciprocal(*diagVec);
454  }
455 
456  for (int iter = 0; iter < niters; ++iter) {
457  z->norm2(norms); // Compute 2-norm of z
458  q->update(one/norms[0], *z, zero); // Set q = z / normz
459  A.apply(*q, *z); // Compute z = A*q
460  if (scaleByDiag)
461  z->elementWiseMultiply(one, *diagInvVec, *z, zero);
462  lambda = q->dot(*z); // Approximate maximum eigenvalue: lamba = dot(q,z)
463 
464  if (iter % 100 == 0 || iter + 1 == niters) {
465  r->update(1.0, *z, -lambda, *q, zero); // Compute A*q - lambda*q
466  r->norm2(norms);
467  residual = STS::magnitude(norms[0] / lambda);
468  if (verbose) {
469  std::cout << "Iter = " << iter
470  << " Lambda = " << lambda
471  << " Residual of A*q - lambda*q = " << residual
472  << std::endl;
473  }
474  }
475  if (residual < tolerance)
476  break;
477  }
478  return lambda;
479  }
480 
481 
482 
483  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) {
484  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
485  return fancy;
486  }
487 
492  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>> &v, LocalOrdinal i0, LocalOrdinal i1) {
493  const size_t numVectors = v.size();
494 
496  for (size_t j = 0; j < numVectors; j++) {
497  d += (v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
498  }
500  }
501 
514  static Teuchos::ArrayRCP<const bool> DetectDirichletRows(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero(), bool count_twos_as_dirichlet=false) {
515  LocalOrdinal numRows = A.getNodeNumRows();
516  typedef Teuchos::ScalarTraits<Scalar> STS;
517  ArrayRCP<bool> boundaryNodes(numRows, true);
518  if (count_twos_as_dirichlet) {
519  for (LocalOrdinal row = 0; row < numRows; row++) {
522  A.getLocalRowView(row, indices, vals);
523  size_t nnz = A.getNumEntriesInLocalRow(row);
524  if (nnz > 2) {
525  size_t col;
526  for (col = 0; col < nnz; col++)
527  if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
528  if (!boundaryNodes[row])
529  break;
530  boundaryNodes[row] = false;
531  }
532  if (col == nnz)
533  boundaryNodes[row] = true;
534  }
535  }
536  } else {
537  for (LocalOrdinal row = 0; row < numRows; row++) {
540  A.getLocalRowView(row, indices, vals);
541  size_t nnz = A.getNumEntriesInLocalRow(row);
542  if (nnz > 1)
543  for (size_t col = 0; col < nnz; col++)
544  if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
545  boundaryNodes[row] = false;
546  break;
547  }
548  }
549  }
550  return boundaryNodes;
551  }
552 
553 
566  static Teuchos::ArrayRCP<const bool> DetectDirichletRowsExt(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, bool & bHasZeroDiagonal, const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
567 
568  // assume that there is no zero diagonal in matrix
569  bHasZeroDiagonal = false;
570 
571  Teuchos::RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
572  A.getLocalDiagCopy(*diagVec);
573  Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
574 
575  LocalOrdinal numRows = A.getNodeNumRows();
576  typedef Teuchos::ScalarTraits<Scalar> STS;
577  ArrayRCP<bool> boundaryNodes(numRows, false);
578  for (LocalOrdinal row = 0; row < numRows; row++) {
581  A.getLocalRowView(row, indices, vals);
582  size_t nnz = 0; // collect nonzeros in row (excluding the diagonal)
583  bool bHasDiag = false;
584  for (decltype(indices.size()) col = 0; col < indices.size(); col++) {
585  if ( indices[col] != row) {
586  if (STS::magnitude(vals[col] / sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col])) ) > tol) {
587  nnz++;
588  }
589  } else bHasDiag = true; // found a diagonal entry
590  }
591  if (bHasDiag == false) bHasZeroDiagonal = true; // we found at least one row without a diagonal
592  else if(nnz == 0) boundaryNodes[row] = true;
593  }
594  return boundaryNodes;
595  }
596 
597 
608  static Teuchos::ArrayRCP<const bool> DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
609  const Teuchos::ArrayRCP<const bool>& dirichletRows) {
610  Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
611  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
614  Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,1);
615  myColsToZero->putScalar(zero);
616  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
617  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
618  if (dirichletRows[i]) {
621  A.getLocalRowView(i,indices,values);
622  for(size_t j=0; j<static_cast<size_t>(indices.size()); j++)
623  myColsToZero->replaceLocalValue(indices[j],0,one);
624  }
625  }
626 
627  Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,1);
628  globalColsToZero->putScalar(zero);
629  Teuchos::RCP<Xpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
630  // export to domain map
631  globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
632  // import to column map
633  myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
634  Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
635  Teuchos::ArrayRCP<bool> dirichletCols(colMap->getNodeNumElements(), true);
637  for(size_t i=0; i<colMap->getNodeNumElements(); i++) {
638  dirichletCols[i] = Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i])>2.0*eps;
639  }
640  return dirichletCols;
641  }
642 
643 
648  static Scalar Frobenius(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
649  // We check only row maps. Column may be different. One would hope that they are the same, as we typically
650  // calculate frobenius norm of the specified sparsity pattern with an updated matrix from the previous step,
651  // but matrix addition, even when one is submatrix of the other, changes column map (though change may be as
652  // simple as couple of elements swapped)
653  TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()), Exceptions::Incompatible, "MueLu::CGSolver::Frobenius: row maps are incompatible");
654  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(), Exceptions::RuntimeError, "Matrices must be fill completed");
655 
656  const Map& AColMap = *A.getColMap();
657  const Map& BColMap = *B.getColMap();
658 
661  size_t nnzA = 0, nnzB = 0;
662 
663  // We use a simple algorithm
664  // for each row we fill valBAll array with the values in the corresponding row of B
665  // as such, it serves as both sorted array and as storage, so we don't need to do a
666  // tricky problem: "find a value in the row of B corresponding to the specific GID"
667  // Once we do that, we translate LID of entries of row of A to LID of B, and multiply
668  // corresponding entries.
669  // The algorithm should be reasonably cheap, as it does not sort anything, provided
670  // that getLocalElement and getGlobalElement functions are reasonably effective. It
671  // *is* possible that the costs are hidden in those functions, but if maps are close
672  // to linear maps, we should be fine
673  Teuchos::Array<Scalar> valBAll(BColMap.getNodeNumElements());
674 
675  LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
676  Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
677  size_t numRows = A.getNodeNumRows();
678  for (size_t i = 0; i < numRows; i++) {
679  A.getLocalRowView(i, indA, valA);
680  B.getLocalRowView(i, indB, valB);
681  nnzA = indA.size();
682  nnzB = indB.size();
683 
684  // Set up array values
685  for (size_t j = 0; j < nnzB; j++)
686  valBAll[indB[j]] = valB[j];
687 
688  for (size_t j = 0; j < nnzA; j++) {
689  // The cost of the whole Frobenius dot product function depends on the
690  // cost of the getLocalElement and getGlobalElement functions here.
691  LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
692  if (ind != invalid)
693  f += valBAll[ind] * valA[j];
694  }
695 
696  // Clean up array values
697  for (size_t j = 0; j < nnzB; j++)
698  valBAll[indB[j]] = zero;
699  }
700 
701  MueLu_sumAll(AColMap.getComm(), f, gf);
702 
703  return gf;
704  }
705 
715  static void SetRandomSeed(const Teuchos::Comm<int> &comm) {
716  // Distribute the seeds evenly in [1,maxint-1]. This guarantees nothing
717  // about where in random number stream we are, but avoids overflow situations
718  // in parallel when multiplying by a PID. It would be better to use
719  // a good parallel random number generator.
720  double one = 1.0;
721  int maxint = INT_MAX; //= 2^31-1 = 2147483647 for 32-bit integers
722  int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.getRank()+1)/(comm.getSize()+one)) );
723  if (mySeed < 1 || mySeed == maxint) {
724  std::ostringstream errStr;
725  errStr << "Error detected with random seed = " << mySeed << ". It should be in the interval [1,2^31-2].";
726  throw Exceptions::RuntimeError(errStr.str());
727  }
728  std::srand(mySeed);
729  // For Tpetra, we could use Kokkos' random number generator here.
731  // Epetra
732  // MultiVector::Random() -> Epetra_Util::RandomDouble() -> Epetra_Utils::RandomInt()
733  // Its own random number generator, based on Seed_. Seed_ is initialized in Epetra_Util constructor with std::rand()
734  // So our setting std::srand() affects that too
735  }
736 
737 
738 
739  // Finds the OAZ Dirichlet rows for this matrix
740  // so far only used in IntrepidPCoarsenFactory
741  // TODO check whether we can use DetectDirichletRows instead
742  static void FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
743  std::vector<LocalOrdinal>& dirichletRows, bool count_twos_as_dirichlet=false) {
745  dirichletRows.resize(0);
746  for(size_t i=0; i<A->getNodeNumRows(); i++) {
749  A->getLocalRowView(i,indices,values);
750  int nnz=0;
751  for (size_t j=0; j<(size_t)indices.size(); j++) {
753  nnz++;
754  }
755  }
756  if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
757  dirichletRows.push_back(i);
758  }
759  }
760  }
761 
762  // Applies Ones-and-Zeros to matrix rows
763  // Takes a vector of row indices
764  static void ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
765  const std::vector<LocalOrdinal>& dirichletRows) {
766  RCP<const Map> Rmap = A->getColMap();
767  RCP<const Map> Cmap = A->getColMap();
770 
771  for(size_t i=0; i<dirichletRows.size(); i++) {
772  GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
773 
776  A->getLocalRowView(dirichletRows[i],indices,values);
777  // NOTE: This won't work with fancy node types.
778  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
779  for(size_t j=0; j<(size_t)indices.size(); j++) {
780  if(Cmap->getGlobalElement(indices[j])==row_gid)
781  valuesNC[j]=one;
782  else
783  valuesNC[j]=zero;
784  }
785  }
786  }
787 
788  // Applies Ones-and-Zeros to matrix rows
789  // Takes a Boolean array.
790  static void ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
791  const Teuchos::ArrayRCP<const bool>& dirichletRows) {
792  RCP<const Map> Rmap = A->getColMap();
793  RCP<const Map> Cmap = A->getColMap();
796 
797  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
798  if (dirichletRows[i]){
799  GlobalOrdinal row_gid = Rmap->getGlobalElement(i);
800 
803  A->getLocalRowView(i,indices,values);
804  // NOTE: This won't work with fancy node types.
805  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
806  for(size_t j=0; j<(size_t)indices.size(); j++) {
807  if(Cmap->getGlobalElement(indices[j])==row_gid)
808  valuesNC[j]=one;
809  else
810  valuesNC[j]=zero;
811  }
812  }
813  }
814  }
815 
816  // Zeros out rows
817  // Takes a vector containg Dirichlet row indices
818  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
819  const std::vector<LocalOrdinal>& dirichletRows,
820  Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
821  for(size_t i=0; i<dirichletRows.size(); i++) {
824  A->getLocalRowView(dirichletRows[i],indices,values);
825  // NOTE: This won't work with fancy node types.
826  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
827  for(size_t j=0; j<(size_t)indices.size(); j++)
828  valuesNC[j]=replaceWith;
829  }
830  }
831 
832  // Zeros out rows
833  // Takes a Boolean ArrayRCP
834  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
835  const Teuchos::ArrayRCP<const bool>& dirichletRows,
836  Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
837  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
838  if (dirichletRows[i]) {
841  A->getLocalRowView(i,indices,values);
842  // NOTE: This won't work with fancy node types.
843  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
844  for(size_t j=0; j<(size_t)indices.size(); j++)
845  valuesNC[j]=replaceWith;
846  }
847  }
848  }
849 
850  // Zeros out rows
851  // Takes a Boolean ArrayRCP
852  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& X,
853  const Teuchos::ArrayRCP<const bool>& dirichletRows,
854  Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
855  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
856  if (dirichletRows[i]) {
857  for(size_t j=0; j<X->getNumVectors(); j++)
858  X->replaceLocalValue(i,j,replaceWith);
859  }
860  }
861  }
862 
863  // Zeros out columns
864  // Takes a Boolean vector
866  const Teuchos::ArrayRCP<const bool>& dirichletCols,
867  Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
868  for(size_t i=0; i<A->getNodeNumRows(); i++) {
871  A->getLocalRowView(i,indices,values);
872  // NOTE: This won't work with fancy node types.
873  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
874  for(size_t j=0; j<static_cast<size_t>(indices.size()); j++)
875  if (dirichletCols[indices[j]])
876  valuesNC[j] = replaceWith;
877  }
878  }
879 
880  // Finds the OAZ Dirichlet rows for this matrix
881  static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
882  Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletRow,
883  Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletCol) {
884 
885  // Make sure A's RowMap == DomainMap
886  if(!A->getRowMap()->isSameAs(*A->getDomainMap())) {
887  throw std::runtime_error("UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
888  }
889  RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A->getCrsGraph()->getImporter();
890  bool has_import = !importer.is_null();
891 
892  // Find the Dirichlet rows
893  std::vector<LocalOrdinal> dirichletRows;
894  FindDirichletRows(A,dirichletRows);
895 
896 #if 0
897  printf("[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
898  for(size_t i=0; i<(size_t) dirichletRows.size(); i++)
899  printf("%d ",dirichletRows[i]);
900  printf("\n");
901  fflush(stdout);
902 #endif
903  // Allocate all as non-Dirichlet
904  isDirichletRow = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getRowMap(),true);
905  isDirichletCol = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getColMap(),true);
906 
907  Teuchos::ArrayRCP<int> dr_rcp = isDirichletRow->getDataNonConst(0);
908  Teuchos::ArrayView<int> dr = dr_rcp();
909  Teuchos::ArrayRCP<int> dc_rcp = isDirichletCol->getDataNonConst(0);
910  Teuchos::ArrayView<int> dc = dc_rcp();
911  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
912  dr[dirichletRows[i]] = 1;
913  if(!has_import) dc[dirichletRows[i]] = 1;
914  }
915 
916  if(has_import)
917  isDirichletCol->doImport(*isDirichletRow,*importer,Xpetra::CombineMode::ADD);
918 
919  }
920 
921  // This routine takes a BlockedMap and an Importer (assuming that the BlockedMap matches the source of the importer) and generates a BlockedMap corresponding
922  // to the Importer's target map. We assume that the targetMap is unique (which, is not a strict requirement of an Importer, but is here and no, we don't check)
923  // This is largely intended to be used in repartitioning of blocked matrices
924  static RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > GeneratedBlockedTargetMap(const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> & sourceBlockedMap,
925  const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> & Importer) {
926  typedef Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> IntVector;
927  Xpetra::UnderlyingLib lib = sourceBlockedMap.lib();
928 
929  // De-stride the map if we have to (might regret this later)
930  RCP<const Map> fullMap = sourceBlockedMap.getMap();
931  RCP<const Map> stridedMap = Teuchos::rcp_dynamic_cast<const Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >(fullMap);
932  if(!stridedMap.is_null()) fullMap = stridedMap->getMap();
933 
934  // Initial sanity checking for map compatibil
935  const size_t numSubMaps = sourceBlockedMap.getNumMaps();
936  if(!Importer.getSourceMap()->isCompatible(*fullMap))
937  throw std::runtime_error("GenerateBlockedTargetMap(): Map compatibility error");
938 
939  // Build an indicator vector
940  RCP<IntVector> block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(fullMap);
941 
942  for(size_t i=0; i<numSubMaps; i++) {
943  RCP<const Map> map = sourceBlockedMap.getMap(i);
944 
945  for(size_t j=0; j<map->getNodeNumElements(); j++) {
946  LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
947  block_ids->replaceLocalValue(jj,(int)i);
948  }
949  }
950 
951  // Get the block ids for the new map
952  RCP<const Map> targetMap = Importer.getTargetMap();
953  RCP<IntVector> new_block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(targetMap);
954  new_block_ids->doImport(*block_ids,Importer,Xpetra::CombineMode::ADD);
955  Teuchos::ArrayRCP<const int> dataRCP = new_block_ids->getData(0);
956  Teuchos::ArrayView<const int> data = dataRCP();
957 
958 
959  // Get the GIDs for each subblock
960  Teuchos::Array<Teuchos::Array<GlobalOrdinal> > elementsInSubMap(numSubMaps);
961  for(size_t i=0; i<targetMap->getNodeNumElements(); i++) {
962  elementsInSubMap[data[i]].push_back(targetMap->getGlobalElement(i));
963  }
964 
965  // Generate the new submaps
966  std::vector<RCP<const Map> > subMaps(numSubMaps);
967  for(size_t i=0; i<numSubMaps; i++) {
968  subMaps[i] = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(lib,Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),elementsInSubMap[i](),targetMap->getIndexBase(),targetMap->getComm(),targetMap->getNode());
969  }
970 
971  // Build the BlockedMap
972  return rcp(new BlockedMap(targetMap,subMaps));
973 
974  }
975 
976  }; // class Utils
977 
978 
979 
981 
982 } //namespace MueLu
983 
984 #define MUELU_UTILITIESBASE_SHORT
985 #endif // MUELU_UTILITIESBASE_DECL_HPP
Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedMultiVector
static Teuchos::ArrayRCP< const bool > DetectDirichletCols(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
static magnitudeType eps()
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.
static void Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS, MultiVector &Resid)
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MultiVector
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows)
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 void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
size_type size() const
magnitude_type tolerance
Namespace for MueLu classes and methods.
static Teuchos::RCP< const Comm< OrdinalType > > getComm()
Exception throws to report incompatible objects (like maps).
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows (extended version)
#define MueLu_sumAll(rcpComm, in, out)
bool is_null() const
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > > &isDirichletRow, Teuchos::RCP< Xpetra::Vector< int, LocalOrdinal, GlobalOrdinal, Node > > &isDirichletCol)
virtual int getRank() const=0
Xpetra::BlockedVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedVector
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 void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
static magnitudeType magnitude(T a)
Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > BlockedMap
static void ZeroDirichletCols(Teuchos::RCP< Matrix > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
void push_back(const value_type &x)
T * getRawPtr() const
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
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)
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::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS, MultiVector &Resid)
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Teuchos::RCP< const Matrix > rcpA)
Extract Matrix Diagonal of lumped matrix.
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &X, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
virtual int getSize() const=0