MueLu  Version of the Day
MueLu_Utilities_kokkos_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_UTILITIES_KOKKOS_DECL_HPP
47 #define MUELU_UTILITIES_KOKKOS_DECL_HPP
48 
49 #include "MueLu_ConfigDefs.hpp"
50 #if defined(HAVE_MUELU_KOKKOS_REFACTOR)
51 
52 #include <string>
53 
54 #include "Teuchos_ScalarTraits.hpp"
55 #include "Teuchos_ParameterList.hpp"
56 
57 #include "Xpetra_BlockedCrsMatrix_fwd.hpp"
58 #include "Xpetra_CrsMatrix_fwd.hpp"
59 #include "Xpetra_CrsMatrixWrap_fwd.hpp"
60 #include "Xpetra_ExportFactory.hpp"
61 #include "Xpetra_ImportFactory_fwd.hpp"
62 #include "Xpetra_MapFactory_fwd.hpp"
63 #include "Xpetra_Map_fwd.hpp"
64 #include "Xpetra_MatrixFactory_fwd.hpp"
65 #include "Xpetra_Matrix_fwd.hpp"
66 #include "Xpetra_MultiVectorFactory_fwd.hpp"
67 #include "Xpetra_MultiVector_fwd.hpp"
68 #include "Xpetra_Operator_fwd.hpp"
69 #include "Xpetra_VectorFactory_fwd.hpp"
70 #include "Xpetra_Vector_fwd.hpp"
71 
72 #include "Xpetra_IO.hpp"
73 
74 #include "Kokkos_ArithTraits.hpp"
75 
76 #ifdef HAVE_MUELU_EPETRA
77 #include "Epetra_MultiVector.h"
78 #include "Epetra_CrsMatrix.h"
79 #include "Xpetra_EpetraCrsMatrix_fwd.hpp"
80 #include "Xpetra_EpetraMultiVector_fwd.hpp"
81 #endif
82 
83 #include "MueLu_Exceptions.hpp"
84 #include "MueLu_Utilities.hpp"
85 #include "MueLu_UtilitiesBase.hpp"
86 
87 #ifdef HAVE_MUELU_TPETRA
88 #include "Tpetra_CrsMatrix.hpp"
89 #include "Tpetra_Map.hpp"
90 #include "Tpetra_MultiVector.hpp"
91 #include "Xpetra_TpetraCrsMatrix_fwd.hpp"
92 #include "Xpetra_TpetraMultiVector_fwd.hpp"
93 #endif
94 
95 
96 namespace MueLu {
97 
105  template <class Scalar,
108  class Node = DefaultNode>
109  class Utilities_kokkos : public MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
110 #undef MUELU_UTILITIES_KOKKOS_SHORT
111 #include "MueLu_UseShortNames.hpp"
112 
113  public:
114  using TST = Teuchos::ScalarTraits<SC>;
115  using Magnitude = typename TST::magnitudeType;
116  using RealValuedMultiVector = Xpetra::MultiVector<Magnitude,LO,GO,NO>;
117 
118 #ifdef HAVE_MUELU_EPETRA
119  // @{
121  static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2EpetraMV(vec); }
122  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstEpetraMV(vec); }
123 
124  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec) { return Utilities::MV2EpetraMV(vec); }
125  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec) { return Utilities::MV2NonConstEpetraMV(vec); }
126 
127  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2EpetraCrs(Op); }
128  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
129 
130  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) { return Utilities::Op2EpetraCrs(Op); }
131  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
132 
133  static const Epetra_Map& Map2EpetraMap(const Map& map) { return Utilities::Map2EpetraMap(map); }
134  // @}
135 #endif
136 
137 #ifdef HAVE_MUELU_TPETRA
138  // @{
140  static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2TpetraMV(vec); }
141  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstTpetraMV(vec); }
142  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV2(vec); }
143 
144  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec) { return Utilities::MV2TpetraMV(vec); }
145  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV(vec); }
146 
147  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2TpetraCrs(Op); }
148  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
149 
150  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op) { return Utilities::Op2TpetraCrs(Op); }
151  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
152 
153  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op) { return Utilities::Op2TpetraRow(Op); }
154  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraRow(Op); }
155 
156  static const RCP<const Tpetra::Map<LO, GO, NO> > Map2TpetraMap(const Map& map) { return Utilities::Map2TpetraMap(map); }
157 #endif
158 
159  static RCP<Xpetra::Matrix<SC,LO,GO,NO> > Crs2Op(RCP<CrsMatrix> Op) { return Utilities::Crs2Op(Op); }
160 
167  static RCP<Vector> GetMatrixDiagonal(const Matrix& A); // FIXME
168 
175  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = TST::eps()*100, const bool doLumped = false); // FIXME
176 
177 
178 
186  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A); // FIXME
187 
195  static Teuchos::RCP<Vector> GetInverse(Teuchos::RCP<const Vector> v, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar tolReplacement = Teuchos::ScalarTraits<Scalar>::zero()) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetInverse(v,tol,tolReplacement); }
196 
197  // TODO: should NOT return an Array. Definition must be changed to:
198  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
199  // or
200  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
201  static Teuchos::Array<Magnitude> ResidualNorm(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
202  return Utilities::ResidualNorm(Op, X, RHS);
203  }
204 
205  static RCP<MultiVector> Residual(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
206  return Utilities::Residual(Op, X, RHS);
207  }
208 
224  static SC PowerMethod(const Matrix& A, bool scaleByDiag = true,
225  LO niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
226  return Utilities::PowerMethod(A, scaleByDiag, niters, tolerance, verbose, seed);
227  }
228 
229  static SC PowerMethod(const Matrix& A, const Teuchos::RCP<Vector> &invDiag,
230  LO niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
231  return Utilities::PowerMethod(A, invDiag, niters, tolerance, verbose, seed);
232  }
233 
234  static void MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector, bool doInverse = true,
235  bool doFillComplete = true, bool doOptimizeStorage = true); // FIXME
236 
237  static void MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
238  bool doFillComplete, bool doOptimizeStorage); // FIXME
239 
240  static void MyOldScaleMatrix_Epetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
241  bool doFillComplete, bool doOptimizeStorage); // FIXME
242 
243  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) { return Utilities::MakeFancy(os); }
244 
252  static Kokkos::View<bool*, typename NO::device_type> DetectDirichletRows(const Matrix& A, const Magnitude& tol = Teuchos::ScalarTraits<SC>::zero(), const bool count_twos_as_dirichlet=false);
253 
254 
261  static void FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
263 
264 
276 
277 
285  static void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
289 
290 
291  static void ZeroDirichletRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows, SC replaceWith=Teuchos::ScalarTraits<SC>::zero());
292 
293  static void ZeroDirichletRows(RCP<MultiVector>& X, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows, SC replaceWith=Teuchos::ScalarTraits<SC>::zero());
294 
295  static void ZeroDirichletCols(RCP<Matrix>& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletCols, SC replaceWith=Teuchos::ScalarTraits<SC>::zero());
296 
297  static void ApplyRowSumCriterion(const Matrix& A,
298  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
300 
301  static RCP<MultiVector> RealValuedToScalarMultiVector(RCP<RealValuedMultiVector> X);
302 
312  static void SetRandomSeed(const Teuchos::Comm<int> &comm) { return Utilities::SetRandomSeed(comm); }
313 
319  static RCP<Matrix> Transpose(Matrix& Op, bool optimizeTranspose = false, const std::string & label = std::string()) {
320  return Utilities::Transpose(Op, optimizeTranspose, label);
321  }
322 
323  static RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LocalOrdinal,GlobalOrdinal,Node> > ExtractCoordinatesFromParameterList(ParameterList& paramList) {
325  }
326 
330  static RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > CuthillMcKee(const Matrix &Op);
331 
334  static RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > ReverseCuthillMcKee(const Matrix &Op);
335 
336  static void ApplyOAZToMatrixRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows);
337 
338  }; // class Utils
339 
340 
350  template <class Node>
351  class Utilities_kokkos<double,int,int,Node> : public UtilitiesBase<double,int,int,Node> {
352  public:
353  typedef double Scalar;
354  typedef int LocalOrdinal;
355  typedef int GlobalOrdinal;
356  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
357  typedef Xpetra::MultiVector<Magnitude,LocalOrdinal,GlobalOrdinal,Node> RealValuedMultiVector;
358 
359  private:
360 #undef MUELU_UTILITIES_KOKKOS_SHORT
361 #include "MueLu_UseShortNames.hpp"
362 
363  public:
364 
365 #ifdef HAVE_MUELU_EPETRA
366  // @{
368  static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2EpetraMV(vec); }
369  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstEpetraMV(vec); }
370 
371  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec) { return Utilities::MV2EpetraMV(vec); }
372  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec) { return Utilities::MV2NonConstEpetraMV(vec); }
373 
374  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2EpetraCrs(Op); }
375  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
376 
377  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) { return Utilities::Op2EpetraCrs(Op); }
378  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
379 
380  static const Epetra_Map& Map2EpetraMap(const Map& map) { return Utilities::Map2EpetraMap(map); }
381  // @}
382 #endif
383 
384 #ifdef HAVE_MUELU_TPETRA
385  // @{
387  static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2TpetraMV(vec); }
388  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstTpetraMV(vec); }
389  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV2(vec); }
390 
391  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec) { return Utilities::MV2TpetraMV(vec); }
392  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV(vec); }
393 
394  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2TpetraCrs(Op); }
395  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
396 
397  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op) { return Utilities::Op2TpetraCrs(Op); }
398  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
399 
400  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op) { return Utilities::Op2TpetraRow(Op); }
401  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraRow(Op); }
402 
403  static const RCP<const Tpetra::Map<LO, GO, NO> > Map2TpetraMap(const Map& map) { return Utilities::Map2TpetraMap(map); }
404 #endif
405  static RCP<Xpetra::Matrix<SC,LO,GO,NO> > Crs2Op(RCP<CrsMatrix> Op) { return Utilities::Crs2Op(Op); }
406 
407  static RCP<Vector> GetMatrixDiagonal(const Matrix& A) {
408  const auto rowMap = A.getRowMap();
409  auto diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
410 
411  A.getLocalDiagCopy(*diag);
412 
413  return diag;
414  }
415  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<SC>::eps()*100, const bool doLumped=false) {
416  return UtilitiesBase::GetMatrixDiagonalInverse(A, tol, doLumped);
417  }
418  static RCP<Vector> GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100, Scalar tolReplacement = Teuchos::ScalarTraits<Scalar>::zero(), const bool replaceSingleEntryRowWithZero = false) {
419  return UtilitiesBase::GetLumpedMatrixDiagonal(A, doReciprocal, tol, tolReplacement, replaceSingleEntryRowWithZero);
420  }
421  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A) {
423  }
424  static RCP<Vector> GetInverse(RCP<const Vector> v, Magnitude tol = Teuchos::ScalarTraits<SC>::eps()*100, SC tolReplacement = Teuchos::ScalarTraits<SC>::zero()) {
425  return UtilitiesBase::GetInverse(v,tol,tolReplacement);
426  }
427  static Teuchos::Array<Magnitude> ResidualNorm(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
428  return UtilitiesBase::ResidualNorm(Op,X,RHS);
429  }
430  static RCP<MultiVector> Residual(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
431  return UtilitiesBase::Residual(Op,X,RHS);
432  }
433  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) {
434  return UtilitiesBase::MakeFancy(os);
435  }
436  static void SetRandomSeed(const Teuchos::Comm<int> &comm) {
438  }
439 
440  // todo: move this to UtilitiesBase::kokkos
441  static Kokkos::View<bool*, typename Node::device_type> DetectDirichletRows(const Matrix& A, const Magnitude& tol = Teuchos::ScalarTraits<SC>::zero(), const bool count_twos_as_dirichlet=false);
442 
444 
445  static void FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
447 
448  static void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
452 
453  static void ZeroDirichletRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows, SC replaceWith=Teuchos::ScalarTraits<SC>::zero());
454 
455  static void ZeroDirichletRows(RCP<MultiVector>& X, const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows, SC replaceWith=Teuchos::ScalarTraits<SC>::zero());
456 
457  static void ZeroDirichletCols(RCP<Matrix>& A, const Kokkos::View<const bool*, typename Node::device_type>& dirichletCols, SC replaceWith=Teuchos::ScalarTraits<SC>::zero());
458 
459  static void ApplyRowSumCriterion(const Matrix& A,
460  const typename Teuchos::ScalarTraits<Scalar>::magnitudeType rowSumTol,
462 
463  static RCP<MultiVector> RealValuedToScalarMultiVector(RCP<RealValuedMultiVector> X);
464 
465  static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true, LO niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
466  return UtilitiesBase::PowerMethod(A,scaleByDiag,niters,tolerance,verbose,seed);
467  }
468 
469  static Scalar PowerMethod(const Matrix& A, const Teuchos::RCP<Vector> &invDiag, LO niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
470  return UtilitiesBase::PowerMethod(A, invDiag, niters, tolerance, verbose, seed);
471  }
472 
473  static void MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector, bool doInverse = true,
474  bool doFillComplete = true, bool doOptimizeStorage = true) {
475  SC one = Teuchos::ScalarTraits<SC>::one();
476  Teuchos::ArrayRCP<SC> sv(scalingVector.size());
477  if (doInverse) {
478  for (int i = 0; i < scalingVector.size(); ++i)
479  sv[i] = one / scalingVector[i];
480  } else {
481  for (int i = 0; i < scalingVector.size(); ++i)
482  sv[i] = scalingVector[i];
483  }
484 
485  switch (Op.getRowMap()->lib()) {
486  case Xpetra::UseTpetra:
487  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
488  break;
489 
490  case Xpetra::UseEpetra:
491  MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
492  break;
493 
494  default:
495  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
496 #ifndef __NVCC__ //prevent nvcc warning
497  break;
498 #endif
499  }
500  }
501 
502  // TODO This is the <double,int,int> specialization
503  static void MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector,
504  bool doFillComplete, bool doOptimizeStorage) {
505  #ifdef HAVE_MUELU_TPETRA
506  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
507  try {
508  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpOp = Op2NonConstTpetraCrs(Op);
509 
510  const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
511  const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
512  const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
513 
514  size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
515  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
516  maxRowSize = 20;
517 
518  std::vector<SC> scaledVals(maxRowSize);
519  if (tpOp.isFillComplete())
520  tpOp.resumeFill();
521 
522  if (Op.isLocallyIndexed() == true) {
523  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::local_inds_host_view_type cols;
524  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::values_host_view_type vals;
525 
526  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
527  tpOp.getLocalRowView(i, cols, vals);
528  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
529  if (nnz > maxRowSize) {
530  maxRowSize = nnz;
531  scaledVals.resize(maxRowSize);
532  }
533  for (size_t j = 0; j < nnz; ++j)
534  scaledVals[j] = vals[j]*scalingVector[i];
535 
536  if (nnz > 0) {
537  Teuchos::ArrayView<const LocalOrdinal> cols_view(cols.data(), nnz);
538  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
539  tpOp.replaceLocalValues(i, cols_view, valview);
540  }
541  } //for (size_t i=0; ...
542 
543  } else {
544  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::global_inds_host_view_type cols;
545  typename Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::values_host_view_type vals;
546 
547  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
548  GO gid = rowMap->getGlobalElement(i);
549  tpOp.getGlobalRowView(gid, cols, vals);
550  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
551  if (nnz > maxRowSize) {
552  maxRowSize = nnz;
553  scaledVals.resize(maxRowSize);
554  }
555  // FIXME FIXME FIXME FIXME FIXME FIXME
556  for (size_t j = 0; j < nnz; ++j)
557  scaledVals[j] = vals[j]*scalingVector[i]; //FIXME i or gid?
558 
559  if (nnz > 0) {
560  Teuchos::ArrayView<const GlobalOrdinal> cols_view(cols.data(), nnz);
561  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
562  tpOp.replaceGlobalValues(gid, cols_view, valview);
563  }
564  } //for (size_t i=0; ...
565  }
566 
567  if (doFillComplete) {
568  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
569  throw Exceptions::RuntimeError("In Utilities::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
570 
571  RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
572  params->set("Optimize Storage", doOptimizeStorage);
573  params->set("No Nonlocal Changes", true);
574  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
575  }
576  } catch(...) {
577  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
578  }
579 #else
580  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been compiled with support for LO=GO=int.");
581 #endif
582 #else
583  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
584 #endif
585  }
586 
587  static void MyOldScaleMatrix_Epetra(Matrix& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector, bool doFillComplete, bool doOptimizeStorage) {
588 #ifdef HAVE_MUELU_EPETRA
589  try {
590  //const Epetra_CrsMatrix& epOp = Utilities<double,int,int>::Op2NonConstEpetraCrs(Op);
591  const Epetra_CrsMatrix& epOp = Op2NonConstEpetraCrs(Op);
592 
593  Epetra_Map const &rowMap = epOp.RowMap();
594  int nnz;
595  double *vals;
596  int *cols;
597 
598  for (int i = 0; i < rowMap.NumMyElements(); ++i) {
599  epOp.ExtractMyRowView(i, nnz, vals, cols);
600  for (int j = 0; j < nnz; ++j)
601  vals[j] *= scalingVector[i];
602  }
603 
604  } catch (...){
605  throw Exceptions::RuntimeError("Only Epetra_CrsMatrix types can be scaled");
606  }
607 #else
608  throw Exceptions::RuntimeError("Matrix scaling is not possible because Epetra has not been enabled.");
609 #endif // HAVE_MUELU_EPETRA
610  }
611 
617  static RCP<Matrix> Transpose(Matrix& Op, bool /* optimizeTranspose */ = false, const std::string & label = std::string(),const Teuchos::RCP<Teuchos::ParameterList> &params=Teuchos::null) {
618  switch (Op.getRowMap()->lib()) {
619  case Xpetra::UseTpetra:
620  {
621 #ifdef HAVE_MUELU_TPETRA
622 #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
623  try {
624  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpetraOp = Utilities::Op2TpetraCrs(Op);
625 
626  // Compute the transpose A of the Tpetra matrix tpetraOp.
627  RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
628  Tpetra::RowMatrixTransposer<SC,LO,GO,NO> transposer(rcpFromRef(tpetraOp),label);
629  {
630  using Teuchos::ParameterList;
631  using Teuchos::rcp;
632  RCP<ParameterList> transposeParams = params.is_null () ?
633  rcp (new ParameterList) :
634  rcp (new ParameterList (*params));
635  transposeParams->set ("sort", false);
636  A = transposer.createTranspose(transposeParams);
637  }
638 
639  RCP<Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > AA = rcp(new Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>(A));
640  RCP<CrsMatrix> AAA = rcp_implicit_cast<CrsMatrix>(AA);
641  RCP<CrsMatrixWrap> AAAA = rcp(new CrsMatrixWrap(AAA));
642 
643  return AAAA;
644  }
645  catch (std::exception& e) {
646  std::cout << "threw exception '" << e.what() << "'" << std::endl;
647  throw Exceptions::RuntimeError("Utilities::Transpose failed, perhaps because matrix is not a Crs matrix");
648  }
649 #else
650  throw Exceptions::RuntimeError("Utilities::Transpose: Tpetra is not compiled with LO=GO=int. Add TPETRA_INST_INT_INT:BOOL=ON to your configuration!");
651 #endif
652 #else
653  throw Exceptions::RuntimeError("Utilities::Transpose: Tpetra is not compiled!");
654 #endif
655 #ifndef __NVCC__ //prevent nvcc warning
656  break;
657 #endif
658  }
659  case Xpetra::UseEpetra:
660  {
661 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
662  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("ZZ Entire Transpose"));
663  // Epetra case
666  Epetra_CrsMatrix * A = dynamic_cast<Epetra_CrsMatrix*>(&transposer(epetraOp));
667  transposer.ReleaseTranspose(); // So we can keep A in Muelu...
668 
669  RCP<Epetra_CrsMatrix> rcpA(A);
670  RCP<Xpetra::EpetraCrsMatrixT<GO,NO> > AA = rcp(new Xpetra::EpetraCrsMatrixT<GO,NO> (rcpA));
671  RCP<CrsMatrix> AAA = rcp_implicit_cast<CrsMatrix>(AA);
672  RCP<CrsMatrixWrap> AAAA = rcp( new CrsMatrixWrap(AAA));
673  AAAA->fillComplete(Op.getRangeMap(), Op.getDomainMap());
674 
675  return AAAA;
676 #else
677  throw Exceptions::RuntimeError("Epetra (Err. 2)");
678 #endif
679 #ifndef __NVCC__ //prevent nvcc warning
680  break;
681 #endif
682  }
683  default:
684  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be transposed.");
685 #ifndef __NVCC__ //prevent nvcc warning
686  break;
687 #endif
688  }
689 
690 #ifndef __NVCC__ //prevent nvcc warning
691  return Teuchos::null;
692 #endif
693  }
694 
697  static RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > ExtractCoordinatesFromParameterList(ParameterList& paramList) {
698  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > coordinates = Teuchos::null;
699 
700  // check whether coordinates are contained in parameter list
701  if(paramList.isParameter ("Coordinates") == false)
702  return coordinates;
703 
704 #if defined(HAVE_MUELU_TPETRA)
705 #if ( defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT)) || \
706  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))
707 
708  // define Tpetra::MultiVector type with Scalar=float only if
709  // * ETI is turned off, since then the compiler will instantiate it automatically OR
710  // * Tpetra is instantiated on Scalar=float
711 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
712  typedef Tpetra::MultiVector<float, LO,GO,NO> tfMV;
713  RCP<tfMV> floatCoords = Teuchos::null;
714 #endif
715 
716  // define Tpetra::MultiVector type with Scalar=double only if
717  // * ETI is turned off, since then the compiler will instantiate it automatically OR
718  // * Tpetra is instantiated on Scalar=double
719  typedef Tpetra::MultiVector<SC,LO,GO,NO> tdMV;
720  RCP<tdMV> doubleCoords = Teuchos::null;
721  if (paramList.isType<RCP<tdMV> >("Coordinates")) {
722  // Coordinates are stored as a double vector
723  doubleCoords = paramList.get<RCP<tdMV> >("Coordinates");
724  paramList.remove("Coordinates");
725  }
726 #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
727  else if (paramList.isType<RCP<tfMV> >("Coordinates")) {
728  // check if coordinates are stored as a float vector
729  floatCoords = paramList.get<RCP<tfMV> >("Coordinates");
730  paramList.remove("Coordinates");
731  doubleCoords = rcp(new tdMV(floatCoords->getMap(), floatCoords->getNumVectors()));
732  deep_copy(*doubleCoords, *floatCoords);
733  }
734 #endif
735  // We have the coordinates in a Tpetra double vector
736  if(doubleCoords != Teuchos::null) {
737  coordinates = Teuchos::rcp(new Xpetra::TpetraMultiVector<SC,LO,GO,NO>(doubleCoords));
738  TEUCHOS_TEST_FOR_EXCEPT(doubleCoords->getNumVectors() != coordinates->getNumVectors());
739  }
740 #endif // Tpetra instantiated on GO=int and EpetraNode
741 #endif // endif HAVE_TPETRA
742 
743 #if defined(HAVE_MUELU_EPETRA)
744  RCP<Epetra_MultiVector> doubleEpCoords;
745  if (paramList.isType<RCP<Epetra_MultiVector> >("Coordinates")) {
746  doubleEpCoords = paramList.get<RCP<Epetra_MultiVector> >("Coordinates");
747  paramList.remove("Coordinates");
748  RCP<Xpetra::EpetraMultiVectorT<GO,NO> > epCoordinates = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GO,NO>(doubleEpCoords));
749  coordinates = rcp_dynamic_cast<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> >(epCoordinates);
750  TEUCHOS_TEST_FOR_EXCEPT(doubleEpCoords->NumVectors() != Teuchos::as<int>(coordinates->getNumVectors()));
751  }
752 #endif
753  TEUCHOS_TEST_FOR_EXCEPT(Teuchos::is_null(coordinates));
754  return coordinates;
755  }
756 
760  static RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > ReverseCuthillMcKee(const Matrix &Op);
761 
765  static RCP<Xpetra::Vector<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node> > CuthillMcKee(const Matrix &Op);
766 
767  static void ApplyOAZToMatrixRows(RCP<Matrix>& A, const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows);
768 
769  }; // class Utilities (specialization SC=double LO=GO=int)
770 
771 
772 
773  // Useful Kokkos conversions
774  template < class View, unsigned AppendValue >
775  struct AppendTrait {
776  // static_assert(false, "Error: NOT a Kokkos::View");
777  };
778 
779  template < class MT, unsigned T >
780  struct CombineMemoryTraits {
781  // empty
782  };
783 
784  template < unsigned U, unsigned T>
785  struct CombineMemoryTraits<Kokkos::MemoryTraits<U>, T> {
786  typedef Kokkos::MemoryTraits<U|T> type;
787  };
788 
789  template < class DataType, unsigned T, class... Pack >
790  struct AppendTrait< Kokkos::View< DataType, Pack... >, T> {
791  typedef Kokkos::View< DataType, Pack... > view_type;
793  };
794 
795 } //namespace MueLu
796 
797 #define MUELU_UTILITIES_KOKKOS_SHORT
798 
799 #endif // #if defined(HAVE_MUELU_KOKKOS_REFACTOR)
800 
801 #endif // MUELU_UTILITIES_KOKKOS_DECL_HPP
const Epetra_Map & RowMap() const
View
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
static const RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > Map2TpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
MueLu::DefaultLocalOrdinal LocalOrdinal
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
static RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
Teuchos::RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
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 RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
static RCP< const Epetra_MultiVector > MV2EpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
void ZeroDirichletRows(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, Scalar replaceWith)
Teuchos::RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op)
Namespace for MueLu classes and methods.
void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows)
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
MueLu::DefaultNode Node
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
void DetectDirichletColsAndDomains(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Kokkos::View< bool *, typename Node::device_type > &dirichletRows, Kokkos::View< bool *, typename Node::device_type > dirichletCols, Kokkos::View< bool *, typename Node::device_type > dirichletDomain)
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
void deep_copy(const DynRankView< DT, DP... > &dst, typename ViewTraits< DT, DP... >::const_value_type &value, typename std::enable_if< std::is_same< typename ViewTraits< DT, DP... >::specialize, void >::value >::type *=nullptr)
int NumMyElements() const
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.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
static RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Extract Matrix Diagonal.
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Crs2Op(RCP< Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
void FindNonZeros(const typename Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >::dual_view_type::t_dev_const_um vals, Kokkos::View< bool *, typename Node::device_type > nonzeros)
void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, typename Node::device_type > &dirichletRows)
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
Kokkos::View< bool *, typename NO::device_type > DetectDirichletCols(const Xpetra::Matrix< SC, LO, GO, NO > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows)
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static Scalar PowerMethod(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.