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>
59 #include <Teuchos_ParameterList.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 namespace MueLu {
88 
89 // MPI helpers
90 #define MueLu_sumAll(rcpComm, in, out) \
91  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_SUM, in, Teuchos::outArg(out))
92 #define MueLu_minAll(rcpComm, in, out) \
93  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MIN, in, Teuchos::outArg(out))
94 #define MueLu_maxAll(rcpComm, in, out) \
95  Teuchos::reduceAll(*rcpComm, Teuchos::REDUCE_MAX, in, Teuchos::outArg(out))
96 
104  template <class Scalar,
107  class Node = DefaultNode>
108  class UtilitiesBase {
109  public:
110 #undef MUELU_UTILITIESBASE_SHORT
111 //#include "MueLu_UseShortNames.hpp"
112  private:
113  typedef Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> CrsMatrixWrap;
114  typedef Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> CrsMatrix;
115  typedef Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Matrix;
116  typedef Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Vector;
117  typedef Xpetra::BlockedVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> BlockedVector;
118  typedef Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MultiVector;
119  typedef Xpetra::BlockedMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> BlockedMultiVector;
120  typedef Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> BlockedMap;
121  typedef Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
122  public:
123  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
124 
125 
126  static RCP<Matrix> Crs2Op(RCP<CrsMatrix> Op) {
127  if (Op.is_null())
128  return Teuchos::null;
129  return rcp(new CrsMatrixWrap(Op));
130  }
131 
138  static Teuchos::ArrayRCP<Scalar> GetMatrixDiagonal(const Matrix& A) {
139  size_t numRows = A.getRowMap()->getNodeNumElements();
140  Teuchos::ArrayRCP<Scalar> diag(numRows);
141  Teuchos::ArrayView<const LocalOrdinal> cols;
142  Teuchos::ArrayView<const Scalar> vals;
143  for (size_t i = 0; i < numRows; ++i) {
144  A.getLocalRowView(i, cols, vals);
145  LocalOrdinal j = 0;
146  for (; j < cols.size(); ++j) {
147  if (Teuchos::as<size_t>(cols[j]) == i) {
148  diag[i] = vals[j];
149  break;
150  }
151  }
152  if (j == cols.size()) {
153  // Diagonal entry is absent
154  diag[i] = Teuchos::ScalarTraits<Scalar>::zero();
155  }
156  }
157  return diag;
158  }
159 
166  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
167  Teuchos::TimeMonitor MM = *Teuchos::TimeMonitor::getNewTimer("UtilitiesBase::GetMatrixDiagonalInverse");
168 
169  RCP<const Map> rowMap = A.getRowMap();
170  RCP<Vector> diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
171 
172  A.getLocalDiagCopy(*diag);
173 
174  RCP<Vector> inv = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(diag, tol, valReplacement);
175 
176  return inv;
177  }
178 
185  static Teuchos::RCP<Vector> GetLumpedMatrixDiagonal(Matrix const & A, const bool doReciprocal = false,
186  Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100,
187  Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero(),
188  const bool replaceSingleEntryRowWithZero = false) {
189 
190  typedef Teuchos::ScalarTraits<Scalar> TST;
191 
192  RCP<Vector> diag = Teuchos::null;
193  const Scalar zero = TST::zero();
194  const Scalar one = TST::one();
195  const Scalar two = one + one;
196 
197 tol = 0.;
198 
199  Teuchos::RCP<const Matrix> rcpA = Teuchos::rcpFromRef(A);
200 
201  RCP<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bA =
202  Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(rcpA);
203  if(bA == Teuchos::null) {
204  RCP<const Map> rowMap = rcpA->getRowMap();
205  diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
206  ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
207  Teuchos::Array<Scalar> regSum(diag->getLocalLength());
208  Teuchos::ArrayView<const LocalOrdinal> cols;
209  Teuchos::ArrayView<const Scalar> vals;
210 
211  std::vector<int> nnzPerRow(rowMap->getNodeNumElements());
212 
213  const Magnitude zeroMagn = TST::magnitude(zero);
214  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
215  nnzPerRow[i] = 0;
216  rcpA->getLocalRowView(i, cols, vals);
217  diagVals[i] = zero;
218  for (LocalOrdinal j = 0; j < cols.size(); ++j) {
219  regSum[i] += vals[j];
220  const Magnitude rowEntryMagn = TST::magnitude(vals[j]);
221  if (rowEntryMagn > zeroMagn)
222  nnzPerRow[i]++;
223  diagVals[i] += rowEntryMagn;
224  }
225  }
226  if (doReciprocal) {
227  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
228  if (replaceSingleEntryRowWithZero && nnzPerRow[i] <= static_cast<int>(1))
229  diagVals[i] = zero;
230  else if (replaceSingleEntryRowWithZero && diagVals[i] != zero && TST::magnitude(diagVals[i]) < TST::magnitude(two*regSum[i]))
231  diagVals[i] = one / (two*regSum[i]);
232  else {
233  if(TST::magnitude(diagVals[i]) > tol)
234  diagVals[i] = one / diagVals[i];
235  else {
236  diagVals[i] = valReplacement;
237  }
238  }
239  }
240  }
241  } else {
242  TEUCHOS_TEST_FOR_EXCEPTION(doReciprocal, Xpetra::Exceptions::RuntimeError,
243  "UtilitiesBase::GetLumpedMatrixDiagonal(): extracting reciprocal of diagonal of a blocked matrix is not supported");
244  diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(bA->getRangeMapExtractor()->getFullMap(),true);
245 
246  for (size_t row = 0; row < bA->Rows(); ++row) {
247  for (size_t col = 0; col < bA->Cols(); ++col) {
248  if (!bA->getMatrix(row,col).is_null()) {
249  // 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!
250  bool bThyraMode = bA->getRangeMapExtractor()->getThyraMode() && (Teuchos::rcp_dynamic_cast<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA->getMatrix(row,col)) == Teuchos::null);
251  RCP<Vector> ddtemp = bA->getRangeMapExtractor()->ExtractVector(diag,row,bThyraMode);
252  RCP<const Vector> dd = GetLumpedMatrixDiagonal(*(bA->getMatrix(row,col)));
253  ddtemp->update(Teuchos::as<Scalar>(1.0),*dd,Teuchos::as<Scalar>(1.0));
254  bA->getRangeMapExtractor()->InsertVector(ddtemp,row,diag,bThyraMode);
255  }
256  }
257  }
258 
259  }
260 
261  return diag;
262  }
263 
270  static Teuchos::ArrayRCP<Magnitude> GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A) {
271  size_t numRows = A.getRowMap()->getNodeNumElements();
272  Magnitude ZERO = Teuchos::ScalarTraits<Magnitude>::zero();
273  Teuchos::ArrayRCP<Magnitude> maxvec(numRows);
274  Teuchos::ArrayView<const LocalOrdinal> cols;
275  Teuchos::ArrayView<const Scalar> vals;
276  for (size_t i = 0; i < numRows; ++i) {
277  A.getLocalRowView(i, cols, vals);
278  Magnitude mymax = ZERO;
279  for (LocalOrdinal j=0; j < cols.size(); ++j) {
280  if (Teuchos::as<size_t>(cols[j]) != i) {
281  mymax = std::max(mymax,-Teuchos::ScalarTraits<Scalar>::real(vals[j]));
282  }
283  }
284  maxvec[i] = mymax;
285  }
286  return maxvec;
287  }
288 
289  static Teuchos::ArrayRCP<Magnitude> GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> &BlockNumber) {
290  TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()),std::runtime_error,"GetMatrixMaxMinusOffDiagonal: BlockNumber must match's A's column map.");
291 
292  Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
293 
294  size_t numRows = A.getRowMap()->getNodeNumElements();
295  Magnitude ZERO = Teuchos::ScalarTraits<Magnitude>::zero();
296  Teuchos::ArrayRCP<Magnitude> maxvec(numRows);
297  Teuchos::ArrayView<const LocalOrdinal> cols;
298  Teuchos::ArrayView<const Scalar> vals;
299  for (size_t i = 0; i < numRows; ++i) {
300  A.getLocalRowView(i, cols, vals);
301  Magnitude mymax = ZERO;
302  for (LocalOrdinal j=0; j < cols.size(); ++j) {
303  if (Teuchos::as<size_t>(cols[j]) != i && block_id[i] == block_id[cols[j]]) {
304  mymax = std::max(mymax,-Teuchos::ScalarTraits<Scalar>::real(vals[j]));
305  }
306  }
307  // printf("A(%d,:) row_scale(block) = %6.4e\n",(int)i,mymax);
308 
309  maxvec[i] = mymax;
310  }
311  return maxvec;
312  }
313 
321  static Teuchos::RCP<Vector> GetInverse(Teuchos::RCP<const Vector> v, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar valReplacement = Teuchos::ScalarTraits<Scalar>::zero()) {
322 
323  RCP<Vector> ret = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(v->getMap(),true);
324 
325  // check whether input vector "v" is a BlockedVector
326  RCP<const BlockedVector> bv = Teuchos::rcp_dynamic_cast<const BlockedVector>(v);
327  if(bv.is_null() == false) {
328  RCP<BlockedVector> bret = Teuchos::rcp_dynamic_cast<BlockedVector>(ret);
329  TEUCHOS_TEST_FOR_EXCEPTION(bret.is_null() == true, MueLu::Exceptions::RuntimeError,"MueLu::UtilitiesBase::GetInverse: return vector should be of type BlockedVector");
330  RCP<const BlockedMap> bmap = bv->getBlockedMap();
331  for(size_t r = 0; r < bmap->getNumMaps(); ++r) {
332  RCP<const MultiVector> submvec = bv->getMultiVector(r,bmap->getThyraMode());
333  RCP<const Vector> subvec = submvec->getVector(0);
334  RCP<Vector> subvecinf = MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(subvec,tol,valReplacement);
335  bret->setMultiVector(r, subvecinf, bmap->getThyraMode());
336  }
337  return ret;
338  }
339 
340  // v is an {Epetra,Tpetra}Vector: work with the underlying raw data
341  ArrayRCP<Scalar> retVals = ret->getDataNonConst(0);
342  ArrayRCP<const Scalar> inputVals = v->getData(0);
343  for (size_t i = 0; i < v->getMap()->getNodeNumElements(); ++i) {
344  if(Teuchos::ScalarTraits<Scalar>::magnitude(inputVals[i]) > tol)
345  retVals[i] = Teuchos::ScalarTraits<Scalar>::one() / inputVals[i];
346  else
347  retVals[i] = valReplacement;
348  }
349  return ret;
350  }
351 
359  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A) {
360  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
361 
362  // Undo block map (if we have one)
363  RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
364  if(!browMap.is_null()) rowMap = browMap->getMap();
365 
366  RCP<Vector> localDiag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap);
367  try {
368  const CrsMatrixWrap* crsOp = dynamic_cast<const CrsMatrixWrap*>(&A);
369  if (crsOp == NULL) {
370  throw Exceptions::RuntimeError("cast to CrsMatrixWrap failed");
371  }
372  Teuchos::ArrayRCP<size_t> offsets;
373  crsOp->getLocalDiagOffsets(offsets);
374  crsOp->getLocalDiagCopy(*localDiag,offsets());
375  }
376  catch (...) {
377  ArrayRCP<Scalar> localDiagVals = localDiag->getDataNonConst(0);
378  Teuchos::ArrayRCP<Scalar> diagVals = GetMatrixDiagonal(A);
379  for (LocalOrdinal i = 0; i < localDiagVals.size(); i++)
380  localDiagVals[i] = diagVals[i];
381  localDiagVals = diagVals = null;
382  }
383 
384  RCP<Vector> diagonal = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap);
385  RCP< const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer;
386  importer = A.getCrsGraph()->getImporter();
387  if (importer == Teuchos::null) {
388  importer = Xpetra::ImportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap, colMap);
389  }
390  diagonal->doImport(*localDiag, *(importer), Xpetra::INSERT);
391  return diagonal;
392  }
393 
394 
402  static RCP<Vector> GetMatrixOverlappedDeletedRowsum(const Matrix& A) {
403  using LO = LocalOrdinal;
404  using GO = GlobalOrdinal;
405  using SC = Scalar;
406  using STS = typename Teuchos::ScalarTraits<SC>;
407 
408  // Undo block map (if we have one)
409  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
410  RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
411  if(!browMap.is_null()) rowMap = browMap->getMap();
412 
413  RCP<Vector> local = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(rowMap);
414  RCP<Vector> ghosted = Xpetra::VectorFactory<SC,LO,GO,Node>::Build(colMap,true);
415  ArrayRCP<SC> localVals = local->getDataNonConst(0);
416 
417  for (LO row = 0; row < static_cast<LO>(A.getRowMap()->getNodeNumElements()); ++row) {
418  size_t nnz = A.getNumEntriesInLocalRow(row);
419  ArrayView<const LO> indices;
420  ArrayView<const SC> vals;
421  A.getLocalRowView(row, indices, vals);
422 
423  SC si = STS::zero();
424 
425  for (LO colID = 0; colID < static_cast<LO>(nnz); colID++) {
426  if(indices[colID] != row) {
427  si += vals[colID];
428  }
429  }
430  localVals[row] = si;
431  }
432 
433  RCP< const Xpetra::Import<LO,GO,Node> > importer;
434  importer = A.getCrsGraph()->getImporter();
435  if (importer == Teuchos::null) {
436  importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
437  }
438  ghosted->doImport(*local, *(importer), Xpetra::INSERT);
439  return ghosted;
440  }
441 
442 
443 
444  static RCP<Xpetra::Vector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> >
446  RCP<const Map> rowMap = A.getRowMap(), colMap = A.getColMap();
447  using STS = typename Teuchos::ScalarTraits<Scalar>;
448  using MTS = typename Teuchos::ScalarTraits<Magnitude>;
449  using MT = Magnitude;
450  using LO = LocalOrdinal;
451  using GO = GlobalOrdinal;
452  using SC = Scalar;
453  using RealValuedVector = Xpetra::Vector<MT,LO,GO,Node>;
454 
455  // Undo block map (if we have one)
456  RCP<const BlockedMap> browMap = Teuchos::rcp_dynamic_cast<const BlockedMap>(rowMap);
457  if(!browMap.is_null()) rowMap = browMap->getMap();
458 
459  RCP<RealValuedVector> local = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(rowMap);
460  RCP<RealValuedVector> ghosted = Xpetra::VectorFactory<MT,LO,GO,Node>::Build(colMap,true);
461  ArrayRCP<MT> localVals = local->getDataNonConst(0);
462 
463  for (LO rowIdx = 0; rowIdx < static_cast<LO>(A.getRowMap()->getNodeNumElements()); ++rowIdx) {
464  size_t nnz = A.getNumEntriesInLocalRow(rowIdx);
465  ArrayView<const LO> indices;
466  ArrayView<const SC> vals;
467  A.getLocalRowView(rowIdx, indices, vals);
468 
469  MT si = MTS::zero();
470 
471  for (LO colID = 0; colID < static_cast<LO>(nnz); ++colID) {
472  if(indices[colID] != rowIdx) {
473  si += STS::magnitude(vals[colID]);
474  }
475  }
476  localVals[rowIdx] = si;
477  }
478 
479  RCP< const Xpetra::Import<LO,GO,Node> > importer;
480  importer = A.getCrsGraph()->getImporter();
481  if (importer == Teuchos::null) {
482  importer = Xpetra::ImportFactory<LO,GO,Node>::Build(rowMap, colMap);
483  }
484  ghosted->doImport(*local, *(importer), Xpetra::INSERT);
485  return ghosted;
486  }
487 
488 
489 
490  // TODO: should NOT return an Array. Definition must be changed to:
491  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
492  // or
493  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
494  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS) {
495  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
496  const size_t numVecs = X.getNumVectors();
497  RCP<MultiVector> RES = Residual(Op, X, RHS);
498  Teuchos::Array<Magnitude> norms(numVecs);
499  RES->norm2(norms);
500  return norms;
501  }
502 
503  static Teuchos::Array<Magnitude> ResidualNorm(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector & Resid) {
504  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
505  const size_t numVecs = X.getNumVectors();
506  Residual(Op,X,RHS,Resid);
507  Teuchos::Array<Magnitude> norms(numVecs);
508  Resid.norm2(norms);
509  return norms;
510  }
511 
512  static RCP<MultiVector> Residual(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS) {
513  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides")
514  const size_t numVecs = X.getNumVectors();
515  // TODO Op.getRangeMap should return a BlockedMap if it is a BlockedCrsOperator
516  RCP<MultiVector> RES = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(RHS.getMap(), numVecs, false); // no need to initialize to zero
517  Op.residual(X,RHS,*RES);
518  return RES;
519  }
520 
521 
522  static void Residual(const Xpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const MultiVector& X, const MultiVector& RHS, MultiVector & Resid) {
523  TEUCHOS_TEST_FOR_EXCEPTION(X.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of solution vectors != number of right-hand sides");
524  TEUCHOS_TEST_FOR_EXCEPTION(Resid.getNumVectors() != RHS.getNumVectors(), Exceptions::RuntimeError, "Number of residual vectors != number of right-hand sides");
525  Op.residual(X,RHS,Resid);
526  }
527 
528 
540  static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true,
541  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
542  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
543  "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
544 
545  // power iteration
546  RCP<Vector> diagInvVec;
547  if (scaleByDiag) {
548  RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
549  A.getLocalDiagCopy(*diagVec);
550  diagInvVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
551  diagInvVec->reciprocal(*diagVec);
552  }
553 
554  Scalar lambda = PowerMethod(A, diagInvVec, niters, tolerance, verbose, seed);
555  return lambda;
556  }
557 
569  static Scalar PowerMethod(const Matrix& A, const RCP<Vector> &diagInvVec,
570  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
571  TEUCHOS_TEST_FOR_EXCEPTION(!(A.getRangeMap()->isSameAs(*(A.getDomainMap()))), Exceptions::Incompatible,
572  "Utils::PowerMethod: operator must have domain and range maps that are equivalent.");
573 
574  // Create three vectors, fill z with random numbers
575  RCP<Vector> q = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getDomainMap());
576  RCP<Vector> r = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
577  RCP<Vector> z = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRangeMap());
578 
579  z->setSeed(seed); // seed random number generator
580  z->randomize(true);// use Xpetra implementation: -> same results for Epetra and Tpetra
581 
582  Teuchos::Array<Magnitude> norms(1);
583 
584  typedef Teuchos::ScalarTraits<Scalar> STS;
585 
586  const Scalar zero = STS::zero(), one = STS::one();
587 
588  Scalar lambda = zero;
589  Magnitude residual = STS::magnitude(zero);
590 
591  // power iteration
592  for (int iter = 0; iter < niters; ++iter) {
593  z->norm2(norms); // Compute 2-norm of z
594  q->update(one/norms[0], *z, zero); // Set q = z / normz
595  A.apply(*q, *z); // Compute z = A*q
596  if (diagInvVec != Teuchos::null)
597  z->elementWiseMultiply(one, *diagInvVec, *z, zero);
598  lambda = q->dot(*z); // Approximate maximum eigenvalue: lamba = dot(q,z)
599 
600  if (iter % 100 == 0 || iter + 1 == niters) {
601  r->update(1.0, *z, -lambda, *q, zero); // Compute A*q - lambda*q
602  r->norm2(norms);
603  residual = STS::magnitude(norms[0] / lambda);
604  if (verbose) {
605  std::cout << "Iter = " << iter
606  << " Lambda = " << lambda
607  << " Residual of A*q - lambda*q = " << residual
608  << std::endl;
609  }
610  }
611  if (residual < tolerance)
612  break;
613  }
614  return lambda;
615  }
616 
617  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) {
618  RCP<Teuchos::FancyOStream> fancy = Teuchos::fancyOStream(Teuchos::rcpFromRef(os));
619  return fancy;
620  }
621 
626  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>> &v, LocalOrdinal i0, LocalOrdinal i1) {
627  const size_t numVectors = v.size();
628 
629  Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
630  for (size_t j = 0; j < numVectors; j++) {
631  d += (v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
632  }
633  return Teuchos::ScalarTraits<Scalar>::magnitude(d);
634  }
635 
640  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Teuchos::ArrayView<double> & weight,const Teuchos::Array<Teuchos::ArrayRCP<const Scalar>> &v, LocalOrdinal i0, LocalOrdinal i1) {
641  const size_t numVectors = v.size();
642 
643  Scalar d = Teuchos::ScalarTraits<Scalar>::zero();
644  for (size_t j = 0; j < numVectors; j++) {
645  d += weight[j]*(v[j][i0] - v[j][i1])*(v[j][i0] - v[j][i1]);
646  }
647  return Teuchos::ScalarTraits<Scalar>::magnitude(d);
648  }
649 
650 
663  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) {
664  LocalOrdinal numRows = A.getNodeNumRows();
665  typedef Teuchos::ScalarTraits<Scalar> STS;
666  ArrayRCP<bool> boundaryNodes(numRows, true);
667  if (count_twos_as_dirichlet) {
668  for (LocalOrdinal row = 0; row < numRows; row++) {
669  ArrayView<const LocalOrdinal> indices;
670  ArrayView<const Scalar> vals;
671  A.getLocalRowView(row, indices, vals);
672  size_t nnz = A.getNumEntriesInLocalRow(row);
673  if (nnz > 2) {
674  size_t col;
675  for (col = 0; col < nnz; col++)
676  if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
677  if (!boundaryNodes[row])
678  break;
679  boundaryNodes[row] = false;
680  }
681  if (col == nnz)
682  boundaryNodes[row] = true;
683  }
684  }
685  } else {
686  for (LocalOrdinal row = 0; row < numRows; row++) {
687  ArrayView<const LocalOrdinal> indices;
688  ArrayView<const Scalar> vals;
689  A.getLocalRowView(row, indices, vals);
690  size_t nnz = A.getNumEntriesInLocalRow(row);
691  if (nnz > 1)
692  for (size_t col = 0; col < nnz; col++)
693  if ( (indices[col] != row) && STS::magnitude(vals[col]) > tol) {
694  boundaryNodes[row] = false;
695  break;
696  }
697  }
698  }
699  return boundaryNodes;
700  }
701 
702 
715  static Teuchos::ArrayRCP<const bool> DetectDirichletRowsExt(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, bool & bHasZeroDiagonal, const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero()) {
716 
717  // assume that there is no zero diagonal in matrix
718  bHasZeroDiagonal = false;
719 
720  Teuchos::RCP<Vector> diagVec = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(A.getRowMap());
721  A.getLocalDiagCopy(*diagVec);
722  Teuchos::ArrayRCP< const Scalar > diagVecData = diagVec->getData(0);
723 
724  LocalOrdinal numRows = A.getNodeNumRows();
725  typedef Teuchos::ScalarTraits<Scalar> STS;
726  ArrayRCP<bool> boundaryNodes(numRows, false);
727  for (LocalOrdinal row = 0; row < numRows; row++) {
728  ArrayView<const LocalOrdinal> indices;
729  ArrayView<const Scalar> vals;
730  A.getLocalRowView(row, indices, vals);
731  size_t nnz = 0; // collect nonzeros in row (excluding the diagonal)
732  bool bHasDiag = false;
733  for (decltype(indices.size()) col = 0; col < indices.size(); col++) {
734  if ( indices[col] != row) {
735  if (STS::magnitude(vals[col] / STS::magnitude(sqrt(STS::magnitude(diagVecData[row]) * STS::magnitude(diagVecData[col]))) ) > tol) {
736  nnz++;
737  }
738  } else bHasDiag = true; // found a diagonal entry
739  }
740  if (bHasDiag == false) bHasZeroDiagonal = true; // we found at least one row without a diagonal
741  else if(nnz == 0) boundaryNodes[row] = true;
742  }
743  return boundaryNodes;
744  }
745 
753  static void FindNonZeros(const Teuchos::ArrayRCP<const Scalar> vals,
754  Teuchos::ArrayRCP<bool> nonzeros) {
755  TEUCHOS_ASSERT(vals.size() == nonzeros.size());
756  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
757  const magnitudeType eps = 2.0*Teuchos::ScalarTraits<magnitudeType>::eps();
758  for(size_t i=0; i<static_cast<size_t>(vals.size()); i++) {
759  nonzeros[i] = (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > eps);
760  }
761  }
762 
771  static void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
772  const Teuchos::ArrayRCP<bool>& dirichletRows,
773  Teuchos::ArrayRCP<bool> dirichletCols,
774  Teuchos::ArrayRCP<bool> dirichletDomain) {
775  const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
776  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A .getDomainMap();
777  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
778  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
779  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == rowMap->getNodeNumElements());
780  TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == colMap->getNodeNumElements());
781  TEUCHOS_ASSERT(static_cast<size_t>(dirichletDomain.size()) == domMap->getNodeNumElements());
782  RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, 1, /*zeroOut=*/true);
783  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
784  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
785  if (dirichletRows[i]) {
786  ArrayView<const LocalOrdinal> indices;
787  ArrayView<const Scalar> values;
788  A.getLocalRowView(i,indices,values);
789  for(size_t j=0; j<static_cast<size_t>(indices.size()); j++)
790  myColsToZero->replaceLocalValue(indices[j],0,one);
791  }
792  }
793 
794  RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
795  RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
796  if (!importer.is_null()) {
797  globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, 1, /*zeroOut=*/true);
798  // export to domain map
799  globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
800  // import to column map
801  myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
802  }
803  else
804  globalColsToZero = myColsToZero;
805 
806  FindNonZeros(globalColsToZero->getData(0),dirichletDomain);
807  FindNonZeros(myColsToZero->getData(0),dirichletCols);
808  }
809 
810 
811 
823  static void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
824  typedef Teuchos::ScalarTraits<Scalar> STS;
825  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
826  typedef Teuchos::ScalarTraits<MT> MTS;
827  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>> rowmap = A.getRowMap();
828  for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getNodeNumElements()); ++row) {
829  size_t nnz = A.getNumEntriesInLocalRow(row);
830  ArrayView<const LocalOrdinal> indices;
831  ArrayView<const Scalar> vals;
832  A.getLocalRowView(row, indices, vals);
833 
834  Scalar rowsum = STS::zero();
835  Scalar diagval = STS::zero();
836 
837  for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
838  LocalOrdinal col = indices[colID];
839  if (row == col)
840  diagval = vals[colID];
841  rowsum += vals[colID];
842  }
843  // printf("A(%d,:) row_sum(point) = %6.4e\n",row,rowsum);
844  if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
845  //printf("Row %d triggers rowsum\n",(int)row);
846  dirichletRows[row] = true;
847  }
848  }
849  }
850 
851  static void ApplyRowSumCriterion(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> &BlockNumber, const Magnitude rowSumTol, Teuchos::ArrayRCP<bool>& dirichletRows) {
852  typedef Teuchos::ScalarTraits<Scalar> STS;
853  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
854  typedef Teuchos::ScalarTraits<MT> MTS;
855  RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowmap = A.getRowMap();
856 
857  TEUCHOS_TEST_FOR_EXCEPTION(!A.getColMap()->isSameAs(*BlockNumber.getMap()),std::runtime_error,"ApplyRowSumCriterion: BlockNumber must match's A's column map.");
858 
859  Teuchos::ArrayRCP<const LocalOrdinal> block_id = BlockNumber.getData(0);
860  for (LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(rowmap->getNodeNumElements()); ++row) {
861  size_t nnz = A.getNumEntriesInLocalRow(row);
862  ArrayView<const LocalOrdinal> indices;
863  ArrayView<const Scalar> vals;
864  A.getLocalRowView(row, indices, vals);
865 
866  Scalar rowsum = STS::zero();
867  Scalar diagval = STS::zero();
868  for (LocalOrdinal colID = 0; colID < Teuchos::as<LocalOrdinal>(nnz); colID++) {
869  LocalOrdinal col = indices[colID];
870  if (row == col)
871  diagval = vals[colID];
872  if(block_id[row] == block_id[col])
873  rowsum += vals[colID];
874  }
875 
876  // printf("A(%d,:) row_sum(block) = %6.4e\n",row,rowsum);
877  if (rowSumTol < MTS::one() && STS::magnitude(rowsum) > STS::magnitude(diagval) * rowSumTol) {
878  //printf("Row %d triggers rowsum\n",(int)row);
879  dirichletRows[row] = true;
880  }
881  }
882  }
883 
884 
885 
896  static Teuchos::ArrayRCP<const bool> DetectDirichletCols(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
897  const Teuchos::ArrayRCP<const bool>& dirichletRows) {
898  Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
899  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
900  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
901  Teuchos::RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
902  Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,1);
903  myColsToZero->putScalar(zero);
904  // Find all local column indices that are in Dirichlet rows, record in myColsToZero as 1.0
905  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
906  if (dirichletRows[i]) {
907  Teuchos::ArrayView<const LocalOrdinal> indices;
908  Teuchos::ArrayView<const Scalar> values;
909  A.getLocalRowView(i,indices,values);
910  for(size_t j=0; j<static_cast<size_t>(indices.size()); j++)
911  myColsToZero->replaceLocalValue(indices[j],0,one);
912  }
913  }
914 
915  Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,1);
916  globalColsToZero->putScalar(zero);
917  Teuchos::RCP<Xpetra::Export<LocalOrdinal,GlobalOrdinal,Node> > exporter = Xpetra::ExportFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,domMap);
918  // export to domain map
919  globalColsToZero->doExport(*myColsToZero,*exporter,Xpetra::ADD);
920  // import to column map
921  myColsToZero->doImport(*globalColsToZero,*exporter,Xpetra::INSERT);
922  Teuchos::ArrayRCP<const Scalar> myCols = myColsToZero->getData(0);
923  Teuchos::ArrayRCP<bool> dirichletCols(colMap->getNodeNumElements(), true);
924  Magnitude eps = Teuchos::ScalarTraits<Magnitude>::eps();
925  for(size_t i=0; i<colMap->getNodeNumElements(); i++) {
926  dirichletCols[i] = Teuchos::ScalarTraits<Scalar>::magnitude(myCols[i])>2.0*eps;
927  }
928  return dirichletCols;
929  }
930 
931 
936  static Scalar Frobenius(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& B) {
937  // We check only row maps. Column may be different. One would hope that they are the same, as we typically
938  // calculate frobenius norm of the specified sparsity pattern with an updated matrix from the previous step,
939  // but matrix addition, even when one is submatrix of the other, changes column map (though change may be as
940  // simple as couple of elements swapped)
941  TEUCHOS_TEST_FOR_EXCEPTION(!A.getRowMap()->isSameAs(*B.getRowMap()), Exceptions::Incompatible, "MueLu::CGSolver::Frobenius: row maps are incompatible");
942  TEUCHOS_TEST_FOR_EXCEPTION(!A.isFillComplete() || !B.isFillComplete(), Exceptions::RuntimeError, "Matrices must be fill completed");
943 
944  const Map& AColMap = *A.getColMap();
945  const Map& BColMap = *B.getColMap();
946 
947  Teuchos::ArrayView<const LocalOrdinal> indA, indB;
948  Teuchos::ArrayView<const Scalar> valA, valB;
949  size_t nnzA = 0, nnzB = 0;
950 
951  // We use a simple algorithm
952  // for each row we fill valBAll array with the values in the corresponding row of B
953  // as such, it serves as both sorted array and as storage, so we don't need to do a
954  // tricky problem: "find a value in the row of B corresponding to the specific GID"
955  // Once we do that, we translate LID of entries of row of A to LID of B, and multiply
956  // corresponding entries.
957  // The algorithm should be reasonably cheap, as it does not sort anything, provided
958  // that getLocalElement and getGlobalElement functions are reasonably effective. It
959  // *is* possible that the costs are hidden in those functions, but if maps are close
960  // to linear maps, we should be fine
961  Teuchos::Array<Scalar> valBAll(BColMap.getNodeNumElements());
962 
963  LocalOrdinal invalid = Teuchos::OrdinalTraits<LocalOrdinal>::invalid();
964  Scalar zero = Teuchos::ScalarTraits<Scalar> ::zero(), f = zero, gf;
965  size_t numRows = A.getNodeNumRows();
966  for (size_t i = 0; i < numRows; i++) {
967  A.getLocalRowView(i, indA, valA);
968  B.getLocalRowView(i, indB, valB);
969  nnzA = indA.size();
970  nnzB = indB.size();
971 
972  // Set up array values
973  for (size_t j = 0; j < nnzB; j++)
974  valBAll[indB[j]] = valB[j];
975 
976  for (size_t j = 0; j < nnzA; j++) {
977  // The cost of the whole Frobenius dot product function depends on the
978  // cost of the getLocalElement and getGlobalElement functions here.
979  LocalOrdinal ind = BColMap.getLocalElement(AColMap.getGlobalElement(indA[j]));
980  if (ind != invalid)
981  f += valBAll[ind] * valA[j];
982  }
983 
984  // Clean up array values
985  for (size_t j = 0; j < nnzB; j++)
986  valBAll[indB[j]] = zero;
987  }
988 
989  MueLu_sumAll(AColMap.getComm(), f, gf);
990 
991  return gf;
992  }
993 
1003  static void SetRandomSeed(const Teuchos::Comm<int> &comm) {
1004  // Distribute the seeds evenly in [1,maxint-1]. This guarantees nothing
1005  // about where in random number stream we are, but avoids overflow situations
1006  // in parallel when multiplying by a PID. It would be better to use
1007  // a good parallel random number generator.
1008  double one = 1.0;
1009  int maxint = INT_MAX; //= 2^31-1 = 2147483647 for 32-bit integers
1010  int mySeed = Teuchos::as<int>((maxint-1) * (one -(comm.getRank()+1)/(comm.getSize()+one)) );
1011  if (mySeed < 1 || mySeed == maxint) {
1012  std::ostringstream errStr;
1013  errStr << "Error detected with random seed = " << mySeed << ". It should be in the interval [1,2^31-2].";
1014  throw Exceptions::RuntimeError(errStr.str());
1015  }
1016  std::srand(mySeed);
1017  // For Tpetra, we could use Kokkos' random number generator here.
1018  Teuchos::ScalarTraits<Scalar>::seedrandom(mySeed);
1019  // Epetra
1020  // MultiVector::Random() -> Epetra_Util::RandomDouble() -> Epetra_Utils::RandomInt()
1021  // Its own random number generator, based on Seed_. Seed_ is initialized in Epetra_Util constructor with std::rand()
1022  // So our setting std::srand() affects that too
1023  }
1024 
1025 
1026 
1027  // Finds the OAZ Dirichlet rows for this matrix
1028  // so far only used in IntrepidPCoarsenFactory
1029  // TODO check whether we can use DetectDirichletRows instead
1030  static void FindDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
1031  std::vector<LocalOrdinal>& dirichletRows, bool count_twos_as_dirichlet=false) {
1032  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType MT;
1033  dirichletRows.resize(0);
1034  for(size_t i=0; i<A->getNodeNumRows(); i++) {
1035  Teuchos::ArrayView<const LocalOrdinal> indices;
1036  Teuchos::ArrayView<const Scalar> values;
1037  A->getLocalRowView(i,indices,values);
1038  int nnz=0;
1039  for (size_t j=0; j<(size_t)indices.size(); j++) {
1040  if (Teuchos::ScalarTraits<Scalar>::magnitude(values[j]) > Teuchos::ScalarTraits<MT>::eps()) {
1041  nnz++;
1042  }
1043  }
1044  if (nnz == 1 || (count_twos_as_dirichlet && nnz == 2)) {
1045  dirichletRows.push_back(i);
1046  }
1047  }
1048  }
1049 
1050  // Applies Ones-and-Zeros to matrix rows
1051  // Takes a vector of row indices
1052  static void ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1053  const std::vector<LocalOrdinal>& dirichletRows) {
1054  RCP<const Map> Rmap = A->getRowMap();
1055  RCP<const Map> Cmap = A->getColMap();
1056  Scalar one =Teuchos::ScalarTraits<Scalar>::one();
1057  Scalar zero =Teuchos::ScalarTraits<Scalar>::zero();
1058 
1059  for(size_t i=0; i<dirichletRows.size(); i++) {
1060  GlobalOrdinal row_gid = Rmap->getGlobalElement(dirichletRows[i]);
1061 
1062  Teuchos::ArrayView<const LocalOrdinal> indices;
1063  Teuchos::ArrayView<const Scalar> values;
1064  A->getLocalRowView(dirichletRows[i],indices,values);
1065  // NOTE: This won't work with fancy node types.
1066  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1067  for(size_t j=0; j<(size_t)indices.size(); j++) {
1068  if(Cmap->getGlobalElement(indices[j])==row_gid)
1069  valuesNC[j]=one;
1070  else
1071  valuesNC[j]=zero;
1072  }
1073  }
1074  }
1075 
1076  // Applies Ones-and-Zeros to matrix rows
1077  // Takes a Boolean array.
1078  static void ApplyOAZToMatrixRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1079  const Teuchos::ArrayRCP<const bool>& dirichletRows) {
1080  TEUCHOS_ASSERT(A->isFillComplete());
1081  RCP<const Map> domMap = A->getDomainMap();
1082  RCP<const Map> ranMap = A->getRangeMap();
1083  RCP<const Map> Rmap = A->getRowMap();
1084  RCP<const Map> Cmap = A->getColMap();
1085  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == Rmap->getNodeNumElements());
1086  const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1087  const Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
1088  A->resumeFill();
1089  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1090  if (dirichletRows[i]){
1091  GlobalOrdinal row_gid = Rmap->getGlobalElement(i);
1092 
1093  Teuchos::ArrayView<const LocalOrdinal> indices;
1094  Teuchos::ArrayView<const Scalar> values;
1095  A->getLocalRowView(i,indices,values);
1096 
1097  Teuchos::ArrayRCP<Scalar> valuesNC(values.size());
1098  for(size_t j=0; j<(size_t)indices.size(); j++) {
1099  if(Cmap->getGlobalElement(indices[j])==row_gid)
1100  valuesNC[j]=one;
1101  else
1102  valuesNC[j]=zero;
1103  }
1104  A->replaceLocalValues(i,indices,valuesNC());
1105  }
1106  }
1107  A->fillComplete(domMap, ranMap);
1108  }
1109 
1110  // Zeros out rows
1111  // Takes a vector containg Dirichlet row indices
1112  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1113  const std::vector<LocalOrdinal>& dirichletRows,
1114  Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1115  for(size_t i=0; i<dirichletRows.size(); i++) {
1116  Teuchos::ArrayView<const LocalOrdinal> indices;
1117  Teuchos::ArrayView<const Scalar> values;
1118  A->getLocalRowView(dirichletRows[i],indices,values);
1119  // NOTE: This won't work with fancy node types.
1120  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1121  for(size_t j=0; j<(size_t)indices.size(); j++)
1122  valuesNC[j]=replaceWith;
1123  }
1124  }
1125 
1126  // Zeros out rows
1127  // Takes a Boolean ArrayRCP
1128  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& A,
1129  const Teuchos::ArrayRCP<const bool>& dirichletRows,
1130  Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1131  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == A->getRowMap()->getNodeNumElements());
1132  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1133  if (dirichletRows[i]) {
1134  Teuchos::ArrayView<const LocalOrdinal> indices;
1135  Teuchos::ArrayView<const Scalar> values;
1136  A->getLocalRowView(i,indices,values);
1137  // NOTE: This won't work with fancy node types.
1138  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1139  for(size_t j=0; j<(size_t)indices.size(); j++)
1140  valuesNC[j]=replaceWith;
1141  }
1142  }
1143  }
1144 
1145  // Zeros out rows
1146  // Takes a Boolean ArrayRCP
1147  static void ZeroDirichletRows(Teuchos::RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >& X,
1148  const Teuchos::ArrayRCP<const bool>& dirichletRows,
1149  Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1150  TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == X->getMap()->getNodeNumElements());
1151  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1152  if (dirichletRows[i]) {
1153  for(size_t j=0; j<X->getNumVectors(); j++)
1154  X->replaceLocalValue(i,j,replaceWith);
1155  }
1156  }
1157  }
1158 
1159  // Zeros out columns
1160  // Takes a Boolean vector
1161  static void ZeroDirichletCols(Teuchos::RCP<Matrix>& A,
1162  const Teuchos::ArrayRCP<const bool>& dirichletCols,
1163  Scalar replaceWith=Teuchos::ScalarTraits<Scalar>::zero()) {
1164  TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == A->getColMap()->getNodeNumElements());
1165  for(size_t i=0; i<A->getNodeNumRows(); i++) {
1166  Teuchos::ArrayView<const LocalOrdinal> indices;
1167  Teuchos::ArrayView<const Scalar> values;
1168  A->getLocalRowView(i,indices,values);
1169  // NOTE: This won't work with fancy node types.
1170  Scalar* valuesNC = const_cast<Scalar*>(values.getRawPtr());
1171  for(size_t j=0; j<static_cast<size_t>(indices.size()); j++)
1172  if (dirichletCols[indices[j]])
1173  valuesNC[j] = replaceWith;
1174  }
1175  }
1176 
1177  // Finds the OAZ Dirichlet rows for this matrix
1178  static void FindDirichletRowsAndPropagateToCols(Teuchos::RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &A,
1179  Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletRow,
1180  Teuchos::RCP<Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> >& isDirichletCol) {
1181 
1182  // Make sure A's RowMap == DomainMap
1183  if(!A->getRowMap()->isSameAs(*A->getDomainMap())) {
1184  throw std::runtime_error("UtilitiesBase::FindDirichletRowsAndPropagateToCols row and domain maps must match.");
1185  }
1186  RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A->getCrsGraph()->getImporter();
1187  bool has_import = !importer.is_null();
1188 
1189  // Find the Dirichlet rows
1190  std::vector<LocalOrdinal> dirichletRows;
1191  FindDirichletRows(A,dirichletRows);
1192 
1193 #if 0
1194  printf("[%d] DirichletRow Ids = ",A->getRowMap()->getComm()->getRank());
1195  for(size_t i=0; i<(size_t) dirichletRows.size(); i++)
1196  printf("%d ",dirichletRows[i]);
1197  printf("\n");
1198  fflush(stdout);
1199 #endif
1200  // Allocate all as non-Dirichlet
1201  isDirichletRow = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getRowMap(),true);
1202  isDirichletCol = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(A->getColMap(),true);
1203 
1204  Teuchos::ArrayRCP<int> dr_rcp = isDirichletRow->getDataNonConst(0);
1205  Teuchos::ArrayView<int> dr = dr_rcp();
1206  Teuchos::ArrayRCP<int> dc_rcp = isDirichletCol->getDataNonConst(0);
1207  Teuchos::ArrayView<int> dc = dc_rcp();
1208  for(size_t i=0; i<(size_t) dirichletRows.size(); i++) {
1209  dr[dirichletRows[i]] = 1;
1210  if(!has_import) dc[dirichletRows[i]] = 1;
1211  }
1212 
1213  if(has_import)
1214  isDirichletCol->doImport(*isDirichletRow,*importer,Xpetra::CombineMode::ADD);
1215 
1216  }
1217 
1218  // This routine takes a BlockedMap and an Importer (assuming that the BlockedMap matches the source of the importer) and generates a BlockedMap corresponding
1219  // 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)
1220  // This is largely intended to be used in repartitioning of blocked matrices
1221  static RCP<const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> > GeneratedBlockedTargetMap(const Xpetra::BlockedMap<LocalOrdinal,GlobalOrdinal,Node> & sourceBlockedMap,
1222  const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> & Importer) {
1223  typedef Xpetra::Vector<int,LocalOrdinal,GlobalOrdinal,Node> IntVector;
1224  Xpetra::UnderlyingLib lib = sourceBlockedMap.lib();
1225 
1226  // De-stride the map if we have to (might regret this later)
1227  RCP<const Map> fullMap = sourceBlockedMap.getMap();
1228  RCP<const Map> stridedMap = Teuchos::rcp_dynamic_cast<const Xpetra::StridedMap<LocalOrdinal,GlobalOrdinal,Node> >(fullMap);
1229  if(!stridedMap.is_null()) fullMap = stridedMap->getMap();
1230 
1231  // Initial sanity checking for map compatibil
1232  const size_t numSubMaps = sourceBlockedMap.getNumMaps();
1233  if(!Importer.getSourceMap()->isCompatible(*fullMap))
1234  throw std::runtime_error("GenerateBlockedTargetMap(): Map compatibility error");
1235 
1236  // Build an indicator vector
1237  RCP<IntVector> block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(fullMap);
1238 
1239  for(size_t i=0; i<numSubMaps; i++) {
1240  RCP<const Map> map = sourceBlockedMap.getMap(i);
1241 
1242  for(size_t j=0; j<map->getNodeNumElements(); j++) {
1243  LocalOrdinal jj = fullMap->getLocalElement(map->getGlobalElement(j));
1244  block_ids->replaceLocalValue(jj,(int)i);
1245  }
1246  }
1247 
1248  // Get the block ids for the new map
1249  RCP<const Map> targetMap = Importer.getTargetMap();
1250  RCP<IntVector> new_block_ids = Xpetra::VectorFactory<int,LocalOrdinal,GlobalOrdinal,Node>::Build(targetMap);
1251  new_block_ids->doImport(*block_ids,Importer,Xpetra::CombineMode::ADD);
1252  Teuchos::ArrayRCP<const int> dataRCP = new_block_ids->getData(0);
1253  Teuchos::ArrayView<const int> data = dataRCP();
1254 
1255 
1256  // Get the GIDs for each subblock
1257  Teuchos::Array<Teuchos::Array<GlobalOrdinal> > elementsInSubMap(numSubMaps);
1258  for(size_t i=0; i<targetMap->getNodeNumElements(); i++) {
1259  elementsInSubMap[data[i]].push_back(targetMap->getGlobalElement(i));
1260  }
1261 
1262  // Generate the new submaps
1263  std::vector<RCP<const Map> > subMaps(numSubMaps);
1264  for(size_t i=0; i<numSubMaps; i++) {
1265  subMaps[i] = Xpetra::MapFactory<LocalOrdinal,GlobalOrdinal,Node>::Build(lib,Teuchos::OrdinalTraits<GlobalOrdinal>::invalid(),elementsInSubMap[i](),targetMap->getIndexBase(),targetMap->getComm());
1266  }
1267 
1268  // Build the BlockedMap
1269  return rcp(new BlockedMap(targetMap,subMaps));
1270 
1271  }
1272 
1273  // Checks to see if the first chunk of the colMap is also the row map. This simiplifies a bunch of
1274  // operation in coarsening
1275  static bool MapsAreNested(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& rowMap, const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node>& colMap) {
1276  ArrayView<const GlobalOrdinal> rowElements = rowMap.getNodeElementList();
1277  ArrayView<const GlobalOrdinal> colElements = colMap.getNodeElementList();
1278 
1279  const size_t numElements = rowElements.size();
1280 
1281  if (size_t(colElements.size()) < numElements)
1282  return false;
1283 
1284  bool goodMap = true;
1285  for (size_t i = 0; i < numElements; i++)
1286  if (rowElements[i] != colElements[i]) {
1287  goodMap = false;
1288  break;
1289  }
1290 
1291  return goodMap;
1292  }
1293 
1294 
1295 
1296 
1297  }; // class Utils
1298 
1299 
1301 
1302 } //namespace MueLu
1303 
1304 #define MUELU_UTILITIESBASE_SHORT
1305 #endif // MUELU_UTILITIESBASE_DECL_HPP
Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
Xpetra::BlockedMultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedMultiVector
MueLu::DefaultLocalOrdinal LocalOrdinal
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.
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
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
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Apply Rowsum Criterion.
static RCP< const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > > GeneratedBlockedTargetMap(const Xpetra::BlockedMap< LocalOrdinal, GlobalOrdinal, Node > &sourceBlockedMap, const Xpetra::Import< LocalOrdinal, GlobalOrdinal, Node > &Importer)
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &BlockNumber)
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows)
static void FindDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, std::vector< LocalOrdinal > &dirichletRows, bool count_twos_as_dirichlet=false)
magnitude_type tolerance
static RCP< Vector > GetMatrixOverlappedDeletedRowsum(const Matrix &A)
Extract Overlapped Matrix Deleted Rowsum.
Namespace for MueLu classes and methods.
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)
MueLu::DefaultNode Node
#define MueLu_sumAll(rcpComm, in, out)
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.
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i=k}(-a_ik), for each for i in the matrix.
static bool MapsAreNested(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &rowMap, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &colMap)
static RCP< Xpetra::Vector< Magnitude, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedAbsDeletedRowsum(const Matrix &A)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
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)
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::ArrayView< double > &weight, const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Weighted squared distance between two rows in a multivector.
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false)
Extract Matrix Diagonal of lumped matrix.
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
Xpetra::BlockedVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > BlockedVector
Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrix
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 void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Extract Matrix Diagonal.
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > &BlockNumber, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
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)
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 void FindNonZeros(const Teuchos::ArrayRCP< const Scalar > vals, Teuchos::ArrayRCP< bool > nonzeros)
Find non-zero values in an ArrayRCP Compares the value to 2 * machine epsilon.
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static Scalar PowerMethod(const Matrix &A, const RCP< Vector > &diagInvVec, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Power method.
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 void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< bool > &dirichletRows, Teuchos::ArrayRCP< bool > dirichletCols, Teuchos::ArrayRCP< bool > dirichletDomain)
Detects Dirichlet columns & domains from a list of Dirichlet rows.
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.