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>
56 
58 #include <Xpetra_CrsMatrix_fwd.hpp>
60 #include <Xpetra_ExportFactory.hpp>
63 #include <Xpetra_Map_fwd.hpp>
65 #include <Xpetra_Matrix_fwd.hpp>
68 #include <Xpetra_Operator_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>
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>
91 #endif
92 
93 
94 namespace MueLu {
95 
103  template <class Scalar,
104  class LocalOrdinal = int,
105  class GlobalOrdinal = LocalOrdinal,
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 
184 
192  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A); // FIXME
193 
202 
203  // TODO: should NOT return an Array. Definition must be changed to:
204  // - ArrayRCP<> ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS)
205  // or
206  // - void ResidualNorm(Matrix const &Op, MultiVector const &X, MultiVector const &RHS, Array &)
207  static Teuchos::Array<Magnitude> ResidualNorm(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
208  return Utilities::ResidualNorm(Op, X, RHS);
209  }
210 
211  static RCP<MultiVector> Residual(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
212  return Utilities::Residual(Op, X, RHS);
213  }
214 
215  static void PauseForDebugger();
216 
232  static SC PowerMethod(const Matrix& A, bool scaleByDiag = true,
233  LO niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
234  return Utilities::PowerMethod(A, scaleByDiag, niters, tolerance, verbose, seed);
235  }
236 
237  static void MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector, bool doInverse = true,
238  bool doFillComplete = true, bool doOptimizeStorage = true); // FIXME
239 
240  static void MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
241  bool doFillComplete, bool doOptimizeStorage); // FIXME
242 
243  static void MyOldScaleMatrix_Epetra(Matrix& Op, const Teuchos::ArrayRCP<SC>& scalingVector,
244  bool doFillComplete, bool doOptimizeStorage); // FIXME
245 
246  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) { return Utilities::MakeFancy(os); }
247 
252  static typename Teuchos::ScalarTraits<Scalar>::magnitudeType Distance2(const MultiVector& v, LocalOrdinal i0, LocalOrdinal i1) { return Utilities::Distance2(v, i0, i1); }
253 
261  static Kokkos::View<const bool*, typename NO::device_type> DetectDirichletRows(const Matrix& A, const Magnitude& tol = Teuchos::ScalarTraits<SC>::zero());
262 
272  static void SetRandomSeed(const Teuchos::Comm<int> &comm) { return Utilities::SetRandomSeed(comm); }
273 
279  static RCP<Matrix> Transpose(Matrix& Op, bool optimizeTranspose = false, const std::string & label = std::string()) {
280  return Utilities::Transpose(Op, optimizeTranspose, label);
281  }
282 
283  static RCP<Xpetra::MultiVector<double,LocalOrdinal,GlobalOrdinal,Node> > ExtractCoordinatesFromParameterList(ParameterList& paramList) {
285  }
286 
287  }; // class Utils
288 
289 
299  template <class Node>
300  class Utilities_kokkos<double,int,int,Node> : public UtilitiesBase<double,int,int,Node> {
301  public:
302  typedef double Scalar;
303  typedef int LocalOrdinal;
304  typedef int GlobalOrdinal;
305  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
306 
307  private:
308 #undef MUELU_UTILITIES_KOKKOS_SHORT
309 #include "MueLu_UseShortNames.hpp"
310 
311  public:
312 
313 #ifdef HAVE_MUELU_EPETRA
314  // @{
316  static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2EpetraMV(vec); }
317  static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstEpetraMV(vec); }
318 
319  static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec) { return Utilities::MV2EpetraMV(vec); }
320  static Epetra_MultiVector& MV2NonConstEpetraMV(MultiVector& vec) { return Utilities::MV2NonConstEpetraMV(vec); }
321 
322  static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2EpetraCrs(Op); }
323  static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
324 
325  static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) { return Utilities::Op2EpetraCrs(Op); }
326  static Epetra_CrsMatrix& Op2NonConstEpetraCrs(Matrix& Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
327 
328  static const Epetra_Map& Map2EpetraMap(const Map& map) { return Utilities::Map2EpetraMap(map); }
329  // @}
330 #endif
331 
332 #ifdef HAVE_MUELU_TPETRA
333  // @{
335  static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2TpetraMV(vec); }
336  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstTpetraMV(vec); }
337  static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV2(vec); }
338 
339  static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec) { return Utilities::MV2TpetraMV(vec); }
340  static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV(vec); }
341 
342  static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2TpetraCrs(Op); }
343  static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
344 
345  static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op) { return Utilities::Op2TpetraCrs(Op); }
346  static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
347 
348  static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op) { return Utilities::Op2TpetraRow(Op); }
349  static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraRow(Op); }
350 
351  static const RCP<const Tpetra::Map<LO, GO, NO> > Map2TpetraMap(const Map& map) { return Utilities::Map2TpetraMap(map); }
352 #endif
353  static RCP<Xpetra::Matrix<SC,LO,GO,NO> > Crs2Op(RCP<CrsMatrix> Op) { return Utilities::Crs2Op(Op); }
354 
355  static ArrayRCP<SC> GetMatrixDiagonal(const Matrix& A) {
357  }
358  static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<SC>::eps()*100) {
360  }
361  static ArrayRCP<SC> GetLumpedMatrixDiagonal(const Matrix& A) {
363  }
364  static RCP<Vector> GetLumpedMatrixDiagonal(RCP<const Matrix > A) {
366  }
367  static RCP<Vector> GetMatrixOverlappedDiagonal(const Matrix& A) {
369  }
370  static RCP<Vector> GetInverse(RCP<const Vector> v, Magnitude tol = Teuchos::ScalarTraits<SC>::eps()*100, SC tolReplacement = Teuchos::ScalarTraits<SC>::zero()) {
371  return UtilitiesBase::GetInverse(v,tol,tolReplacement);
372  }
373  static Teuchos::Array<Magnitude> ResidualNorm(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
374  return UtilitiesBase::ResidualNorm(Op,X,RHS);
375  }
376  static RCP<MultiVector> Residual(const Operator& Op, const MultiVector& X, const MultiVector& RHS) {
377  return UtilitiesBase::Residual(Op,X,RHS);
378  }
379  static void PauseForDebugger() {
381  }
382  static RCP<Teuchos::FancyOStream> MakeFancy(std::ostream& os) {
383  return UtilitiesBase::MakeFancy(os);
384  }
385  static typename Teuchos::ScalarTraits<SC>::magnitudeType Distance2(const MultiVector& v, LO i0, LO i1) {
386  return UtilitiesBase::Distance2(v,i0,i1);
387  }
388  static void SetRandomSeed(const Teuchos::Comm<int> &comm) {
390  }
391 
392  // todo: move this to UtilitiesBase::kokkos
393  static Kokkos::View<const bool*, typename Node::device_type> DetectDirichletRows(const Matrix& A, const Magnitude& tol = Teuchos::ScalarTraits<SC>::zero());
394 
395  static Scalar PowerMethod(const Matrix& A, bool scaleByDiag = true, LO niters = 10, Magnitude tolerance = 1e-2, bool verbose = false, unsigned int seed = 123) {
396  return UtilitiesBase::PowerMethod(A,scaleByDiag,niters,tolerance,verbose,seed);
397  }
398 
399  static void MyOldScaleMatrix(Matrix& Op, const Teuchos::ArrayRCP<const SC>& scalingVector, bool doInverse = true,
400  bool doFillComplete = true, bool doOptimizeStorage = true) {
402  Teuchos::ArrayRCP<SC> sv(scalingVector.size());
403  if (doInverse) {
404  for (int i = 0; i < scalingVector.size(); ++i)
405  sv[i] = one / scalingVector[i];
406  } else {
407  for (int i = 0; i < scalingVector.size(); ++i)
408  sv[i] = scalingVector[i];
409  }
410 
411  switch (Op.getRowMap()->lib()) {
412  case Xpetra::UseTpetra:
413  MyOldScaleMatrix_Tpetra(Op, sv, doFillComplete, doOptimizeStorage);
414  break;
415 
416  case Xpetra::UseEpetra:
417  MyOldScaleMatrix_Epetra(Op, sv, doFillComplete, doOptimizeStorage);
418  break;
419 
420  default:
421  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be scaled.");
422  break;
423  }
424  }
425 
426  // TODO This is the <double,int,int> specialization
427  static void MyOldScaleMatrix_Tpetra(Matrix& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector,
428  bool doFillComplete, bool doOptimizeStorage) {
429  #ifdef HAVE_MUELU_TPETRA
430  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
431  try {
432  Tpetra::CrsMatrix<SC,LO,GO,NO>& tpOp = Op2NonConstTpetraCrs(Op);
433 
434  const RCP<const Tpetra::Map<LO,GO,NO> > rowMap = tpOp.getRowMap();
435  const RCP<const Tpetra::Map<LO,GO,NO> > domainMap = tpOp.getDomainMap();
436  const RCP<const Tpetra::Map<LO,GO,NO> > rangeMap = tpOp.getRangeMap();
437 
438  size_t maxRowSize = tpOp.getNodeMaxNumRowEntries();
439  if (maxRowSize == Teuchos::as<size_t>(-1)) // hasn't been determined yet
440  maxRowSize = 20;
441 
442  std::vector<SC> scaledVals(maxRowSize);
443  if (tpOp.isFillComplete())
444  tpOp.resumeFill();
445 
446  if (Op.isLocallyIndexed() == true) {
449 
450  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
451  tpOp.getLocalRowView(i, cols, vals);
452  size_t nnz = tpOp.getNumEntriesInLocalRow(i);
453  if (nnz > maxRowSize) {
454  maxRowSize = nnz;
455  scaledVals.resize(maxRowSize);
456  }
457  for (size_t j = 0; j < nnz; ++j)
458  scaledVals[j] = vals[j]*scalingVector[i];
459 
460  if (nnz > 0) {
461  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
462  tpOp.replaceLocalValues(i, cols, valview);
463  }
464  } //for (size_t i=0; ...
465 
466  } else {
469 
470  for (size_t i = 0; i < rowMap->getNodeNumElements(); ++i) {
471  GO gid = rowMap->getGlobalElement(i);
472  tpOp.getGlobalRowView(gid, cols, vals);
473  size_t nnz = tpOp.getNumEntriesInGlobalRow(gid);
474  if (nnz > maxRowSize) {
475  maxRowSize = nnz;
476  scaledVals.resize(maxRowSize);
477  }
478  // FIXME FIXME FIXME FIXME FIXME FIXME
479  for (size_t j = 0; j < nnz; ++j)
480  scaledVals[j] = vals[j]*scalingVector[i]; //FIXME i or gid?
481 
482  if (nnz > 0) {
483  Teuchos::ArrayView<const SC> valview(&scaledVals[0], nnz);
484  tpOp.replaceGlobalValues(gid, cols, valview);
485  }
486  } //for (size_t i=0; ...
487  }
488 
489  if (doFillComplete) {
490  if (domainMap == Teuchos::null || rangeMap == Teuchos::null)
491  throw Exceptions::RuntimeError("In Utilities::Scaling: cannot fillComplete because the domain and/or range map hasn't been defined");
492 
494  params->set("Optimize Storage", doOptimizeStorage);
495  params->set("No Nonlocal Changes", true);
496  Op.fillComplete(Op.getDomainMap(), Op.getRangeMap(), params);
497  }
498  } catch(...) {
499  throw Exceptions::RuntimeError("Only Tpetra::CrsMatrix types can be scaled (Err.1)");
500  }
501  #else
502  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been compiled with support for LO=GO=int.");
503  #endif
504  #else
505  throw Exceptions::RuntimeError("Matrix scaling is not possible because Tpetra has not been enabled.");
506  #endif
507  }
508 
509  static void MyOldScaleMatrix_Epetra(Matrix& Op, const Teuchos::ArrayRCP<Scalar>& scalingVector, bool doFillComplete, bool doOptimizeStorage) {
510  #ifdef HAVE_MUELU_EPETRA
511  try {
512  //const Epetra_CrsMatrix& epOp = Utilities<double,int,int>::Op2NonConstEpetraCrs(Op);
513  const Epetra_CrsMatrix& epOp = Op2NonConstEpetraCrs(Op);
514 
515  Epetra_Map const &rowMap = epOp.RowMap();
516  int nnz;
517  double *vals;
518  int *cols;
519 
520  for (int i = 0; i < rowMap.NumMyElements(); ++i) {
521  epOp.ExtractMyRowView(i, nnz, vals, cols);
522  for (int j = 0; j < nnz; ++j)
523  vals[j] *= scalingVector[i];
524  }
525 
526  } catch (...){
527  throw Exceptions::RuntimeError("Only Epetra_CrsMatrix types can be scaled");
528  }
529  #else
530  throw Exceptions::RuntimeError("Matrix scaling is not possible because Epetra has not been enabled.");
531  #endif // HAVE_MUELU_EPETRA
532  }
533 
539  static RCP<Matrix> Transpose(Matrix& Op, bool optimizeTranspose = false, const std::string & label = std::string()) {
540  switch (Op.getRowMap()->lib()) {
541  case Xpetra::UseTpetra:
542  {
543  #ifdef HAVE_MUELU_TPETRA
544  #ifdef HAVE_MUELU_TPETRA_INST_INT_INT
545  try {
546  const Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>& tpetraOp = Utilities::Op2TpetraCrs(Op);
547 
548  // Compute the transpose A of the Tpetra matrix tpetraOp.
549  RCP<Tpetra::CrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > A;
550  Tpetra::RowMatrixTransposer<SC,LO,GO,NO> transposer(rcpFromRef(tpetraOp),label);
551  A = transposer.createTranspose();
552  RCP<Xpetra::TpetraCrsMatrix<SC,LO,GO,NO> > AA = rcp(new Xpetra::TpetraCrsMatrix<SC,LO,GO,NO>(A));
553  RCP<CrsMatrix> AAA = rcp_implicit_cast<CrsMatrix>(AA);
554  RCP<CrsMatrixWrap> AAAA = rcp(new CrsMatrixWrap(AAA));
555 
556  return AAAA;
557  }
558  catch (std::exception& e) {
559  std::cout << "threw exception '" << e.what() << "'" << std::endl;
560  throw Exceptions::RuntimeError("Utilities::Transpose failed, perhaps because matrix is not a Crs matrix");
561  }
562  #else
563  throw Exceptions::RuntimeError("Utilities::Transpose: Tpetra is not compiled with LO=GO=int. Add TPETRA_INST_INT_INT:BOOL=ON to your configuration!");
564  #endif
565  #else
566  throw Exceptions::RuntimeError("Utilities::Transpose: Tpetra is not compiled!");
567  #endif
568  break;
569  }
570  case Xpetra::UseEpetra:
571  {
572  #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_EPETRAEXT)
573  Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer("ZZ Entire Transpose"));
574  // Epetra case
576  EpetraExt::RowMatrix_Transpose transposer;
577  Epetra_CrsMatrix * A = dynamic_cast<Epetra_CrsMatrix*>(&transposer(epetraOp));
578  transposer.ReleaseTranspose(); // So we can keep A in Muelu...
579 
580  RCP<Epetra_CrsMatrix> rcpA(A);
581  RCP<Xpetra::EpetraCrsMatrixT<GO,NO> > AA = rcp(new Xpetra::EpetraCrsMatrixT<GO,NO> (rcpA));
582  RCP<CrsMatrix> AAA = rcp_implicit_cast<CrsMatrix>(AA);
583  RCP<CrsMatrixWrap> AAAA = rcp( new CrsMatrixWrap(AAA));
584  AAAA->fillComplete(Op.getRangeMap(), Op.getDomainMap());
585 
586  return AAAA;
587  #else
588  throw Exceptions::RuntimeError("Epetra (Err. 2)");
589  #endif
590  break;
591  }
592  default:
593  throw Exceptions::RuntimeError("Only Epetra and Tpetra matrices can be transposed.");
594  break;
595  }
596 
597  return Teuchos::null;
598  }
599 
602  static RCP<Xpetra::MultiVector<double,LO,GO,NO> > ExtractCoordinatesFromParameterList(ParameterList& paramList) {
603  RCP<Xpetra::MultiVector<double,LO,GO,NO> > coordinates = Teuchos::null;
604 
605  // check whether coordinates are contained in parameter list
606  if(paramList.isParameter ("Coordinates") == false)
607  return coordinates;
608 
609  #if defined(HAVE_MUELU_TPETRA)
610  #if ( defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_OPENMP) && defined(HAVE_TPETRA_INST_INT_INT)) || \
611  (!defined(EPETRA_HAVE_OMP) && defined(HAVE_TPETRA_INST_SERIAL) && defined(HAVE_TPETRA_INST_INT_INT))
612 
613  // define Tpetra::MultiVector type with Scalar=float only if
614  // * ETI is turned off, since then the compiler will instantiate it automatically OR
615  // * Tpetra is instantiated on Scalar=float
616  #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
617  typedef Tpetra::MultiVector<float, LO,GO,NO> tfMV;
618  RCP<tfMV> floatCoords = Teuchos::null;
619  #endif
620 
621  // define Tpetra::MultiVector type with Scalar=double only if
622  // * ETI is turned off, since then the compiler will instantiate it automatically OR
623  // * Tpetra is instantiated on Scalar=double
624  typedef Tpetra::MultiVector<SC,LO,GO,NO> tdMV;
625  RCP<tdMV> doubleCoords = Teuchos::null;
626  if (paramList.isType<RCP<tdMV> >("Coordinates")) {
627  // Coordinates are stored as a double vector
628  doubleCoords = paramList.get<RCP<tdMV> >("Coordinates");
629  paramList.remove("Coordinates");
630  }
631  #if !defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) || defined(HAVE_TPETRA_INST_FLOAT)
632  else if (paramList.isType<RCP<tfMV> >("Coordinates")) {
633  // check if coordinates are stored as a float vector
634  floatCoords = paramList.get<RCP<tfMV> >("Coordinates");
635  paramList.remove("Coordinates");
636  doubleCoords = rcp(new tdMV(floatCoords->getMap(), floatCoords->getNumVectors()));
637  deep_copy(*doubleCoords, *floatCoords);
638  }
639  #endif
640  // We have the coordinates in a Tpetra double vector
641  if(doubleCoords != Teuchos::null) {
642  coordinates = Teuchos::rcp(new Xpetra::TpetraMultiVector<SC,LO,GO,NO>(doubleCoords));
643  TEUCHOS_TEST_FOR_EXCEPT(doubleCoords->getNumVectors() != coordinates->getNumVectors());
644  }
645  #endif // Tpetra instantiated on GO=int and EpetraNode
646  #endif // endif HAVE_TPETRA
647 
648  #if defined(HAVE_MUELU_EPETRA)
649  RCP<Epetra_MultiVector> doubleEpCoords;
650  if (paramList.isType<RCP<Epetra_MultiVector> >("Coordinates")) {
651  doubleEpCoords = paramList.get<RCP<Epetra_MultiVector> >("Coordinates");
652  paramList.remove("Coordinates");
653  RCP<Xpetra::EpetraMultiVectorT<GO,NO> > epCoordinates = Teuchos::rcp(new Xpetra::EpetraMultiVectorT<GO,NO>(doubleEpCoords));
654  coordinates = rcp_dynamic_cast<Xpetra::MultiVector<double,LO,GO,NO> >(epCoordinates);
655  TEUCHOS_TEST_FOR_EXCEPT(doubleEpCoords->NumVectors() != Teuchos::as<int>(coordinates->getNumVectors()));
656  }
657  #endif
659  return coordinates;
660  }
661  }; // class Utilities (specialization SC=double LO=GO=int)
662 
663 
664 
665  // Useful Kokkos conversions
666  template < class View, unsigned AppendValue >
667  struct AppendTrait {
668  // static_assert(false, "Error: NOT a Kokkos::View");
669  };
670 
671  template < class MT, unsigned T >
672  struct CombineMemoryTraits {
673  // empty
674  };
675 
676  template < unsigned U, unsigned T>
677  struct CombineMemoryTraits<Kokkos::MemoryTraits<U>, T> {
678  typedef Kokkos::MemoryTraits<U|T> type;
679  };
680 
681  template < class DataType, unsigned T, class... Pack >
682  struct AppendTrait< Kokkos::View< DataType, Pack... >, T> {
683  typedef Kokkos::View< DataType, Pack... > view_type;
684  using type = Kokkos::View< DataType, typename view_type::array_layout, typename view_type::device_type, typename CombineMemoryTraits<typename view_type::memory_traits,T>::type >;
685  };
686 
687 } //namespace MueLu
688 
689 #define MUELU_UTILITIES_KOKKOS_SHORT
690 
691 #endif // #if defined(HAVE_MUELU_KOKKOS_REFACTOR)
692 
693 #endif // MUELU_UTILITIES_KOKKOS_DECL_HPP
const Epetra_Map & RowMap() const
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)
size_type size() const
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LocalOrdinal niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Simple transpose for Tpetra::CrsMatrix types.
GlobalOrdinal GO
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool is_null(const std::shared_ptr< T > &p)
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::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
LocalOrdinal LO
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.
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
static RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Time > getNewTimer(const std::string &name)
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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
View
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
int NumMyElements() const
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< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100)
Extract Matrix Diagonal.
Scalar SC
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 void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.