10 #ifndef BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
11 #define BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
13 #include "BelosTpetraAdapter.hpp"
15 #include "KokkosBlas.hpp"
17 #ifdef HAVE_BELOS_TSQR
19 #endif // HAVE_BELOS_TSQR
39 template<
class Storage,
class LO,
class GO,
class Node>
41 Tpetra::MultiVector< Sacado::UQ::PCE<Storage>,
47 typedef Tpetra::MultiVector<Scalar, LO, GO, Node>
MV;
49 typedef typename Tpetra::MultiVector<Scalar,LO,GO,Node>::dot_type
dot_type;
50 typedef typename Tpetra::MultiVector<Scalar,LO,GO,Node>::mag_type
mag_type;
94 #ifdef HAVE_TPETRA_DEBUG
95 const char fnName[] =
"Belos::MultiVecTraits::CloneCopy(mv,index)";
96 const size_t inNumVecs = mv.getNumVectors ();
98 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
99 std::runtime_error, fnName <<
": All indices must be nonnegative.");
102 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
104 fnName <<
": All indices must be strictly less than the number of "
105 "columns " << inNumVecs <<
" of the input multivector mv.");
106 #endif // HAVE_TPETRA_DEBUG
110 for (std::vector<int>::size_type
j = 0;
j < index.size (); ++
j) {
111 columns[
j] = index[
j];
130 const bool validRange = index.
size() > 0 &&
132 index.
ubound() < GetNumberVecs(mv);
134 std::ostringstream os;
135 os <<
"Belos::MultiVecTraits::CloneCopy(mv,index=["
138 index.
size() == 0, std::invalid_argument,
139 os.str() <<
"Empty index range is not allowed.");
141 index.
lbound() < 0, std::invalid_argument,
142 os.str() <<
"Index range includes negative index/ices, which is not "
145 index.
ubound() >= GetNumberVecs(mv), std::invalid_argument,
146 os.str() <<
"Index range exceeds number of vectors "
147 << mv.getNumVectors() <<
" in the input multivector.");
149 os.str() <<
"Should never get here!");
160 #ifdef HAVE_TPETRA_DEBUG
161 const char fnName[] =
"Belos::MultiVecTraits::CloneViewNonConst(mv,index)";
162 const size_t numVecs = mv.getNumVectors ();
164 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
165 std::invalid_argument,
166 fnName <<
": All indices must be nonnegative.");
169 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
170 std::invalid_argument,
171 fnName <<
": All indices must be strictly less than the number of "
172 "columns " << numVecs <<
" in the input MultiVector mv.");
173 #endif // HAVE_TPETRA_DEBUG
177 for (std::vector<int>::size_type
j = 0;
j < index.size (); ++
j) {
178 columns[
j] = index[
j];
194 const int numCols =
static_cast<int> (mv.getNumVectors());
195 const bool validRange = index.
size() > 0 &&
198 std::ostringstream os;
199 os <<
"Belos::MultiVecTraits::CloneViewNonConst(mv,index=["
202 index.
size() == 0, std::invalid_argument,
203 os.str() <<
"Empty index range is not allowed.");
205 index.
lbound() < 0, std::invalid_argument,
206 os.str() <<
"Index range includes negative inde{x,ices}, which is "
209 index.
ubound() >= numCols, std::invalid_argument,
210 os.str() <<
"Index range exceeds number of vectors " << numCols
211 <<
" in the input multivector.");
213 true, std::logic_error,
214 os.str() <<
"Should never get here!");
224 #ifdef HAVE_TPETRA_DEBUG
225 const char fnName[] =
"Belos::MultiVecTraits<Scalar, "
226 "Tpetra::MultiVector<...> >::CloneView(mv,index)";
227 const size_t numVecs = mv.getNumVectors ();
229 *std::min_element (index.begin (), index.end ()) < 0,
230 std::invalid_argument,
231 fnName <<
": All indices must be nonnegative.");
233 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
234 std::invalid_argument,
235 fnName <<
": All indices must be strictly less than the number of "
236 "columns " << numVecs <<
" in the input MultiVector mv.");
237 #endif // HAVE_TPETRA_DEBUG
241 for (std::vector<int>::size_type
j = 0;
j < index.size (); ++
j) {
242 columns[
j] = index[
j];
258 const int numCols =
static_cast<int> (mv.getNumVectors());
259 const bool validRange = index.
size() > 0 &&
262 std::ostringstream os;
263 os <<
"Belos::MultiVecTraits::CloneView(mv, index=["
266 os.str() <<
"Empty index range is not allowed.");
268 os.str() <<
"Index range includes negative index/ices, which is not "
271 index.
ubound() >= numCols, std::invalid_argument,
272 os.str() <<
"Index range exceeds number of vectors " << numCols
273 <<
" in the input multivector.");
275 os.str() <<
"Should never get here!");
283 return static_cast<ptrdiff_t
> (mv.getGlobalLength ());
287 return static_cast<int> (mv.getNumVectors ());
291 return mv.isConstantStride ();
296 const Tpetra::MultiVector<Scalar,LO,GO,Node>&
A,
299 Tpetra::MultiVector<Scalar,LO,GO,Node>&
C)
305 const int numRowsB = B.
numRows ();
306 const int numColsB = B.
numCols ();
307 const int strideB = B.
stride ();
308 if (numRowsB == 1 && numColsB == 1) {
309 C.update (alpha*B(0,0), A, beta);
317 else Atmp =
rcp(&A,
false);
320 else Ctmp =
rcp(&C,
false);
323 typedef Tpetra::MultiVector<dot_type,LO,GO,Node> FMV;
324 typedef Tpetra::MultiVector<const dot_type,LO,GO,Node> CFMV;
325 typedef typename FMV::dual_view_type::t_dev flat_view_type;
326 typedef typename CFMV::dual_view_type::t_dev const_flat_view_type;
328 const_flat_view_type flat_A_view = Atmp->template getLocalView<execution_space>(Tpetra::Access::ReadOnly);
329 flat_view_type flat_C_view = Ctmp->template getLocalView<execution_space>(Tpetra::Access::ReadWrite);
332 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> b_host_view_type;
333 b_host_view_type B_view_host_input( B.
values(), strideB, numColsB);
334 auto B_view_host = Kokkos::subview( B_view_host_input,
335 Kokkos::pair<int,int>(0, numRowsB),
336 Kokkos::pair<int,int>(0, numColsB));
340 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> b_view_type;
341 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> b_1d_view_type;
342 b_1d_view_type B_1d_view_dev(Kokkos::ViewAllocateWithoutInitializing(
"B"), numRowsB*numColsB);
343 b_view_type B_view_dev( B_1d_view_dev.data(), numRowsB, numColsB);
347 for (
int j=0;
j<numColsB; ++
j)
349 Kokkos::subview(B_view_host,Kokkos::ALL,
j));
353 const char ctransA =
'N', ctransB =
'N';
356 alpha, flat_A_view, B_view_dev, beta, flat_C_view);
359 if (C.isConstantStride() ==
false)
372 const Tpetra::MultiVector<Scalar,LO,GO,Node>&
A,
374 const Tpetra::MultiVector<Scalar,LO,GO,Node>&
B,
375 Tpetra::MultiVector<Scalar,LO,GO,Node>& mv)
381 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
388 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
389 const std::vector<BaseScalar>& alphas)
391 std::vector<Scalar> alphas_mp(alphas.size());
392 const size_t sz = alphas.size();
393 for (
size_t i=0; i<sz; ++i)
394 alphas_mp[i] = alphas[i];
395 mv.scale (alphas_mp);
399 MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
400 const std::vector<Scalar>& alphas)
407 const Tpetra::MultiVector<Scalar,LO,GO,Node>&
A,
408 const Tpetra::MultiVector<Scalar,LO,GO,Node>&
B,
415 using Teuchos::reduceAll;
418 const int numRowsC = C.
numRows ();
419 const int numColsC = C.
numCols ();
420 const int strideC = C.
stride ();
421 if (numRowsC == 1 && numColsC == 1) {
435 RCP<const MV> Atmp, Btmp;
437 else Atmp =
rcp(&A,
false);
440 else Btmp =
rcp(&B,
false);
443 typedef Tpetra::MultiVector<const dot_type,LO,GO,Node> FMV;
444 typedef typename FMV::dual_view_type::t_dev flat_view_type;
446 flat_view_type flat_A_view = Atmp->template getLocalView<execution_space>(Tpetra::Access::ReadOnly);
447 flat_view_type flat_B_view = Btmp->template getLocalView<execution_space>(Tpetra::Access::ReadOnly);
450 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> c_host_view_type;
451 c_host_view_type C_view_host_input( C.
values(), strideC, numColsC);
452 auto C_view_host = Kokkos::subview(C_view_host_input,
453 Kokkos::pair<int,int>(0, numRowsC),
454 Kokkos::pair<int,int>(0, numColsC));
458 typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> c_view_type;
459 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> c_1d_view_type;
460 c_1d_view_type C_1d_view_dev(
"C", numRowsC*numColsC);
461 c_view_type C_view_dev( C_1d_view_dev.data(), numRowsC, numColsC);
465 const char ctransA =
'C', ctransB =
'N';
468 alpha, flat_A_view, flat_B_view,
469 #
if KOKKOS_VERSION >= 40799
470 KokkosKernels::ArithTraits<dot_type>::zero(),
472 Kokkos::ArithTraits<dot_type>::zero(),
477 RCP<const Comm<int> > pcomm = A.getMap()->getComm ();
478 if (pcomm->getSize () == 1) {
482 for (
int j=0;
j<numColsC; ++
j)
484 Kokkos::subview(C_view_dev,Kokkos::ALL,
j));
487 typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, Kokkos::HostSpace> c_1d_host_view_type;
488 c_1d_host_view_type C_1d_view_tmp(Kokkos::ViewAllocateWithoutInitializing(
"C_tmp"), strideC*numColsC);
489 c_host_view_type C_view_tmp( C_1d_view_tmp.data(),
491 for (
int j=0;
j<numColsC; ++
j)
493 Kokkos::pair<int,int>(0, numRowsC),
495 Kokkos::subview(C_view_dev, Kokkos::ALL,
j));
496 reduceAll<int> (*pcomm, REDUCE_SUM, strideC*numColsC,
504 MvDot (
const Tpetra::MultiVector<Scalar,LO,GO,Node>&
A,
505 const Tpetra::MultiVector<Scalar,LO,GO,Node>&
B,
506 std::vector<dot_type>& dots)
508 const size_t numVecs = A.getNumVectors ();
511 numVecs != B.getNumVectors (), std::invalid_argument,
512 "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): "
513 "A and B must have the same number of columns. "
514 "A has " << numVecs <<
" column(s), "
515 "but B has " << B.getNumVectors () <<
" column(s).");
516 #ifdef HAVE_TPETRA_DEBUG
518 dots.size() < numVecs, std::invalid_argument,
519 "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): "
520 "The output array 'dots' must have room for all dot products. "
521 "A and B each have " << numVecs <<
" column(s), "
522 "but 'dots' only has " << dots.size() <<
" entry(/ies).");
523 #endif // HAVE_TPETRA_DEBUG
526 A.dot (B, av (0, numVecs));
531 MvNorm (
const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
532 std::vector<mag_type> &normvec,
533 NormType type=TwoNorm)
536 #ifdef HAVE_TPETRA_DEBUG
537 typedef std::vector<int>::size_type size_type;
539 normvec.size () <
static_cast<size_type
> (mv.getNumVectors ()),
540 std::invalid_argument,
541 "Belos::MultiVecTraits::MvNorm(mv,normvec): The normvec output "
542 "argument must have at least as many entries as the number of vectors "
543 "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
544 <<
" < mv.getNumVectors() = " << mv.getNumVectors () <<
".");
550 mv.norm1(av(0,mv.getNumVectors()));
553 mv.norm2(av(0,mv.getNumVectors()));
556 mv.normInf(av(0,mv.getNumVectors()));
563 true, std::logic_error,
564 "Belos::MultiVecTraits::MvNorm: Invalid NormType value " << type
565 <<
". Valid values are OneNorm=" << OneNorm <<
", TwoNorm="
566 << TwoNorm <<
", and InfNorm=" << InfNorm <<
". If you are a Belos "
567 "user and have not modified Belos in any way, and you get this "
568 "message, then this is probably a bug in the Belos solver you were "
569 "using. Please report this to the Belos developers.");
578 const size_t inNumVecs = A.getNumVectors ();
579 #ifdef HAVE_TPETRA_DEBUG
581 inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
582 "Belos::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must "
583 "have no more entries as the number of columns in the input MultiVector"
584 " A. A.getNumVectors() = " << inNumVecs <<
" < index.size () = "
585 << index.size () <<
".");
586 #endif // HAVE_TPETRA_DEBUG
587 RCP<MV> mvsub = CloneViewNonConst (mv, index);
588 if (inNumVecs > static_cast<size_t> (index.size ())) {
589 RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
606 const size_t maxInt =
608 const bool overflow =
609 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
611 std::ostringstream os;
612 os <<
"Belos::MultiVecTraits::SetBlock(A, index=[" << index.
lbound ()
613 <<
", " << index.
ubound () <<
"], mv): ";
615 maxInt < A.getNumVectors (), std::range_error, os.str () <<
"Number "
616 "of columns (size_t) in the input MultiVector 'A' overflows int.");
618 maxInt < mv.getNumVectors (), std::range_error, os.str () <<
"Number "
619 "of columns (size_t) in the output MultiVector 'mv' overflows int.");
622 const int numColsA =
static_cast<int> (A.getNumVectors ());
623 const int numColsMv =
static_cast<int> (mv.getNumVectors ());
625 const bool validIndex =
628 const bool validSource = index.
size () <= numColsA;
630 if (! validIndex || ! validSource) {
631 std::ostringstream os;
632 os <<
"Belos::MultiVecTraits::SetBlock(A, index=[" << index.
lbound ()
633 <<
", " << index.
ubound () <<
"], mv): ";
635 index.
lbound() < 0, std::invalid_argument,
636 os.str() <<
"Range lower bound must be nonnegative.");
638 index.
ubound() >= numColsMv, std::invalid_argument,
639 os.str() <<
"Range upper bound must be less than the number of "
640 "columns " << numColsA <<
" in the 'mv' output argument.");
642 index.
size() > numColsA, std::invalid_argument,
643 os.str() <<
"Range must have no more elements than the number of "
644 "columns " << numColsA <<
" in the 'A' input argument.");
646 true, std::logic_error,
"Should never get here!");
653 if (index.
lbound () == 0 && index.
ubound () + 1 == numColsMv) {
654 mv_view = Teuchos::rcpFromRef (mv);
656 mv_view = CloneViewNonConst (mv, index);
663 if (index.
size () == numColsA) {
664 A_view = Teuchos::rcpFromRef (A);
674 const char errPrefix[] =
"Belos::MultiVecTraits::Assign(A, mv): ";
683 const size_t maxInt =
685 const bool overflow =
686 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
689 maxInt < A.getNumVectors(), std::range_error,
690 errPrefix <<
"Number of columns in the input multivector 'A' "
691 "(a size_t) overflows int.");
693 maxInt < mv.getNumVectors(), std::range_error,
694 errPrefix <<
"Number of columns in the output multivector 'mv' "
695 "(a size_t) overflows int.");
697 true, std::logic_error,
"Should never get here!");
700 const int numColsA =
static_cast<int> (A.getNumVectors ());
701 const int numColsMv =
static_cast<int> (mv.getNumVectors ());
702 if (numColsA > numColsMv) {
704 numColsA > numColsMv, std::invalid_argument,
705 errPrefix <<
"Input multivector 'A' has " << numColsA <<
" columns, "
706 "but output multivector 'mv' has only " << numColsMv <<
" columns.");
709 if (numColsA == numColsMv) {
726 mv.putScalar (alpha);
734 #ifdef HAVE_BELOS_TSQR
735 typedef Tpetra::TsqrAdaptor< Tpetra::MultiVector< Scalar, LO, GO, Node > > tsqr_adaptor_type;
739 #endif // HAVE_BELOS_TSQR
749 template <
class Storage,
class LO,
class GO,
class Node>
751 Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
753 Tpetra::Operator<Sacado::UQ::PCE<Storage>,
759 Apply (
const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
760 const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
761 Tpetra::MultiVector<Scalar,LO,GO,Node>& Y,
762 ETrans trans=NOTRANS)
779 const std::string otName =
"Belos::OperatorTraits<" + scalarName
780 +
"," + loName +
"," + goName +
"," + nodeName +
">";
782 "get here; fell through a switch statement. "
783 "Please report this bug to the Belos developers.");
790 return Op.hasTransposeApply ();
ScalarType * values() const
Tpetra::MultiVector< Scalar, LO, GO, Node > MV
static Teuchos::RCP< MV > CloneCopy(const MV &X)
Create and return a deep copy of X.
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
static bool HasConstantStride(const MV &mv)
Sacado::UQ::PCE< Storage > Scalar
Tpetra::MultiVector< Scalar, LO, GO, Node >::mag_type mag_type
Kokkos::DefaultExecutionSpace execution_space
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void MvPrint(const MV &mv, std::ostream &os)
static Teuchos::RCP< const MV > CloneView(const MV &mv, const Teuchos::Range1D &index)
static void MvTransMv(dot_type alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, Teuchos::SerialDenseMatrix< int, dot_type > &C)
Tpetra::MultiVector< Scalar, LO, GO, Node >::dot_type dot_type
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< Scalar > &alphas)
static bool HasApplyTranspose(const Tpetra::Operator< Scalar, LO, GO, Node > &Op)
Sacado::UQ::PCE< Storage > Scalar
static void SetBlock(const MV &A, const Teuchos::Range1D &index, MV &mv)
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const Scalar &alpha)
static void Apply(const Tpetra::Operator< Scalar, LO, GO, Node > &Op, const Tpetra::MultiVector< Scalar, LO, GO, Node > &X, Tpetra::MultiVector< Scalar, LO, GO, Node > &Y, ETrans trans=NOTRANS)
static void MvDot(const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, std::vector< dot_type > &dots)
For all columns j of A, set dots[j] := A[j]^T * B[j].
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const Teuchos::Range1D &index)
Create and return a deep copy of the given columns of mv.
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
static void MvAddMv(Scalar alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, Scalar beta, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, Tpetra::MultiVector< Scalar, LO, GO, Node > &mv)
mv := alpha*A + beta*B
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< BaseScalar > &alphas)
OrdinalType numCols() const
static void MvRandom(MV &mv)
static int GetNumberVecs(const MV &mv)
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const std::vector< int > &index)
Create and return a deep copy of the given columns of mv.
std::enable_if< Kokkos::is_view_mp_vector< Kokkos::View< DA, PA...> >::value &&Kokkos::is_view_mp_vector< Kokkos::View< DB, PB...> >::value &&Kokkos::is_view_mp_vector< Kokkos::View< DC, PC...> >::value >::type gemm(const char transA[], const char transB[], typename Kokkos::View< DA, PA...>::const_value_type &alpha, const Kokkos::View< DA, PA...> &A, const Kokkos::View< DB, PB...> &B, typename Kokkos::View< DC, PC...>::const_value_type &beta, const Kokkos::View< DC, PC...> &C)
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
static Teuchos::RCP< MV > Clone(const MV &X, const int numVecs)
Create a new MultiVector with numVecs columns.
static void MvNorm(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, std::vector< mag_type > &normvec, NormType type=TwoNorm)
For all columns j of mv, set normvec[j] = norm(mv[j]).
static void Assign(const MV &A, MV &mv)
Storage::value_type BaseScalar
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const Teuchos::Range1D &index)
static ptrdiff_t GetGlobalLength(const MV &mv)
OrdinalType stride() const
Storage::ordinal_type s_ordinal
static std::string name()
OrdinalType numRows() const
static void MvTimesMatAddMv(const dot_type &alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Teuchos::SerialDenseMatrix< int, dot_type > &B, const dot_type &beta, Tpetra::MultiVector< Scalar, LO, GO, Node > &C)