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 #ifdef HAVE_MUELU_EPETRA
75 #include <Epetra_MultiVector.h>
76 #include <Epetra_CrsMatrix.h>
77 #include <Xpetra_EpetraCrsMatrix_fwd.hpp>
78 #include <Xpetra_EpetraMultiVector_fwd.hpp>
79 #endif
80 
81 #include "MueLu_Exceptions.hpp"
82 #include "MueLu_Utilities.hpp"
83 #include "MueLu_UtilitiesBase.hpp"
84 
85 #ifdef HAVE_MUELU_TPETRA
86 #include <Tpetra_CrsMatrix.hpp>
87 #include <Tpetra_Map.hpp>
88 #include <Tpetra_MultiVector.hpp>
89 #include <Xpetra_TpetraCrsMatrix_fwd.hpp>
90 #include <Xpetra_TpetraMultiVector_fwd.hpp>
91 #endif
92 
93 
94 namespace MueLu {
95 
103  template <class Scalar,
104  class LocalOrdinal = int,
105  class GlobalOrdinal = LocalOrdinal,
106  class Node = KokkosClassic::DefaultNode::DefaultNodeType>
107  class Utilities_kokkos : public MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
108 #undef MUELU_UTILITIES_KOKKOS_SHORT
109 #include "MueLu_UseShortNames.hpp"
110 
111  public:
112  typedef typename Teuchos::ScalarTraits<SC>::magnitudeType Magnitude;
113 
114 #ifdef HAVE_MUELU_EPETRA
115  // @{
117  static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2EpetraMV(vec); }
118  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstEpetraMV(vec); }
119 
120  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec) { return Utilities::MV2EpetraMV(vec); }
121  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec) { return Utilities::MV2NonConstEpetraMV(vec); }
122 
123  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2EpetraCrs(Op); }
124  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
125 
126  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) { return Utilities::Op2EpetraCrs(Op); }
127  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
128 
129  static const Epetra_Map& Map2EpetraMap(const Map& map) { return Utilities::Map2EpetraMap(map); }
130  // @}
131 #endif
132 
133 #ifdef HAVE_MUELU_TPETRA
134  // @{
136  static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2TpetraMV(vec); }
137  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstTpetraMV(vec); }
138  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV2(vec); }
139 
140  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec) { return Utilities::MV2TpetraMV(vec); }
141  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV(vec); }
142 
143  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2TpetraCrs(Op); }
144  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
145 
146  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op) { return Utilities::Op2TpetraCrs(Op); }
147  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
148 
149  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op) { return Utilities::Op2TpetraRow(Op); }
150  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraRow(Op); }
151 
152  static const RCP<const Tpetra::Map<LO, GO, NO> > Map2TpetraMap(const Map& map) { return Utilities::Map2TpetraMap(map); }
153 #endif
154 
155  static RCP<Xpetra::Matrix<SC,LO,GO,NO> > Crs2Op(RCP<CrsMatrix> Op) { return Utilities::Crs2Op(Op); }
156 
163  static Teuchos::ArrayRCP<SC> GetMatrixDiagonal(const Matrix& A); // FIXME
164 
171  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<SC>::eps()*100); // FIXME
172 
173 
174 
181  static Teuchos::ArrayRCP<SC> GetLumpedMatrixDiagonal(const Matrix& A); // FIXME
182 
190  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A); // FIXME
191 
192  // TODO: should NOT return an Array. Definition must be changed to:
193  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
194  // or
195  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
196  static Teuchos::Array<Magnitude> ResidualNorm(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
197  return Utilities::ResidualNorm(Op, X, RHS);
198  }
199 
200  static RCP<MultiVector> Residual(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
201  return Utilities::Residual(Op, X, RHS);
202  }
203 
204  static void PauseForDebugger();
205 
221  static SC PowerMethod(const Matrix& A, bool scaleByDiag = true,
222  LO niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
223  return Utilities::PowerMethod(A, scaleByDiag, niters, tolerance, verbose, seed);
224  }
225 
226  static void MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector, bool doInverse = true,
227  bool doFillComplete = true, bool doOptimizeStorage = true); // FIXME
228 
229  static void MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
230  bool doFillComplete, bool doOptimizeStorage); // FIXME
231 
232  static void MyOldScaleMatrix_Epetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
233  bool doFillComplete, bool doOptimizeStorage); // FIXME
234 
235  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) { return Utilities::MakeFancy(os); }
236 
241  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const MultiVector& v, LocalOrdinal i0, LocalOrdinal i1) { return Utilities::Distance2(v, i0, i1); }
242 
250  static Kokkos::View<const bool*, typename NO::device_type> DetectDirichletRows(const Matrix& A, const Magnitude& tol = Teuchos::ScalarTraits<SC>::zero());
251 
261  static void SetRandomSeed(const Teuchos::Comm<int> &comm) { return Utilities::SetRandomSeed(comm); }
262 
268  static RCP<Matrix> Transpose(Matrix& Op, bool optimizeTranspose = false, const std::string & label = std::string()) {
269  return Utilities::Transpose(Op, optimizeTranspose, label);
270  }
271 
272  }; // class Utils
273 
274 
284  template <class Node>
285  class Utilities_kokkos<double,int,int,Node> : public UtilitiesBase<double,int,int,Node> {
286  public:
287  typedef double Scalar;
288  typedef int LocalOrdinal;
289  typedef int GlobalOrdinal;
290  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
291 
292  #ifdef HAVE_MUELU_EPETRA
293  // @{
295  static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > const vec) {
296  RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(vec);
297  if (tmpVec == Teuchos::null)
298  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
299  return tmpVec->getEpetra_MultiVector();
300  }
301  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > vec) {
302  RCP<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> >(vec);
303  if (tmpVec == Teuchos::null)
304  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::EpetraMultiVector failed");
305  return tmpVec->getEpetra_MultiVector();
306  }
307 
308  static const Epetra_MultiVector& MV2EpetraMV(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & vec) {
309  const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> & tmpVec = dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> &>(vec);
310  return *(tmpVec.getEpetra_MultiVector());
311  }
312  static Epetra_MultiVector& MV2NonConstEpetraMV(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & vec) {
313  const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> & tmpVec = dynamic_cast<const Xpetra::EpetraMultiVectorT<GlobalOrdinal,Node> &>(vec);
314  return *(tmpVec.getEpetra_MultiVector());
315  }
316 
317  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
318  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
319  if (crsOp == Teuchos::null)
320  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
321  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >& tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
322  if (tmp_ECrsMtx == Teuchos::null)
323  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
324  return tmp_ECrsMtx->getEpetra_CrsMatrix();
325  }
326  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
327  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
328  if (crsOp == Teuchos::null)
329  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
330  const RCP<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
331  if (tmp_ECrsMtx == Teuchos::null)
332  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
333  return tmp_ECrsMtx->getEpetra_CrsMatrixNonConst();
334  }
335 
336  static const Epetra_CrsMatrix& Op2EpetraCrs(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> & Op) {
337  try {
338  const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
339  try {
340  const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> & tmp_ECrsMtx = dynamic_cast<const Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> &>(*crsOp.getCrsMatrix());
341  return *tmp_ECrsMtx.getEpetra_CrsMatrix();
342  } catch (std::bad_cast) {
343  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
344  }
345  } catch (std::bad_cast) {
346  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
347  }
348  }
349  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> & Op) {
350  try {
351  Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
352  try {
353  Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> & tmp_ECrsMtx = dynamic_cast<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> &>(*crsOp.getCrsMatrix());
354  return *tmp_ECrsMtx.getEpetra_CrsMatrixNonConst();
355  } catch (std::bad_cast) {
356  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::EpetraCrsMatrix failed");
357  }
358  } catch (std::bad_cast) {
359  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
360  }
361  }
362 
363  static const Epetra_Map& Map2EpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> & map) {
364  RCP<const Xpetra::EpetraMapT<GlobalOrdinal,Node> > xeMap = rcp_dynamic_cast<const Xpetra::EpetraMapT<GlobalOrdinal,Node> >(rcpFromRef(map));
365  if (xeMap == Teuchos::null)
366  throw Exceptions::BadCast("Utilities::Map2EpetraMap : Cast from Xpetra::Map to Xpetra::EpetraMap failed");
367  return xeMap->getEpetra_Map();
368  }
369  // @}
370  #endif
371 
372  #ifdef HAVE_MUELU_TPETRA
373  // @{
375  static RCP<const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > MV2TpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > const vec) {
376  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
377  RCP<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vec);
378  if (tmpVec == Teuchos::null)
379  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
380  return tmpVec->getTpetra_MultiVector();
381  #else
382  throw Exceptions::RuntimeError("MV2TpetraMV: Tpetra has not been compiled with support for LO=GO=int.");
383  #endif
384  }
385  static RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > MV2NonConstTpetraMV(RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > vec) {
386  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
387  RCP<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmpVec = rcp_dynamic_cast<Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(vec);
388  if (tmpVec == Teuchos::null)
389  throw Exceptions::BadCast("Cast from Xpetra::MultiVector to Xpetra::TpetraMultiVector failed");
390  return tmpVec->getTpetra_MultiVector();
391  #else
392  throw Exceptions::RuntimeError("MV2NonConstTpetraMV: Tpetra has not been compiled with support for LO=GO=int.");
393  #endif
394 
395  }
396  static RCP< Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > MV2NonConstTpetraMV2(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & vec) {
397  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
398  const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
399  return tmpVec.getTpetra_MultiVector();
400  #else
401  throw Exceptions::RuntimeError("MV2NonConstTpetraMV2: Tpetra has not been compiled with support for LO=GO=int.");
402  #endif
403  }
404 
405  static const Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& MV2TpetraMV(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & vec) {
406  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
407  const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
408  return *(tmpVec.getTpetra_MultiVector());
409  #else
410  throw Exceptions::RuntimeError("MV2TpetraMV: Tpetra has not been compiled with support for LO=GO=int.");
411  #endif
412  }
413  static Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& MV2NonConstTpetraMV(Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> & vec) {
414  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
415  const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmpVec = dynamic_cast<const Xpetra::TpetraMultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(vec);
416  return *(tmpVec.getTpetra_MultiVector());
417  #else
418  throw Exceptions::RuntimeError("MV2NonConstTpetraMV: Tpetra has not been compiled with support for LO=GO=int.");
419  #endif
420  }
421 
422  static RCP<const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op2TpetraCrs(RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
423  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
424  // Get the underlying Tpetra Mtx
425  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
426  if (crsOp == Teuchos::null)
427  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
428  const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
429  if (tmp_ECrsMtx == Teuchos::null)
430  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
431  return tmp_ECrsMtx->getTpetra_CrsMatrix();
432  #else
433  throw Exceptions::RuntimeError("Op2TpetraCrs: Tpetra has not been compiled with support for LO=GO=int.");
434  #endif
435  }
436  static RCP< Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op2NonConstTpetraCrs(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op){
437  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
438  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
439  if (crsOp == Teuchos::null)
440  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
441  const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > &tmp_ECrsMtx = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsOp->getCrsMatrix());
442  if (tmp_ECrsMtx == Teuchos::null)
443  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
444  return tmp_ECrsMtx->getTpetra_CrsMatrixNonConst();
445  #else
446  throw Exceptions::RuntimeError("Op2NonConstTpetraCrs: Tpetra has not been compiled with support for LO=GO=int.");
447  #endif
448  };
449 
450  static const Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op2TpetraCrs(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> & Op) {
451  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
452  try {
453  const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
454  try {
455  const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
456  return *tmp_ECrsMtx.getTpetra_CrsMatrix();
457  } catch (std::bad_cast) {
458  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
459  }
460  } catch (std::bad_cast) {
461  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
462  }
463  #else
464  throw Exceptions::RuntimeError("Op2TpetraCrs: Tpetra has not been compiled with support for LO=GO=int.");
465  #endif
466  }
467  static Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op2NonConstTpetraCrs(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> & Op) {
468  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
469  try {
470  Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>& crsOp = dynamic_cast<Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(Op);
471  try {
472  Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& tmp_ECrsMtx = dynamic_cast<Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>&>(*crsOp.getCrsMatrix());
473  return *tmp_ECrsMtx.getTpetra_CrsMatrixNonConst();
474  } catch (std::bad_cast) {
475  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix failed");
476  }
477  } catch (std::bad_cast) {
478  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
479  }
480  #else
481  throw Exceptions::RuntimeError("Op2NonConstTpetraCrs: Tpetra has not been compiled with support for LO=GO=int.");
482  #endif
483  }
484 
485  static RCP<const Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op2TpetraRow(RCP<const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
486  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
487  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
488  if (crsOp == Teuchos::null)
489  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
490 
491  RCP<const Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsMat = crsOp->getCrsMatrix();
492  const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
493  RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_BlockCrs;
494  if(!tmp_Crs.is_null()) {
495  return tmp_Crs->getTpetra_CrsMatrixNonConst();
496  }
497  else {
498  tmp_BlockCrs= rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
499  if (tmp_BlockCrs.is_null())
500  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
501  return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
502  }
503  #else
504  throw Exceptions::RuntimeError("Op2TpetraRow: Tpetra has not been compiled with support for LO=GO=int.");
505  #endif
506  }
507  static RCP< Tpetra::RowMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op2NonConstTpetraRow(RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) {
508  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
509  RCP<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsOp = rcp_dynamic_cast<const Xpetra::CrsMatrixWrap<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(Op);
510  if (crsOp == Teuchos::null)
511  throw Exceptions::BadCast("Cast from Xpetra::Matrix to Xpetra::CrsMatrixWrap failed");
512 
513  RCP<const Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > crsMat = crsOp->getCrsMatrix();
514  const RCP<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_Crs = rcp_dynamic_cast<const Xpetra::TpetraCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
515  RCP<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > tmp_BlockCrs;
516  if(!tmp_Crs.is_null()) {
517  return tmp_Crs->getTpetra_CrsMatrixNonConst();
518  }
519  else {
520  tmp_BlockCrs= rcp_dynamic_cast<const Xpetra::TpetraBlockCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(crsMat);
521  if (tmp_BlockCrs.is_null())
522  throw Exceptions::BadCast("Cast from Xpetra::CrsMatrix to Xpetra::TpetraCrsMatrix and Xpetra::TpetraBlockCrsMatrix failed");
523  return tmp_BlockCrs->getTpetra_BlockCrsMatrixNonConst();
524  }
525  #else
526  throw Exceptions::RuntimeError("Op2NonConstTpetraRow: Tpetra has not been compiled with support for LO=GO=int.");
527  #endif
528  };
529 
530 
531  static const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > Map2TpetraMap(const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> & map) {
532  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
533  const RCP<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node>>& tmp_TMap = rcp_dynamic_cast<const Xpetra::TpetraMap<LocalOrdinal,GlobalOrdinal,Node> >(rcpFromRef(map));
534  if (tmp_TMap == Teuchos::null)
535  throw Exceptions::BadCast("Utilities::Map2TpetraMap : Cast from Xpetra::Map to Xpetra::TpetraMap failed");
536  return tmp_TMap->getTpetra_Map();
537  #else
538  throw Exceptions::RuntimeError("Map2TpetraMap: Tpetra has not been compiled with support for LO=GO=int.");
539  #endif
540  };
541  #endif
542 
543  static RCP<Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Crs2Op(RCP<Xpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > Op) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Crs2Op(Op); }
544  static Teuchos::ArrayRCP<Scalar> GetMatrixDiagonal(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetMatrixDiagonal(A); }
545  static RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > GetMatrixDiagonalInverse(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, Magnitude tol = Teuchos::ScalarTraits<Scalar>::eps()*100) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetMatrixDiagonalInverse(A,tol); }
546  static Teuchos::ArrayRCP<Scalar> GetLumpedMatrixDiagonal(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetLumpedMatrixDiagonal(A); }
547  static RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > GetMatrixOverlappedDiagonal(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::GetMatrixOverlappedDiagonal(A); }
548  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) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::ResidualNorm(Op,X,RHS); }
549  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) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Residual(Op,X,RHS); }
551  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::MakeFancy(os); }
552  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>& v, LocalOrdinal i0, LocalOrdinal i1) { return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Distance2(v,i0,i1); }
554 
555  // todo: move this to UtilitiesBase::kokkos
556  static Kokkos::View<const bool*, typename Node::device_type> DetectDirichletRows(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, const Magnitude& tol = Teuchos::ScalarTraits<Scalar>::zero());
557 
558  static Scalar PowerMethod(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A, bool scaleByDiag = true,
559  LocalOrdinal niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
560  return MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node>::PowerMethod(A,scaleByDiag,niters,tolerance,verbose,seed);
561  }
562 
563  static void MyOldScaleMatrix(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const Teuchos::ArrayRCP<const Scalar>& scalingVector, bool doInverse = true,
564  bool doFillComplete = true, bool doOptimizeStorage = true) {
565  Scalar one = Teuchos::ScalarTraits<Scalar>::one();
566  Teuchos::ArrayRCP<Scalar> sv(scalingVector.size());
567  if (doInverse) {
568  for (int i = 0; i < scalingVector.size(); ++i)
569  sv[i] = one / scalingVector[i];
570  } else {
571  for (int i = 0; i < scalingVector.size(); ++i)
572  sv[i] = scalingVector[i];
573  }
574 
575  switch (Op.getRowMap()->lib()) {
576  case Xpetra::UseTpetra:
577  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
578  break;
579 
580  case Xpetra::UseEpetra:
581  MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
582  break;
583 
584  default:
585  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
586  break;
587  }
588  }
589 
590  // TODO This is the <double,int,int> specialization
591  static void MyOldScaleMatrix_Tpetra(Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector,
592  bool doFillComplete, bool doOptimizeStorage) {
593  #ifdef HAVE_MUELU_TPETRA
594  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
595  try {
596  Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpOp = Op2NonConstTpetraCrs(Op);
597 
598  const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = tpOp.getRowMap();
599  const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domainMap = tpOp.getDomainMap();
600  const RCP<const Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rangeMap = tpOp.getRangeMap();
601 
602  size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
603  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
604  maxRowSize = 20;
605 
606  std::vector<Scalar> scaledVals(maxRowSize);
607  if (tpOp.isFillComplete())
608  tpOp.resumeFill();
609 
610  if (Op.isLocallyIndexed() == true) {
611  Teuchos::ArrayView<const LocalOrdinal> cols;
612  Teuchos::ArrayView<const Scalar> vals;
613 
614  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
615  tpOp.getLocalRowView(i, cols, vals);
616  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
617  if (nnz > maxRowSize) {
618  maxRowSize = nnz;
619  scaledVals.resize(maxRowSize);
620  }
621  for (size_t j = 0; j < nnz; ++j)
622  scaledVals[j] = vals[j]*scalingVector[i];
623 
624  if (nnz > 0) {
625  Teuchos::ArrayView<const Scalar> valview(&scaledVals[0], nnz);
626  tpOp.replaceLocalValues(i, cols, valview);
627  }
628  } //for (size_t i=0; ...
629 
630  } else {
631  Teuchos::ArrayView<const GlobalOrdinal> cols;
632  Teuchos::ArrayView<const Scalar> vals;
633 
634  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
635  GlobalOrdinal gid = rowMap->getGlobalElement(i);
636  tpOp.getGlobalRowView(gid, cols, vals);
637  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
638  if (nnz > maxRowSize) {
639  maxRowSize = nnz;
640  scaledVals.resize(maxRowSize);
641  }
642  // FIXME FIXME FIXME FIXME FIXME FIXME
643  for (size_t j = 0; j < nnz; ++j)
644  scaledVals[j] = vals[j]*scalingVector[i]; //FIXME i or gid?
645 
646  if (nnz > 0) {
647  Teuchos::ArrayView<const Scalar> valview(&scaledVals[0], nnz);
648  tpOp.replaceGlobalValues(gid, cols, valview);
649  }
650  } //for (size_t i=0; ...
651  }
652 
653  if (doFillComplete) {
654  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
655  throw Exceptions::RuntimeError("In Utilities::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
656 
657  RCP<Teuchos::ParameterList> params = rcp(new Teuchos::ParameterList());
658  params->set("Optimize Storage", doOptimizeStorage);
659  params->set("No Nonlocal Changes", true);
660  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
661  }
662  } catch(...) {
663  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
664  }
665  #else
666  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been compiled with support for LO=GO=int.");
667  #endif
668  #else
669  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
670  #endif
671  }
672 
673  static void MyOldScaleMatrix_Epetra (Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector, bool doFillComplete, bool doOptimizeStorage) {
674  #ifdef HAVE_MUELU_EPETRA
675  try {
676  //const Epetra_CrsMatrix& epOp = Utilities<double,int,int>::Op2NonConstEpetraCrs(Op);
677  const Epetra_CrsMatrix& epOp = Op2NonConstEpetraCrs(Op);
678 
679  Epetra_Map const &rowMap = epOp.RowMap();
680  int nnz;
681  double *vals;
682  int *cols;
683 
684  for (int i = 0; i < rowMap.NumMyElements(); ++i) {
685  epOp.ExtractMyRowView(i, nnz, vals, cols);
686  for (int j = 0; j < nnz; ++j)
687  vals[j] *= scalingVector[i];
688  }
689 
690  } catch (...){
691  throw Exceptions::RuntimeError("Only Epetra_CrsMatrix types can be scaled");
692  }
693  #else
694  throw Exceptions::RuntimeError("Matrix scaling is not possible because Epetra has not been enabled.");
695  #endif // HAVE_MUELU_EPETRA
696  }
697 
703  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()) {
704  switch (Op.getRowMap()->lib()) {
705  case Xpetra::UseTpetra:
706  {
707  #ifdef HAVE_MUELU_TPETRA
708  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
709  try {
710  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpetraOp = Utilities<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Op2TpetraCrs(Op);
711 
712  // Compute the transpose A of the Tpetra matrix tpetraOp.
713  RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
714  Tpetra::RowMatrixTransposer<Scalar, LocalOrdinal, GlobalOrdinal, Node> transposer(rcpFromRef(tpetraOp),label);
715  A = transposer.createTranspose();
716  RCP<Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AA = rcp(new Xpetra::TpetraCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>(A));
717  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAA = rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(AA);
718  RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAAA = rcp( new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> (AAA));
719 
720  return AAAA;
721  }
722  catch (std::exception& e) {
723  std::cout << "threw exception '" << e.what() << "'" << std::endl;
724  throw Exceptions::RuntimeError("Utilities::Transpose failed, perhaps because matrix is not a Crs matrix");
725  }
726  #else
727  throw Exceptions::RuntimeError("Utilities::Transpose: Tpetra is not compiled with LO=GO=int. Add TPETRA_INST_INT_INT:BOOL=ON to your configuration!");
728  #endif
729  #else
730  throw Exceptions::RuntimeError("Utilities::Transpose: Tpetra is not compiled!");
731  #endif
732  break;
733  }
734  case Xpetra::UseEpetra:
735  {
736  #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
737  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("ZZ Entire Transpose"));
738  // Epetra case
741  Epetra_CrsMatrix * A = dynamic_cast<Epetra_CrsMatrix*>(&transposer(epetraOp));
742  transposer.ReleaseTranspose(); // So we can keep A in Muelu...
743 
744  RCP<Epetra_CrsMatrix> rcpA(A);
745  RCP<Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> > AA = rcp(new Xpetra::EpetraCrsMatrixT<GlobalOrdinal,Node> (rcpA));
746  RCP<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAA = rcp_implicit_cast<Xpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(AA);
747  RCP<Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> > AAAA = rcp( new Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node>(AAA));
748  AAAA->fillComplete(Op.getRangeMap(), Op.getDomainMap());
749 
750  return AAAA;
751  #else
752  throw Exceptions::RuntimeError("Epetra (Err. 2)");
753  #endif
754  break;
755  }
756  default:
757  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be transposed.");
758  break;
759  }
760 
761  return Teuchos::null;
762  }
763 
764 
765  }; // class Utilities (specialization SC=double LO=GO=int)
766 
767 
768 
769  // Useful Kokkos conversions
770  template < class View, unsigned AppendValue >
771  struct AppendTrait {
772  // static_assert(false, "Error: NOT a Kokkos::View");
773  };
774 
775  // Arg3 == MemoryTraits
776  template < class DataType, class Arg1, class Arg2, unsigned U, unsigned T >
777  struct AppendTrait< Kokkos::View< DataType, Arg1, Arg2, Kokkos::MemoryTraits<U> >, T> {
779  };
780 
781  // Arg2 == MemoryTraits
782  template < class DataType, class Arg1, unsigned U, unsigned T >
783  struct AppendTrait< Kokkos::View< DataType, Arg1, Kokkos::MemoryTraits<U>, void >, T> {
785  };
786 
787 
788  // Arg1 == MemoryTraits
789  template < class DataType, unsigned U, unsigned T >
790  struct AppendTrait< Kokkos::View< DataType, Kokkos::MemoryTraits<U>, void, void >, T> {
791  using type = Kokkos::View< DataType, Kokkos::MemoryTraits<U|T>, void, void >;
792  };
793 
794  // 2 arguments -- no traits
795  template < class DataType, class Arg1, class Arg2, unsigned T >
796  struct AppendTrait< Kokkos::View< DataType, Arg1, Arg2, void >, T> {
798  };
799 
800  // 1 arguments -- no traits
801  template < class DataType, class Arg1, unsigned T >
802  struct AppendTrait< Kokkos::View< DataType, Arg1, void, void >, T> {
804  };
805 
806  // 0 arguments
807  template < class DataType, unsigned T >
808  struct AppendTrait< Kokkos::View< DataType, void, void, void >, T> {
809  using type = Kokkos::View< DataType, Kokkos::MemoryTraits<T>, void, void >;
810  };
811 
812 
813 } //namespace MueLu
814 
815 #define MUELU_UTILITIES_KOKKOS_SHORT
816 
817 #endif // #if defined(HAVE_MUELU_KOKKOS_REFACTOR)
818 
819 #endif // MUELU_UTILITIES_KOKKOS_DECL_HPP
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)
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)
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())
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)
Simple transpose for Tpetra::CrsMatrix types.
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.
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &v, LocalOrdinal i0, LocalOrdinal i1)
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
Namespace for MueLu classes and methods.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
Kokkos::View< const bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol)
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)
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)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
static RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static Teuchos::ArrayRCP< Scalar > GetLumpedMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal of lumped matrix.
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)
static RCP< Matrix > Crs2Op(RCP< CrsMatrix > Op)
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.
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 RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100)
Extract Matrix Diagonal.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.