Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Belos_TpetraAdapter_UQ_PCE.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
11 #define BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
12 
13 #include "BelosTpetraAdapter.hpp"
15 #include "KokkosBlas.hpp"
16 
17 #ifdef HAVE_BELOS_TSQR
19 #endif // HAVE_BELOS_TSQR
20 
21 namespace Belos {
22 
24  //
25  // Implementation of Belos::MultiVecTraits for Tpetra::MultiVector.
26  //
28 
39  template<class Storage, class LO, class GO, class Node>
40  class MultiVecTraits<typename Storage::value_type,
41  Tpetra::MultiVector< Sacado::UQ::PCE<Storage>,
42  LO, GO, Node > > {
43  public:
44  typedef typename Storage::ordinal_type s_ordinal;
45  typedef typename Storage::value_type BaseScalar;
47  typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
48  public:
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;
51 
57  static Teuchos::RCP<MV> Clone (const MV& X, const int numVecs) {
58  Teuchos::RCP<MV> Y (new MV (X.getMap (), numVecs, false));
59  Y->setCopyOrView (Teuchos::View);
60  return Y;
61  }
62 
64  static Teuchos::RCP<MV> CloneCopy (const MV& X)
65  {
66  // Make a deep copy of X. The one-argument copy constructor
67  // does a shallow copy by default; the second argument tells it
68  // to do a deep copy.
69  Teuchos::RCP<MV> X_copy (new MV (X, Teuchos::Copy));
70  // Make Tpetra::MultiVector use the new view semantics. This is
71  // a no-op for the Kokkos refactor version of Tpetra; it only
72  // does something for the "classic" version of Tpetra. This
73  // shouldn't matter because Belos only handles MV through RCP
74  // and through this interface anyway, but it doesn't hurt to set
75  // it and make sure that it works.
76  X_copy->setCopyOrView (Teuchos::View);
77  return X_copy;
78  }
79 
91  static Teuchos::RCP<MV>
92  CloneCopy (const MV& mv, const std::vector<int>& index)
93  {
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.");
101  index.size () > 0 &&
102  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
103  std::runtime_error,
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
107 
108  // Tpetra wants an array of size_t, not of int.
109  Teuchos::Array<size_t> columns (index.size ());
110  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
111  columns[j] = index[j];
112  }
113  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
114  // continuous column index range in MultiVector::subCopy, so we
115  // don't have to check here.
116  Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
117  X_copy->setCopyOrView (Teuchos::View);
118  return X_copy;
119  }
120 
127  static Teuchos::RCP<MV>
128  CloneCopy (const MV& mv, const Teuchos::Range1D& index)
129  {
130  const bool validRange = index.size() > 0 &&
131  index.lbound() >= 0 &&
132  index.ubound() < GetNumberVecs(mv);
133  if (! validRange) { // invalid range; generate error message
134  std::ostringstream os;
135  os << "Belos::MultiVecTraits::CloneCopy(mv,index=["
136  << index.lbound() << "," << index.ubound() << "]): ";
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 "
143  "allowed.");
145  index.ubound() >= GetNumberVecs(mv), std::invalid_argument,
146  os.str() << "Index range exceeds number of vectors "
147  << mv.getNumVectors() << " in the input multivector.");
148  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
149  os.str() << "Should never get here!");
150  }
151  Teuchos::RCP<MV> X_copy = mv.subCopy (index);
152  X_copy->setCopyOrView (Teuchos::View);
153  return X_copy;
154  }
155 
156 
157  static Teuchos::RCP<MV>
158  CloneViewNonConst (MV& mv, const std::vector<int>& index)
159  {
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.");
168  index.size () > 0 &&
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
174 
175  // Tpetra wants an array of size_t, not of int.
176  Teuchos::Array<size_t> columns (index.size ());
177  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
178  columns[j] = index[j];
179  }
180  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
181  // continuous column index range in
182  // MultiVector::subViewNonConst, so we don't have to check here.
183  Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
184  X_view->setCopyOrView (Teuchos::View);
185  return X_view;
186  }
187 
188  static Teuchos::RCP<MV>
190  {
191  // NOTE (mfh 11 Jan 2011) We really should check for possible
192  // overflow of int here. However, the number of columns in a
193  // multivector typically fits in an int.
194  const int numCols = static_cast<int> (mv.getNumVectors());
195  const bool validRange = index.size() > 0 &&
196  index.lbound() >= 0 && index.ubound() < numCols;
197  if (! validRange) {
198  std::ostringstream os;
199  os << "Belos::MultiVecTraits::CloneViewNonConst(mv,index=["
200  << index.lbound() << ", " << index.ubound() << "]): ";
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 "
207  "not allowed.");
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!");
215  }
216  Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
217  X_view->setCopyOrView (Teuchos::View);
218  return X_view;
219  }
220 
222  CloneView (const MV& mv, const std::vector<int>& index)
223  {
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
238 
239  // Tpetra wants an array of size_t, not of int.
240  Teuchos::Array<size_t> columns (index.size ());
241  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
242  columns[j] = index[j];
243  }
244  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
245  // continuous column index range in MultiVector::subView, so we
246  // don't have to check here.
247  Teuchos::RCP<const MV> X_view = mv.subView (columns);
248  Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
249  return X_view;
250  }
251 
253  CloneView (const MV& mv, const Teuchos::Range1D& index)
254  {
255  // NOTE (mfh 11 Jan 2011) We really should check for possible
256  // overflow of int here. However, the number of columns in a
257  // multivector typically fits in an int.
258  const int numCols = static_cast<int> (mv.getNumVectors());
259  const bool validRange = index.size() > 0 &&
260  index.lbound() >= 0 && index.ubound() < numCols;
261  if (! validRange) {
262  std::ostringstream os;
263  os << "Belos::MultiVecTraits::CloneView(mv, index=["
264  << index.lbound () << ", " << index.ubound() << "]): ";
265  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
266  os.str() << "Empty index range is not allowed.");
267  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
268  os.str() << "Index range includes negative index/ices, which is not "
269  "allowed.");
271  index.ubound() >= numCols, std::invalid_argument,
272  os.str() << "Index range exceeds number of vectors " << numCols
273  << " in the input multivector.");
274  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
275  os.str() << "Should never get here!");
276  }
277  Teuchos::RCP<const MV> X_view = mv.subView (index);
278  Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
279  return X_view;
280  }
281 
282  static ptrdiff_t GetGlobalLength (const MV& mv) {
283  return static_cast<ptrdiff_t> (mv.getGlobalLength ());
284  }
285 
286  static int GetNumberVecs (const MV& mv) {
287  return static_cast<int> (mv.getNumVectors ());
288  }
289 
290  static bool HasConstantStride (const MV& mv) {
291  return mv.isConstantStride ();
292  }
293 
294  static void
295  MvTimesMatAddMv (const dot_type& alpha,
296  const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
298  const dot_type& beta,
299  Tpetra::MultiVector<Scalar,LO,GO,Node>& C)
300  {
301  using Teuchos::RCP;
302  using Teuchos::rcp;
303 
304  // Check if numRowsB == numColsB == 1, in which case we can call update()
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);
310  return;
311  }
312 
313  // Ensure A and C have constant stride
314  RCP<const MV> Atmp;
315  RCP<MV> Ctmp;
316  if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
317  else Atmp = rcp(&A,false);
318 
319  if (C.isConstantStride() == false) Ctmp = rcp (new MV (C, Teuchos::Copy));
320  else Ctmp = rcp(&C,false);
321 
322  // Create flattened view's
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);
330 
331  // Create a view for B on the host
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));
337 
338  // Create view for B on the device -- need to be careful to get the
339  // right stride to match B
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);
344  //Kokkos::deep_copy(B_view_dev, B_view_host);
345  // Device-to-host copies requires dest to be contiguous, which
346  // C_view_host may not be. So do 1 column at a time.
347  for (int j=0; j<numColsB; ++j)
348  Kokkos::deep_copy(Kokkos::subview(B_view_dev,Kokkos::ALL,j),
349  Kokkos::subview(B_view_host,Kokkos::ALL,j));
350 
351  // Do local multiply
352  {
353  const char ctransA = 'N', ctransB = 'N';
355  &ctransA, &ctransB,
356  alpha, flat_A_view, B_view_dev, beta, flat_C_view);
357  }
358  // Copy back to C if we made a copy
359  if (C.isConstantStride() == false)
360  C.assign(*Ctmp);
361  }
362 
370  static void
371  MvAddMv (Scalar alpha,
372  const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
373  Scalar beta,
374  const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
375  Tpetra::MultiVector<Scalar,LO,GO,Node>& mv)
376  {
377  mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
378  }
379 
380  static void
381  MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
382  const Scalar& alpha)
383  {
384  mv.scale (alpha);
385  }
386 
387  static void
388  MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
389  const std::vector<BaseScalar>& alphas)
390  {
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);
396  }
397 
398  static void
399  MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
400  const std::vector<Scalar>& alphas)
401  {
402  mv.scale (alphas);
403  }
404 
405  static void
407  const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
408  const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
410  {
411  using Teuchos::Comm;
412  using Teuchos::RCP;
413  using Teuchos::rcp;
414  using Teuchos::REDUCE_SUM;
415  using Teuchos::reduceAll;
416 
417  // Check if numRowsC == numColsC == 1, in which case we can call dot()
418  const int numRowsC = C.numRows ();
419  const int numColsC = C.numCols ();
420  const int strideC = C.stride ();
421  if (numRowsC == 1 && numColsC == 1) {
422  if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
423  // Short-circuit, as required by BLAS semantics.
424  C(0,0) = alpha;
425  return;
426  }
427  A.dot (B, Teuchos::ArrayView<dot_type> (C.values (), 1));
428  if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
429  C(0,0) *= alpha;
430  }
431  return;
432  }
433 
434  // Ensure A and B have constant stride
435  RCP<const MV> Atmp, Btmp;
436  if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
437  else Atmp = rcp(&A,false);
438 
439  if (B.isConstantStride() == false) Btmp = rcp (new MV (B, Teuchos::Copy));
440  else Btmp = rcp(&B,false);
441 
442  // Create flattened Kokkos::MultiVector's
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);
448 
449  // Create a view for C on the host
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));
455 
456  // Create view for C on the device -- need to be careful to get the
457  // right stride to match C (allow setting to 0 for first-touch)
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);
462 
463  // Do local multiply
464  {
465  const char ctransA = 'C', ctransB = 'N';
467  &ctransA, &ctransB,
468  alpha, flat_A_view, flat_B_view,
469 #if KOKKOS_VERSION >= 40799
470  KokkosKernels::ArithTraits<dot_type>::zero(),
471 #else
472  Kokkos::ArithTraits<dot_type>::zero(),
473 #endif
474  C_view_dev);
475  }
476  // reduce across processors -- could check for RDMA
477  RCP<const Comm<int> > pcomm = A.getMap()->getComm ();
478  if (pcomm->getSize () == 1) {
479  //Kokkos::deep_copy(C_view_host, C_view_dev);
480  // Device-to-host copies requires dest to be contiguous, which
481  // C_view_host may not be. So do 1 column at a time.
482  for (int j=0; j<numColsC; ++j)
483  Kokkos::deep_copy(Kokkos::subview(C_view_host,Kokkos::ALL,j),
484  Kokkos::subview(C_view_dev,Kokkos::ALL,j));
485  }
486  else {
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(),
490  strideC, numColsC);
491  for (int j=0; j<numColsC; ++j)
492  Kokkos::deep_copy(Kokkos::subview(C_view_tmp,
493  Kokkos::pair<int,int>(0, numRowsC),
494  j),
495  Kokkos::subview(C_view_dev, Kokkos::ALL, j));
496  reduceAll<int> (*pcomm, REDUCE_SUM, strideC*numColsC,
497  C_view_tmp.data(),
498  C_view_host.data());
499  }
500  }
501 
503  static void
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)
507  {
508  const size_t numVecs = A.getNumVectors ();
509 
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
524 
526  A.dot (B, av (0, numVecs));
527  }
528 
530  static void
531  MvNorm (const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
532  std::vector<mag_type> &normvec,
533  NormType type=TwoNorm)
534  {
535 
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 () << ".");
545 #endif
546 
547  Teuchos::ArrayView<mag_type> av(normvec);
548  switch (type) {
549  case OneNorm:
550  mv.norm1(av(0,mv.getNumVectors()));
551  break;
552  case TwoNorm:
553  mv.norm2(av(0,mv.getNumVectors()));
554  break;
555  case InfNorm:
556  mv.normInf(av(0,mv.getNumVectors()));
557  break;
558  default:
559  // Throw logic_error rather than invalid_argument, because if
560  // we get here, it's probably the fault of a Belos solver,
561  // rather than a user giving Belos an invalid input.
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.");
570  }
571  }
572 
573  static void
574  SetBlock (const MV& A, const std::vector<int>& index, MV& mv)
575  {
576  using Teuchos::Range1D;
577  using Teuchos::RCP;
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));
590  ::Tpetra::deep_copy (*mvsub, *Asub);
591  } else {
592  ::Tpetra::deep_copy (*mvsub, A);
593  }
594  }
595 
596  static void
597  SetBlock (const MV& A, const Teuchos::Range1D& index, MV& mv)
598  {
599  // Range1D bounds are signed; size_t is unsigned.
600  // Assignment of Tpetra::MultiVector is a deep copy.
601 
602  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
603  // fair to assume that the number of vectors won't overflow int,
604  // since the typical use case of multivectors involves few
605  // columns, but it's friendly to check just in case.
606  const size_t maxInt =
607  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
608  const bool overflow =
609  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
610  if (overflow) {
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.");
620  }
621  // We've already validated the static casts above.
622  const int numColsA = static_cast<int> (A.getNumVectors ());
623  const int numColsMv = static_cast<int> (mv.getNumVectors ());
624  // 'index' indexes into mv; it's the index set of the target.
625  const bool validIndex =
626  index.lbound () >= 0 && index.ubound () < numColsMv;
627  // We can't take more columns out of A than A has.
628  const bool validSource = index.size () <= numColsA;
629 
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!");
647  }
648 
649  // View of the relevant column(s) of the target multivector mv.
650  // We avoid view creation overhead by only creating a view if
651  // the index range is different than [0, (# columns in mv) - 1].
652  Teuchos::RCP<MV> mv_view;
653  if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
654  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
655  } else {
656  mv_view = CloneViewNonConst (mv, index);
657  }
658 
659  // View of the relevant column(s) of the source multivector A.
660  // If A has fewer columns than mv_view, then create a view of
661  // the first index.size() columns of A.
662  Teuchos::RCP<const MV> A_view;
663  if (index.size () == numColsA) {
664  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
665  } else {
666  A_view = CloneView (A, Teuchos::Range1D (0, index.size () - 1));
667  }
668 
669  ::Tpetra::deep_copy (*mv_view, *A_view);
670  }
671 
672  static void Assign (const MV& A, MV& mv)
673  {
674  const char errPrefix[] = "Belos::MultiVecTraits::Assign(A, mv): ";
675 
676  // Range1D bounds are signed; size_t is unsigned.
677  // Assignment of Tpetra::MultiVector is a deep copy.
678 
679  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
680  // fair to assume that the number of vectors won't overflow int,
681  // since the typical use case of multivectors involves few
682  // columns, but it's friendly to check just in case.
683  const size_t maxInt =
684  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
685  const bool overflow =
686  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
687  if (overflow) {
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!");
698  }
699  // We've already validated the static casts above.
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.");
707  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
708  }
709  if (numColsA == numColsMv) {
710  ::Tpetra::deep_copy (mv, A);
711  } else {
712  Teuchos::RCP<MV> mv_view =
713  CloneViewNonConst (mv, Teuchos::Range1D (0, numColsA - 1));
714  ::Tpetra::deep_copy (*mv_view, A);
715  }
716  }
717 
718 
719  static void MvRandom (MV& mv) {
720  mv.randomize ();
721  }
722 
723  static void
724  MvInit (MV& mv, const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
725  {
726  mv.putScalar (alpha);
727  }
728 
729  static void MvPrint (const MV& mv, std::ostream& os) {
730  Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
731  mv.describe (fos, Teuchos::VERB_EXTREME);
732  }
733 
734 #ifdef HAVE_BELOS_TSQR
735  typedef Tpetra::TsqrAdaptor< Tpetra::MultiVector< Scalar, LO, GO, Node > > tsqr_adaptor_type;
739 #endif // HAVE_BELOS_TSQR
740  };
741 
743  //
744  // Implementation of the Belos::OperatorTraits for Tpetra::Operator.
745  //
747 
749  template <class Storage, class LO, class GO, class Node>
750  class OperatorTraits <typename Storage::value_type,
751  Tpetra::MultiVector<Sacado::UQ::PCE<Storage>,
752  LO,GO,Node>,
753  Tpetra::Operator<Sacado::UQ::PCE<Storage>,
754  LO,GO,Node> >
755  {
756  public:
758  static void
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)
763  {
764  switch (trans) {
765  case NOTRANS:
766  Op.apply(X,Y,Teuchos::NO_TRANS);
767  break;
768  case TRANS:
769  Op.apply(X,Y,Teuchos::TRANS);
770  break;
771  case CONJTRANS:
772  Op.apply(X,Y,Teuchos::CONJ_TRANS);
773  break;
774  default:
775  const std::string scalarName = Teuchos::TypeNameTraits<Scalar>::name();
776  const std::string loName = Teuchos::TypeNameTraits<LO>::name();
777  const std::string goName = Teuchos::TypeNameTraits<GO>::name();
778  const std::string nodeName = Teuchos::TypeNameTraits<Node>::name();
779  const std::string otName = "Belos::OperatorTraits<" + scalarName
780  + "," + loName + "," + goName + "," + nodeName + ">";
781  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, otName << ": Should never "
782  "get here; fell through a switch statement. "
783  "Please report this bug to the Belos developers.");
784  }
785  }
786 
787  static bool
788  HasApplyTranspose (const Tpetra::Operator<Scalar,LO,GO,Node>& Op)
789  {
790  return Op.hasTransposeApply ();
791  }
792  };
793 
794 } // end of Belos namespace
795 
796 #endif
797 // end of file BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
ScalarType * values() const
Kokkos::DefaultExecutionSpace execution_space
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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)
Ordinal ubound() const
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< Scalar > &alphas)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)
Ordinal lbound() const
OrdinalType numCols() const
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)
Ordinal size() const
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]).
OrdinalType stride() const
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)