Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Amesos2_Solver_UQ_PCE.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef AMESOS2_SOLVER_UQ_PCE_HPP
43#define AMESOS2_SOLVER_UQ_PCE_HPP
44
45#include "Amesos2_Solver.hpp"
46#include "Amesos2_Factory.hpp"
49
50namespace Amesos2 {
51
52 template <class S, class LO, class GO, class D>
55 const Teuchos::RCP<const Tpetra::CrsMatrix<Sacado::UQ::PCE<S>, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode<D> > >& A = Teuchos::null,
56 const Teuchos::RCP<Tpetra::MultiVector<Sacado::UQ::PCE<S>, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode<D> > >& X = Teuchos::null,
57 const Teuchos::RCP<const Tpetra::MultiVector<Sacado::UQ::PCE<S>, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode<D> > >& B = Teuchos::null)
58 {
59 if (A != Teuchos::null) {
60 return Kokkos::cijk(A->getLocalValuesDevice(Tpetra::Access::ReadOnly));
61 }
62 else if (X != Teuchos::null) {
63 return Kokkos::cijk(X->getLocalViewDevice(Tpetra::Access::ReadOnly));
64 }
65 else if (B != Teuchos::null) {
66 return Kokkos::cijk(B->getLocalViewDevice(Tpetra::Access::ReadOnly));
67 }
68 return typename Sacado::UQ::PCE<S>::cijk_type();
69 }
70
77 template <class Storage, class LocalOrdinal, class GlobalOrdinal,
78 class Device, template<class,class> class ConcreteSolver>
80 public Solver< Tpetra::CrsMatrix<Sacado::UQ::PCE<Storage>,
81 LocalOrdinal,
82 GlobalOrdinal,
83 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >,
84 Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
85 LocalOrdinal,
86 GlobalOrdinal,
87 Kokkos::Compat::KokkosDeviceWrapperNode<Device> >
88 >
89 {
90 public:
91
93 typedef Kokkos::Compat::KokkosDeviceWrapperNode<Device> Node;
94 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Matrix;
95 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Vector;
96
97 typedef typename Scalar::value_type BaseScalar;
98 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Map;
99 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> FlatGraph;
100 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatMatrix;
101 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FlatVector;
102 typedef ConcreteSolver<FlatMatrix,FlatVector> FlatConcreteSolver;
103 typedef Solver<FlatMatrix,FlatVector> FlatSolver;
104
105 typedef Solver<Matrix,Vector> solver_type;
106 typedef typename solver_type::type type;
107 typedef typename Scalar::cijk_type cijk_type;
108
111 const Teuchos::RCP<const Matrix>& A_,
112 const Teuchos::RCP<Vector>& X_,
113 const Teuchos::RCP<const Vector>& B_) :
114 A(A_), X(X_), B(B_) {
115 cijk = get_pce_cijk(A, X, B);
116 const size_t pce_size =
117 Kokkos::dimension_scalar(A->getLocalMatrixDevice().values);
118 flat_graph =
119 Stokhos::create_flat_pce_graph(*(A->getCrsGraph()),
120 cijk,
124 pce_size);
125 if (A != Teuchos::null)
127 if (X != Teuchos::null)
129 if (B != Teuchos::null)
132 create_solver_with_supported_type<ConcreteSolver,FlatMatrix,FlatVector>::apply(flat_A, flat_X, flat_B);
133 }
134
136
137
145 virtual type& preOrdering( void ) {
146 flat_solver->preOrdering();
147 return *this;
148 }
149
150
156 virtual type& symbolicFactorization( void ) {
157 flat_solver->symbolicFactorization();
158 return *this;
159 }
160
161
174 virtual type& numericFactorization( void ) {
175 flat_solver->numericFactorization();
176 return *this;
177 }
178
179
192 virtual void solve( void ) {
193 flat_solver->solve();
194 }
195
196
211 virtual void solve(const Teuchos::Ptr<Vector> XX,
212 const Teuchos::Ptr<const Vector> BB) const {
213 flat_solver->solve(
216 }
217
218
233 virtual void solve(Vector* XX, const Vector* BB) const {
234 flat_solver->solve(
237 }
238
240
241
259 const Teuchos::RCP<Teuchos::ParameterList> & parameterList ) {
260 flat_solver->setParameters(parameterList);
261 return *this;
262 }
263
264
269 virtual Teuchos::RCP<const Teuchos::ParameterList>
270 getValidParameters( void ) const {
271 return flat_solver->getValidParameters();
272 }
273
275
276
300 virtual void setA( const Teuchos::RCP<const Matrix> a,
301 EPhase keep_phase = CLEAN ) {
302 A = a;
303
304 // Rebuild flat matrix/graph
306 if (keep_phase <= CLEAN) {
307 flat_X_map = Teuchos::null;
308 flat_B_map = Teuchos::null;
309 flat_graph = Teuchos::null;
310 flat_graph =
311 Stokhos::create_flat_pce_graph(*(A->getCrsGraph()),
312 cijk,
316 Kokkos::dimension_scalar(A->getLocalMatrixDevice().values));
317 }
318 if (keep_phase <= SYMBFACT) // should this by NUMFACT???
320
321 flat_solver->setA(flat_A, keep_phase);
322 }
323
343 virtual void setA( const Matrix* a, EPhase keep_phase = CLEAN ) {
344 this->setA(Teuchos::rcp(a,false), keep_phase);
345 }
346
347
349 virtual bool matrixShapeOK( void ) {
350 return flat_solver->matrixShapeOK();
351 }
352
353
355 virtual void setX( const Teuchos::RCP<Vector> x ) {
356 X = x;
357 if (x != Teuchos::null)
359 else
360 flat_X = Teuchos::null;
361 flat_solver->setX(flat_X);
362 }
363
364
366 virtual void setX( Vector* x ) {
367 if (x != 0) {
368 X = Teuchos::rcp(x, false);
370 }
371 else {
372 X = Teuchos::null;
373 flat_X = Teuchos::null;
374 }
375 flat_solver->setX(flat_X);
376 }
377
378
380 virtual const Teuchos::RCP<Vector> getX( void ) {
381 return X;
382 }
383
384
386 virtual Vector* getXRaw( void ) {
387 return X.get();
388 }
389
390
392 virtual void setB( const Teuchos::RCP<const Vector> b ) {
393 B = b;
394 if (b != Teuchos::null)
396 else
397 flat_B = Teuchos::null;
398 flat_solver->setB(flat_B);
399 }
400
401
403 virtual void setB( const Vector* b ) {
404 if (b != 0) {
405 B = Teuchos::rcp(b, false);
407 }
408 else {
409 B = Teuchos::null;
410 flat_B = Teuchos::null;
411 }
412 flat_solver->setB(flat_B);
413 }
414
415
417 virtual const Teuchos::RCP<const Vector> getB( void ) {
418 return B;
419 }
420
421
423 virtual const Vector* getBRaw( void ) {
424 return B.get();
425 }
426
427
429 virtual Teuchos::RCP<const Teuchos::Comm<int> > getComm( void ) const {
430 return flat_solver->getComm();
431 }
432
433
435 virtual Status& getStatus() const {
436 return flat_solver->getStatus();
437 }
438
439
441 virtual std::string name( void ) const {
442 return flat_solver->name();
443 }
444
446
447
453 virtual std::string description( void ) const {
454 return flat_solver->description();
455 }
456
457
460 virtual void describe( Teuchos::FancyOStream &out,
461 const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default ) const {
462 flat_solver->describe(out, verbLevel);
463 }
464
466
467
473 virtual void printTiming( Teuchos::FancyOStream &out,
474 const Teuchos::EVerbosityLevel verbLevel = Teuchos::Describable::verbLevel_default ) const{
475 flat_solver->printTiming(out, verbLevel);
476 }
477
478
487 virtual void getTiming( Teuchos::ParameterList& timingParameterList ) const{
488 flat_solver->getTiming(timingParameterList);
489 }
490
492
493 protected:
494
495 Teuchos::RCP<const Matrix> A;
496 Teuchos::RCP<Vector> X;
497 Teuchos::RCP<const Vector> B;
498 Teuchos::RCP<const Map> flat_X_map, flat_B_map;
499 Teuchos::RCP<const FlatGraph> flat_graph, cijk_graph;
500 Teuchos::RCP<const FlatMatrix> flat_A;
501 Teuchos::RCP<FlatVector> flat_X;
502 Teuchos::RCP<const FlatVector> flat_B;
503 Teuchos::RCP<FlatSolver> flat_solver;
505
506 };
507
508 // Specialization of create_solver_with_supported_type for
509 // Sacado::UQ::PCE where we create PCESolverAdapter wrapping
510 // each solver
511 template < template <class,class> class ConcreteSolver,
512 class ST, class LO, class GO, class D >
513 struct create_solver_with_supported_type<
514 ConcreteSolver,
515 Tpetra::CrsMatrix<Sacado::UQ::PCE<ST>,LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<D> >,
516 Tpetra::MultiVector<Sacado::UQ::PCE<ST>,LO,GO,Kokkos::Compat::KokkosDeviceWrapperNode<D> > > {
518 typedef Kokkos::Compat::KokkosDeviceWrapperNode<D> NO;
519 typedef Tpetra::CrsMatrix<SC,LO,GO,NO> Matrix;
520 typedef Tpetra::MultiVector<SC,LO,GO,NO> Vector;
521 static Teuchos::RCP<Solver<Matrix,Vector> >
522 apply(Teuchos::RCP<const Matrix> A,
523 Teuchos::RCP<Vector> X,
524 Teuchos::RCP<const Vector> B )
525 {
526 ctassert<
527 Meta::is_same<
528 typename MatrixTraits<Matrix>::scalar_t,
529 typename MultiVecAdapter<Vector>::scalar_t
530 >::value
531 > same_scalar_assertion;
532 (void)same_scalar_assertion; // This stops the compiler from warning about unused declared variables
533
534 // If our assertion did not fail, then create and return a new solver
535 return Teuchos::rcp( new PCESolverAdapter<ST,LO,GO,D,ConcreteSolver>(A, X, B) );
536 }
537 };
538
539 // Specialization for solver_supports_scalar for Sacado::UQ::PCE<Storage>
540 // value == true if and only if
541 // solver_supprts_scalar<ConcreteSolver,Storage::value_type> == true
542 template <template <class,class> class ConcreteSolver,
543 typename Storage>
544 struct solver_supports_scalar<ConcreteSolver, Sacado::UQ::PCE<Storage> > {
546 typedef typename Scalar::value_type BaseScalar;
547 typedef typename solver_traits<ConcreteSolver>::supported_scalars supported_scalars;
548 static const bool value =
549 Meta::if_then_else<Meta::is_same<supported_scalars, Meta::nil_t>::value,
550 Meta::true_type,
551 Meta::type_list_contains<supported_scalars,
552 BaseScalar> >::type::value;
553 };
554
555
556} // namespace Amesos2
557
558#endif // AMESOS2_SOLVER_UQ_PCE_HPP
Amesos2 solver adapter for UQ::PCE scalar type.
virtual void setX(const Teuchos::RCP< Vector > x)
Sets the LHS vector X.
virtual std::string name(void) const
Return the name of this solver.
virtual Teuchos::RCP< const Teuchos::ParameterList > getValidParameters(void) const
Return a const parameter list of all of the valid parameters that this->setParameterList(....
Tpetra::MultiVector< BaseScalar, LocalOrdinal, GlobalOrdinal, Node > FlatVector
Teuchos::RCP< const Vector > B
Teuchos::RCP< const FlatGraph > cijk_graph
Tpetra::CrsMatrix< BaseScalar, LocalOrdinal, GlobalOrdinal, Node > FlatMatrix
virtual void getTiming(Teuchos::ParameterList &timingParameterList) const
Extracts timing information from the current solver.
virtual type & numericFactorization(void)
Performs numeric factorization on the matrix.
virtual void setB(const Teuchos::RCP< const Vector > b)
Sets the RHS vector B.
Teuchos::RCP< const Map > flat_B_map
virtual Vector * getXRaw(void)
Returns a raw pointer to the LHS of the linear system.
virtual void solve(Vector *XX, const Vector *BB) const
Solve using the given X and B vectors.
virtual bool matrixShapeOK(void)
Returns true if the solver can handle the matrix shape.
Teuchos::RCP< const FlatGraph > flat_graph
virtual void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Teuchos::RCP< const FlatMatrix > flat_A
Teuchos::RCP< const Matrix > A
ConcreteSolver< FlatMatrix, FlatVector > FlatConcreteSolver
Teuchos::RCP< const FlatVector > flat_B
virtual type & setParameters(const Teuchos::RCP< Teuchos::ParameterList > &parameterList)
Set/update internal variables and solver options.
virtual void solve(void)
Solves (or )
virtual type & symbolicFactorization(void)
Performs symbolic factorization on the matrix.
virtual const Vector * getBRaw(void)
Returns a raw pointer to the RHS of the linear system.
virtual void printTiming(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Prints timing information about the current solver.
virtual Status & getStatus() const
Returns a reference to this solver's internal status object.
virtual type & preOrdering(void)
Pre-orders the matrix.
Teuchos::RCP< FlatVector > flat_X
virtual void setA(const Matrix *a, EPhase keep_phase=CLEAN)
Sets the matrix A of this solver.
Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > Matrix
Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > Map
PCESolverAdapter(const Teuchos::RCP< const Matrix > &A_, const Teuchos::RCP< Vector > &X_, const Teuchos::RCP< const Vector > &B_)
Constructor.
virtual void solve(const Teuchos::Ptr< Vector > XX, const Teuchos::Ptr< const Vector > BB) const
Solve using the given X and B vectors.
Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > FlatGraph
Sacado::UQ::PCE< Storage > Scalar
virtual std::string description(void) const
Returns a short description of this Solver.
virtual void setX(Vector *x)
Sets the LHS vector X using a raw pointer.
virtual Teuchos::RCP< const Teuchos::Comm< int > > getComm(void) const
Returns a pointer to the Teuchos::Comm communicator with this matrix.
Teuchos::RCP< const Map > flat_X_map
virtual void setA(const Teuchos::RCP< const Matrix > a, EPhase keep_phase=CLEAN)
Sets the matrix A of this solver.
virtual const Teuchos::RCP< const Vector > getB(void)
Returns the vector that is the RHS of the linear system.
Solver< FlatMatrix, FlatVector > FlatSolver
virtual const Teuchos::RCP< Vector > getX(void)
Returns the vector that is the LHS of the linear system.
virtual void setB(const Vector *b)
Sets the RHS vector B using a raw pointer.
Solver< Matrix, Vector > solver_type
Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > Vector
Kokkos::Compat::KokkosDeviceWrapperNode< Device > Node
Teuchos::RCP< FlatSolver > flat_solver
Stokhos::StandardStorage< int, double > Storage
Sacado::UQ::PCE< S >::cijk_type get_pce_cijk(const Teuchos::RCP< const Tpetra::CrsMatrix< Sacado::UQ::PCE< S >, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< D > > > &A=Teuchos::null, const Teuchos::RCP< Tpetra::MultiVector< Sacado::UQ::PCE< S >, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< D > > > &X=Teuchos::null, const Teuchos::RCP< const Tpetra::MultiVector< Sacado::UQ::PCE< S >, LO, GO, Kokkos::Compat::KokkosDeviceWrapperNode< D > > > &B=Teuchos::null)
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< View< T, P... > >::value, unsigned >::type dimension_scalar(const View< T, P... > &view)
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typenameCijkType< view_type >::type >::type cijk(const view_type &view)
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_pce_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &graph, const CijkType &cijk, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_range_map, Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &cijk_graph, const size_t matrix_pce_size)
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_map)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &cijk_graph, const CijkType &cijk_dev)