MueLu Version of the Day
Loading...
Searching...
No Matches
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
51#include <string>
52
53#include "Teuchos_ScalarTraits.hpp"
54#include "Teuchos_ParameterList.hpp"
55
56#include "Xpetra_BlockedCrsMatrix_fwd.hpp"
57#include "Xpetra_CrsMatrix_fwd.hpp"
58#include "Xpetra_CrsMatrixWrap_fwd.hpp"
59#include "Xpetra_ExportFactory.hpp"
60#include "Xpetra_ImportFactory_fwd.hpp"
61#include "Xpetra_MapFactory_fwd.hpp"
62#include "Xpetra_Map_fwd.hpp"
63#include "Xpetra_MatrixFactory_fwd.hpp"
64#include "Xpetra_Matrix_fwd.hpp"
65#include "Xpetra_MultiVectorFactory_fwd.hpp"
66#include "Xpetra_MultiVector_fwd.hpp"
67#include "Xpetra_Operator_fwd.hpp"
68#include "Xpetra_VectorFactory_fwd.hpp"
69#include "Xpetra_Vector_fwd.hpp"
70
71#include "Xpetra_IO.hpp"
72
73#include "Kokkos_ArithTraits.hpp"
74
75#ifdef HAVE_MUELU_EPETRA
76#include "Epetra_MultiVector.h"
77#include "Epetra_CrsMatrix.h"
78#include "Xpetra_EpetraCrsMatrix_fwd.hpp"
79#include "Xpetra_EpetraMultiVector_fwd.hpp"
80#endif
81
82#include "MueLu_Exceptions.hpp"
83#include "MueLu_Utilities.hpp"
84#include "MueLu_UtilitiesBase.hpp"
85
86#ifdef HAVE_MUELU_TPETRA
87#include "Tpetra_CrsMatrix.hpp"
88#include "Tpetra_Map.hpp"
89#include "Tpetra_MultiVector.hpp"
90#include "Xpetra_TpetraCrsMatrix_fwd.hpp"
91#include "Xpetra_TpetraMultiVector_fwd.hpp"
92#endif
93
94
95namespace MueLu {
96
104 template <class Scalar,
107 class Node = DefaultNode>
108 class Utilities_kokkos : public MueLu::UtilitiesBase<Scalar,LocalOrdinal,GlobalOrdinal,Node> {
109#undef MUELU_UTILITIES_KOKKOS_SHORT
111
112 public:
113 using TST = Teuchos::ScalarTraits<SC>;
114 using Magnitude = typename TST::magnitudeType;
115 using CoordinateType = typename TST::coordinateType;
116 using RealValuedMultiVector = Xpetra::MultiVector<CoordinateType,LO,GO,NO>;
117
118#ifdef HAVE_MUELU_EPETRA
120 // @{
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); }
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); }
132
133 static const Epetra_Map& Map2EpetraMap(const Map& map) { return Utilities::Map2EpetraMap(map); }
134 // @}
135#endif
136
137#ifdef HAVE_MUELU_TPETRA
139 // @{
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<typename Teuchos::ScalarTraits<SC>::magnitudeType>::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,
262 Kokkos::View<bool*, typename Node::device_type> nonzeros);
263
264
275 static Kokkos::View<bool*, typename NO::device_type> DetectDirichletCols(const Matrix& A, const Kokkos::View<const bool*, typename NO::device_type>& dirichletRows);
276
277
285 static void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
286 const Kokkos::View<bool*, typename Node::device_type> & dirichletRows,
287 Kokkos::View<bool*, typename Node::device_type> dirichletCols,
288 Kokkos::View<bool*, typename Node::device_type> dirichletDomain);
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,
299 Kokkos::View<bool*, typename NO::device_type> & dirichletRows);
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 using CoordinateType = typename Teuchos::ScalarTraits<Scalar>::coordinateType;
358 using RealValuedMultiVector = Xpetra::MultiVector<CoordinateType,LocalOrdinal,GlobalOrdinal,Node>;
359
360 private:
361#undef MUELU_UTILITIES_KOKKOS_SHORT
363
364 public:
365
366#ifdef HAVE_MUELU_EPETRA
368 // @{
369 static RCP<const Epetra_MultiVector> MV2EpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2EpetraMV(vec); }
370 static RCP< Epetra_MultiVector> MV2NonConstEpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstEpetraMV(vec); }
371
372 static const Epetra_MultiVector& MV2EpetraMV(const MultiVector& vec) { return Utilities::MV2EpetraMV(vec); }
374
375 static RCP<const Epetra_CrsMatrix> Op2EpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2EpetraCrs(Op); }
376 static RCP< Epetra_CrsMatrix> Op2NonConstEpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstEpetraCrs(Op); }
377
378 static const Epetra_CrsMatrix& Op2EpetraCrs(const Matrix& Op) { return Utilities::Op2EpetraCrs(Op); }
380
381 static const Epetra_Map& Map2EpetraMap(const Map& map) { return Utilities::Map2EpetraMap(map); }
382 // @}
383#endif
384
385#ifdef HAVE_MUELU_TPETRA
387 // @{
388 static RCP<const Tpetra::MultiVector<SC,LO,GO,NO> > MV2TpetraMV(RCP<MultiVector> const vec) { return Utilities::MV2TpetraMV(vec); }
389 static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV(RCP<MultiVector> vec) { return Utilities::MV2NonConstTpetraMV(vec); }
390 static RCP< Tpetra::MultiVector<SC,LO,GO,NO> > MV2NonConstTpetraMV2(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV2(vec); }
391
392 static const Tpetra::MultiVector<SC,LO,GO,NO>& MV2TpetraMV(const MultiVector& vec) { return Utilities::MV2TpetraMV(vec); }
393 static Tpetra::MultiVector<SC,LO,GO,NO>& MV2NonConstTpetraMV(MultiVector& vec) { return Utilities::MV2NonConstTpetraMV(vec); }
394
395 static RCP<const Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2TpetraCrs(RCP<const Matrix> Op) { return Utilities::Op2TpetraCrs(Op); }
396 static RCP< Tpetra::CrsMatrix<SC,LO,GO,NO> > Op2NonConstTpetraCrs(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
397
398 static const Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2TpetraCrs(const Matrix& Op) { return Utilities::Op2TpetraCrs(Op); }
399 static Tpetra::CrsMatrix<SC,LO,GO,NO>& Op2NonConstTpetraCrs(Matrix& Op) { return Utilities::Op2NonConstTpetraCrs(Op); }
400
401 static RCP<const Tpetra::RowMatrix<SC,LO,GO,NO> > Op2TpetraRow(RCP<const Matrix> Op) { return Utilities::Op2TpetraRow(Op); }
402 static RCP< Tpetra::RowMatrix<SC,LO,GO,NO> > Op2NonConstTpetraRow(RCP<Matrix> Op) { return Utilities::Op2NonConstTpetraRow(Op); }
403
404 static const RCP<const Tpetra::Map<LO, GO, NO> > Map2TpetraMap(const Map& map) { return Utilities::Map2TpetraMap(map); }
405#endif
406 static RCP<Xpetra::Matrix<SC,LO,GO,NO> > Crs2Op(RCP<CrsMatrix> Op) { return Utilities::Crs2Op(Op); }
407
408 static RCP<Vector> GetMatrixDiagonal(const Matrix& A) {
409 const auto rowMap = A.getRowMap();
410 auto diag = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(rowMap,true);
411
412 A.getLocalDiagCopy(*diag);
413
414 return diag;
415 }
416 static RCP<Vector> GetMatrixDiagonalInverse(const Matrix& A, Magnitude tol = Teuchos::ScalarTraits<SC>::eps()*100, const bool doLumped=false);
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, const bool useAverageAbsDiagVal = false) {
419 return UtilitiesBase::GetLumpedMatrixDiagonal(A, doReciprocal, tol, tolReplacement, replaceSingleEntryRowWithZero, useAverageAbsDiagVal);
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
443 static Kokkos::View<bool*, typename Node::device_type> DetectDirichletCols(const Matrix& A, const Kokkos::View<const bool*, typename Node::device_type>& dirichletRows);
444
445 static void FindNonZeros(const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_const_um vals,
446 Kokkos::View<bool*, typename Node::device_type> nonzeros);
447
448 static void DetectDirichletColsAndDomains(const Xpetra::Matrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>& A,
449 const Kokkos::View<bool*, typename Node::device_type> & dirichletRows,
450 Kokkos::View<bool*, typename Node::device_type> dirichletCols,
451 Kokkos::View<bool*, typename Node::device_type> dirichletDomain);
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,
461 Kokkos::View<bool*, typename NO::device_type> & dirichletRows);
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.getLocalMaxNumRowEntries();
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->getLocalNumElements(); ++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->getLocalNumElements(); ++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 >
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;
792 using type = Kokkos::View< DataType, typename view_type::array_layout, typename view_type::device_type, typename CombineMemoryTraits<typename view_type::memory_traits,T>::type >;
793 };
794
795} //namespace MueLu
796
797#define MUELU_UTILITIES_KOKKOS_SHORT
798
799#endif // MUELU_UTILITIES_KOKKOS_DECL_HPP
View
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
int NumMyElements() const
const Epetra_Map & RowMap() const
int ExtractMyRowView(int MyRow, int &NumEntries, double *&Values, int *&Indices) const
Exception throws to report errors in the internal logical of the program.
Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > MultiVector
Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > Map
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
static Teuchos::Array< Magnitude > ResidualNorm(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
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 Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
Xpetra::CrsMatrixWrap< Scalar, LocalOrdinal, GlobalOrdinal, Node > CrsMatrixWrap
Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
static RCP< MultiVector > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const MultiVector &X, const MultiVector &RHS)
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
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< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< const Matrix > Op)
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LO, GO, NO > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
Extract coordinates from parameter list and return them in a Xpetra::MultiVector.
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< MultiVector > vec)
static Teuchos::Array< Magnitude > ResidualNorm(const Operator &Op, const MultiVector &X, const MultiVector &RHS)
static void MyOldScaleMatrix(Matrix &Op, const Teuchos::ArrayRCP< const SC > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
static RCP< Vector > GetInverse(RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< SC >::eps() *100, SC tolReplacement=Teuchos::ScalarTraits< SC >::zero())
static void ZeroDirichletRows(RCP< MultiVector > &X, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static const Tpetra::MultiVector< SC, LO, GO, NO > & MV2TpetraMV(const MultiVector &vec)
static const Epetra_MultiVector & MV2EpetraMV(const MultiVector &vec)
static RCP< Xpetra::Matrix< SC, LO, GO, NO > > Crs2Op(RCP< CrsMatrix > Op)
static RCP< Tpetra::MultiVector< SC, LO, GO, NO > > MV2NonConstTpetraMV2(MultiVector &vec)
Teuchos::ScalarTraits< Scalar >::magnitudeType Magnitude
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Matrix > Op)
Xpetra::MultiVector< CoordinateType, LocalOrdinal, GlobalOrdinal, Node > RealValuedMultiVector
static void ZeroDirichletRows(RCP< Matrix > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static void MyOldScaleMatrix_Epetra(Matrix &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
static Epetra_MultiVector & MV2NonConstEpetraMV(MultiVector &vec)
static void MyOldScaleMatrix_Tpetra(Matrix &Op, const Teuchos::ArrayRCP< Scalar > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
static RCP< const Tpetra::MultiVector< SC, LO, GO, NO > > MV2TpetraMV(RCP< MultiVector > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
static RCP< Tpetra::MultiVector< SC, LO, GO, NO > > MV2NonConstTpetraMV(RCP< MultiVector > vec)
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
static RCP< Tpetra::RowMatrix< SC, LO, GO, NO > > Op2NonConstTpetraRow(RCP< Matrix > Op)
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, const bool useAverageAbsDiagVal=false)
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)
static RCP< const Epetra_MultiVector > MV2EpetraMV(RCP< MultiVector > const vec)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
static RCP< Matrix > Transpose(Matrix &Op, bool=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Transpose a Xpetra::Matrix.
static RCP< MultiVector > Residual(const Operator &Op, const MultiVector &X, const MultiVector &RHS)
static Tpetra::MultiVector< SC, LO, GO, NO > & MV2NonConstTpetraMV(MultiVector &vec)
static Scalar PowerMethod(const Matrix &A, bool scaleByDiag=true, LO niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
static Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraCrs(Matrix &Op)
typename Teuchos::ScalarTraits< Scalar >::coordinateType CoordinateType
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
static const RCP< const Tpetra::Map< LO, GO, NO > > Map2TpetraMap(const Map &map)
static RCP< const Tpetra::RowMatrix< SC, LO, GO, NO > > Op2TpetraRow(RCP< const Matrix > Op)
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
static const Epetra_Map & Map2EpetraMap(const Map &map)
static const Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2TpetraCrs(const Matrix &Op)
static void ApplyOAZToMatrixRows(RCP< Matrix > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows)
static RCP< Tpetra::RowMatrix< SC, LO, GO, NO > > Op2NonConstTpetraRow(RCP< Matrix > Op)
static RCP< const Tpetra::RowMatrix< SC, LO, GO, NO > > Op2TpetraRow(RCP< const Matrix > Op)
static Teuchos::Array< Magnitude > ResidualNorm(const Operator &Op, const MultiVector &X, const MultiVector &RHS)
static RCP< const Tpetra::MultiVector< SC, LO, GO, NO > > MV2TpetraMV(RCP< MultiVector > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
static Epetra_CrsMatrix & Op2NonConstEpetraCrs(Matrix &Op)
static void MyOldScaleMatrix_Tpetra(Matrix &Op, const Teuchos::ArrayRCP< SC > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
static RCP< const Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Matrix > Op)
static RCP< Vector > GetMatrixDiagonal(const Matrix &A)
Extract Matrix Diagonal.
Teuchos::ScalarTraits< SC > TST
Xpetra::MultiVector< CoordinateType, LO, GO, NO > RealValuedMultiVector
static const Tpetra::MultiVector< SC, LO, GO, NO > & MV2TpetraMV(const MultiVector &vec)
static RCP< const Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2TpetraCrs(RCP< const Matrix > Op)
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > ReverseCuthillMcKee(const Matrix &Op)
static const Epetra_CrsMatrix & Op2EpetraCrs(const Matrix &Op)
static SC PowerMethod(const Matrix &A, bool scaleByDiag=true, LO niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
Simple transpose for Tpetra::CrsMatrix types.
static Tpetra::MultiVector< SC, LO, GO, NO > & MV2NonConstTpetraMV(MultiVector &vec)
static RCP< MultiVector > RealValuedToScalarMultiVector(RCP< RealValuedMultiVector > X)
static RCP< Vector > GetMatrixDiagonalInverse(const Matrix &A, Magnitude tol=TST::eps() *100, const bool doLumped=false)
Extract Matrix Diagonal.
static Epetra_MultiVector & MV2NonConstEpetraMV(MultiVector &vec)
static void ZeroDirichletRows(RCP< Matrix > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static void MyOldScaleMatrix_Epetra(Matrix &Op, const Teuchos::ArrayRCP< SC > &scalingVector, bool doFillComplete, bool doOptimizeStorage)
static RCP< MultiVector > Residual(const Operator &Op, const MultiVector &X, const MultiVector &RHS)
static 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)
Detects Dirichlet columns & domains from a list of Dirichlet rows.
static RCP< const Epetra_MultiVector > MV2EpetraMV(RCP< MultiVector > const vec)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
Set seed for random number generator.
static RCP< Xpetra::Vector< LocalOrdinal, LocalOrdinal, GlobalOrdinal, Node > > CuthillMcKee(const Matrix &Op)
static SC PowerMethod(const Matrix &A, const Teuchos::RCP< Vector > &invDiag, LO niters=10, Magnitude tolerance=1e-2, bool verbose=false, unsigned int seed=123)
static RCP< Matrix > Transpose(Matrix &Op, bool optimizeTranspose=false, const std::string &label=std::string())
Transpose a Xpetra::Matrix.
static const RCP< const Tpetra::Map< LO, GO, NO > > Map2TpetraMap(const Map &map)
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
static Kokkos::View< bool *, typename NO::device_type > DetectDirichletCols(const Matrix &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
static RCP< Tpetra::MultiVector< SC, LO, GO, NO > > MV2NonConstTpetraMV2(MultiVector &vec)
static void ZeroDirichletRows(RCP< MultiVector > &X, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< MultiVector > vec)
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Matrix > Op)
static 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)
Find non-zero values in an ArrayRCP Compares the value to 2 * machine epsilon.
static const Epetra_Map & Map2EpetraMap(const Map &map)
static void ApplyRowSumCriterion(const Matrix &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, typename NO::device_type > &dirichletRows)
static void MyOldScaleMatrix(Matrix &Op, const Teuchos::ArrayRCP< const SC > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
static Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
static RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > ExtractCoordinatesFromParameterList(ParameterList &paramList)
static RCP< Xpetra::Matrix< SC, LO, GO, NO > > Crs2Op(RCP< CrsMatrix > Op)
static Tpetra::CrsMatrix< SC, LO, GO, NO > & Op2NonConstTpetraCrs(Matrix &Op)
typename TST::magnitudeType Magnitude
static void ZeroDirichletCols(RCP< Matrix > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletCols, SC replaceWith=Teuchos::ScalarTraits< SC >::zero())
static RCP< Tpetra::MultiVector< SC, LO, GO, NO > > MV2NonConstTpetraMV(RCP< MultiVector > vec)
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
static const Epetra_MultiVector & MV2EpetraMV(const MultiVector &vec)
typename TST::coordinateType CoordinateType
static RCP< Tpetra::CrsMatrix< SC, LO, GO, NO > > Op2NonConstTpetraCrs(RCP< Matrix > Op)
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Teuchos::FancyOStream > MakeFancy(std::ostream &os)
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 Epetra_CrsMatrix > Op2EpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
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 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 const Epetra_Map & Map2EpetraMap(const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &map)
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< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
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_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraCrs(RCP< const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< const Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2TpetraRow(RCP< const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
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.
static void SetRandomSeed(const Teuchos::Comm< int > &comm)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Crs2Op(RCP< Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Namespace for MueLu classes and methods.
KokkosClassic::DefaultNode::DefaultNodeType DefaultNode
magnitude_type tolerance
Kokkos::View< DataType, typename view_type::array_layout, typename view_type::device_type, typename CombineMemoryTraits< typename view_type::memory_traits, T >::type > type