Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Belos_TpetraAdapter_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 BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
43 #define BELOS_TPETRA_ADAPTER_UQ_PCE_HPP
44 
45 #include "BelosTpetraAdapter.hpp"
47 #include "Kokkos_MV_GEMM.hpp"
48 
49 #ifdef HAVE_BELOS_TSQR
51 #endif // HAVE_BELOS_TSQR
52 
53 namespace Belos {
54 
56  //
57  // Implementation of Belos::MultiVecTraits for Tpetra::MultiVector.
58  //
60 
71  template<class Storage, class LO, class GO, class Node>
72  class MultiVecTraits<typename Storage::value_type,
73  Tpetra::MultiVector< Sacado::UQ::PCE<Storage>,
74  LO, GO, Node > > {
75  public:
76  typedef typename Storage::ordinal_type s_ordinal;
77  typedef typename Storage::value_type BaseScalar;
80  public:
83 
89  static Teuchos::RCP<MV> Clone (const MV& X, const int numVecs) {
90  Teuchos::RCP<MV> Y (new MV (X.getMap (), numVecs, false));
91  Y->setCopyOrView (Teuchos::View);
92  return Y;
93  }
94 
96  static Teuchos::RCP<MV> CloneCopy (const MV& X)
97  {
98  // Make a deep copy of X. The one-argument copy constructor
99  // does a shallow copy by default; the second argument tells it
100  // to do a deep copy.
101  Teuchos::RCP<MV> X_copy (new MV (X, Teuchos::Copy));
102  // Make Tpetra::MultiVector use the new view semantics. This is
103  // a no-op for the Kokkos refactor version of Tpetra; it only
104  // does something for the "classic" version of Tpetra. This
105  // shouldn't matter because Belos only handles MV through RCP
106  // and through this interface anyway, but it doesn't hurt to set
107  // it and make sure that it works.
108  X_copy->setCopyOrView (Teuchos::View);
109  return X_copy;
110  }
111 
123  static Teuchos::RCP<MV>
124  CloneCopy (const MV& mv, const std::vector<int>& index)
125  {
126 #ifdef HAVE_TPETRA_DEBUG
127  const char fnName[] = "Belos::MultiVecTraits::CloneCopy(mv,index)";
128  const size_t inNumVecs = mv.getNumVectors ();
129  TEUCHOS_TEST_FOR_EXCEPTION(
130  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
131  std::runtime_error, fnName << ": All indices must be nonnegative.");
132  TEUCHOS_TEST_FOR_EXCEPTION(
133  index.size () > 0 &&
134  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
135  std::runtime_error,
136  fnName << ": All indices must be strictly less than the number of "
137  "columns " << inNumVecs << " of the input multivector mv.");
138 #endif // HAVE_TPETRA_DEBUG
139 
140  // Tpetra wants an array of size_t, not of int.
141  Teuchos::Array<size_t> columns (index.size ());
142  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
143  columns[j] = index[j];
144  }
145  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
146  // continuous column index range in MultiVector::subCopy, so we
147  // don't have to check here.
148  Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
149  X_copy->setCopyOrView (Teuchos::View);
150  return X_copy;
151  }
152 
159  static Teuchos::RCP<MV>
160  CloneCopy (const MV& mv, const Teuchos::Range1D& index)
161  {
162  const bool validRange = index.size() > 0 &&
163  index.lbound() >= 0 &&
164  index.ubound() < GetNumberVecs(mv);
165  if (! validRange) { // invalid range; generate error message
166  std::ostringstream os;
167  os << "Belos::MultiVecTraits::CloneCopy(mv,index=["
168  << index.lbound() << "," << index.ubound() << "]): ";
169  TEUCHOS_TEST_FOR_EXCEPTION(
170  index.size() == 0, std::invalid_argument,
171  os.str() << "Empty index range is not allowed.");
172  TEUCHOS_TEST_FOR_EXCEPTION(
173  index.lbound() < 0, std::invalid_argument,
174  os.str() << "Index range includes negative index/ices, which is not "
175  "allowed.");
176  TEUCHOS_TEST_FOR_EXCEPTION(
177  index.ubound() >= GetNumberVecs(mv), std::invalid_argument,
178  os.str() << "Index range exceeds number of vectors "
179  << mv.getNumVectors() << " in the input multivector.");
180  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
181  os.str() << "Should never get here!");
182  }
183  Teuchos::RCP<MV> X_copy = mv.subCopy (index);
184  X_copy->setCopyOrView (Teuchos::View);
185  return X_copy;
186  }
187 
188 
189  static Teuchos::RCP<MV>
190  CloneViewNonConst (MV& mv, const std::vector<int>& index)
191  {
192 #ifdef HAVE_TPETRA_DEBUG
193  const char fnName[] = "Belos::MultiVecTraits::CloneViewNonConst(mv,index)";
194  const size_t numVecs = mv.getNumVectors ();
195  TEUCHOS_TEST_FOR_EXCEPTION(
196  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
197  std::invalid_argument,
198  fnName << ": All indices must be nonnegative.");
199  TEUCHOS_TEST_FOR_EXCEPTION(
200  index.size () > 0 &&
201  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
202  std::invalid_argument,
203  fnName << ": All indices must be strictly less than the number of "
204  "columns " << numVecs << " in the input MultiVector mv.");
205 #endif // HAVE_TPETRA_DEBUG
206 
207  // Tpetra wants an array of size_t, not of int.
208  Teuchos::Array<size_t> columns (index.size ());
209  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
210  columns[j] = index[j];
211  }
212  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
213  // continuous column index range in
214  // MultiVector::subViewNonConst, so we don't have to check here.
215  Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
216  X_view->setCopyOrView (Teuchos::View);
217  return X_view;
218  }
219 
220  static Teuchos::RCP<MV>
221  CloneViewNonConst (MV& mv, const Teuchos::Range1D& index)
222  {
223  // NOTE (mfh 11 Jan 2011) We really should check for possible
224  // overflow of int here. However, the number of columns in a
225  // multivector typically fits in an int.
226  const int numCols = static_cast<int> (mv.getNumVectors());
227  const bool validRange = index.size() > 0 &&
228  index.lbound() >= 0 && index.ubound() < numCols;
229  if (! validRange) {
230  std::ostringstream os;
231  os << "Belos::MultiVecTraits::CloneViewNonConst(mv,index=["
232  << index.lbound() << ", " << index.ubound() << "]): ";
233  TEUCHOS_TEST_FOR_EXCEPTION(
234  index.size() == 0, std::invalid_argument,
235  os.str() << "Empty index range is not allowed.");
236  TEUCHOS_TEST_FOR_EXCEPTION(
237  index.lbound() < 0, std::invalid_argument,
238  os.str() << "Index range includes negative inde{x,ices}, which is "
239  "not allowed.");
240  TEUCHOS_TEST_FOR_EXCEPTION(
241  index.ubound() >= numCols, std::invalid_argument,
242  os.str() << "Index range exceeds number of vectors " << numCols
243  << " in the input multivector.");
244  TEUCHOS_TEST_FOR_EXCEPTION(
245  true, std::logic_error,
246  os.str() << "Should never get here!");
247  }
248  Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
249  X_view->setCopyOrView (Teuchos::View);
250  return X_view;
251  }
252 
253  static Teuchos::RCP<const MV>
254  CloneView (const MV& mv, const std::vector<int>& index)
255  {
256 #ifdef HAVE_TPETRA_DEBUG
257  const char fnName[] = "Belos::MultiVecTraits<Scalar, "
258  "Tpetra::MultiVector<...> >::CloneView(mv,index)";
259  const size_t numVecs = mv.getNumVectors ();
260  TEUCHOS_TEST_FOR_EXCEPTION(
261  *std::min_element (index.begin (), index.end ()) < 0,
262  std::invalid_argument,
263  fnName << ": All indices must be nonnegative.");
264  TEUCHOS_TEST_FOR_EXCEPTION(
265  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
266  std::invalid_argument,
267  fnName << ": All indices must be strictly less than the number of "
268  "columns " << numVecs << " in the input MultiVector mv.");
269 #endif // HAVE_TPETRA_DEBUG
270 
271  // Tpetra wants an array of size_t, not of int.
272  Teuchos::Array<size_t> columns (index.size ());
273  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
274  columns[j] = index[j];
275  }
276  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
277  // continuous column index range in MultiVector::subView, so we
278  // don't have to check here.
279  Teuchos::RCP<const MV> X_view = mv.subView (columns);
280  Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
281  return X_view;
282  }
283 
284  static Teuchos::RCP<const MV>
285  CloneView (const MV& mv, const Teuchos::Range1D& index)
286  {
287  // NOTE (mfh 11 Jan 2011) We really should check for possible
288  // overflow of int here. However, the number of columns in a
289  // multivector typically fits in an int.
290  const int numCols = static_cast<int> (mv.getNumVectors());
291  const bool validRange = index.size() > 0 &&
292  index.lbound() >= 0 && index.ubound() < numCols;
293  if (! validRange) {
294  std::ostringstream os;
295  os << "Belos::MultiVecTraits::CloneView(mv, index=["
296  << index.lbound () << ", " << index.ubound() << "]): ";
297  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
298  os.str() << "Empty index range is not allowed.");
299  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
300  os.str() << "Index range includes negative index/ices, which is not "
301  "allowed.");
302  TEUCHOS_TEST_FOR_EXCEPTION(
303  index.ubound() >= numCols, std::invalid_argument,
304  os.str() << "Index range exceeds number of vectors " << numCols
305  << " in the input multivector.");
306  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
307  os.str() << "Should never get here!");
308  }
309  Teuchos::RCP<const MV> X_view = mv.subView (index);
310  Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
311  return X_view;
312  }
313 
314  static ptrdiff_t GetGlobalLength (const MV& mv) {
315  return static_cast<ptrdiff_t> (mv.getGlobalLength ());
316  }
317 
318  static int GetNumberVecs (const MV& mv) {
319  return static_cast<int> (mv.getNumVectors ());
320  }
321 
322  static bool HasConstantStride (const MV& mv) {
323  return mv.isConstantStride ();
324  }
325 
326  static void
327  MvTimesMatAddMv (const dot_type& alpha,
329  const Teuchos::SerialDenseMatrix<int,dot_type>& B,
330  const dot_type& beta,
332  {
333  using Teuchos::RCP;
334  using Teuchos::rcp;
335 
336  // Check if numRowsB == numColsB == 1, in which case we can call update()
337  const int numRowsB = B.numRows ();
338  const int numColsB = B.numCols ();
339  const int strideB = B.stride ();
340  if (numRowsB == 1 && numColsB == 1) {
341  C.update (alpha*B(0,0), A, beta);
342  return;
343  }
344 
345  // Ensure A and C have constant stride
346  RCP<const MV> Atmp, Ctmp;
347  if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
348  else Atmp = rcp(&A,false);
349 
350  if (C.isConstantStride() == false) Ctmp = rcp (new MV (C, Teuchos::Copy));
351  else Ctmp = rcp(&C,false);
352 
353  // Create flattened view's
355  typedef typename FMV::dual_view_type::t_dev flat_view_type;
357  flat_view_type flat_A_view = Atmp->template getLocalView<execution_space>();
358  flat_view_type flat_C_view = Ctmp->template getLocalView<execution_space>();
359 
360  // Create a view for B on the host
361  typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> b_host_view_type;
362  b_host_view_type B_view_host( B.values(), strideB, numColsB);
363 
364  // Create view for B on the device -- need to be careful to get the
365  // right stride to match B
366  typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> b_view_type;
367  typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> b_1d_view_type;
368  b_1d_view_type B_1d_view_dev(Kokkos::ViewAllocateWithoutInitializing("B"), strideB*numColsB);
369  b_view_type B_view_dev( B_1d_view_dev.ptr_on_device(), strideB, numColsB);
370  Kokkos::deep_copy(B_view_dev, B_view_host);
371 
372  // Do local multiply
373  Kokkos::DeviceGEMM<dot_type,execution_space>::GEMM(
374  Teuchos::NO_TRANS, Teuchos::NO_TRANS,
375  alpha, flat_A_view, B_view_dev, beta, flat_C_view);
376 
377  // Copy back to C if we made a copy
378  if (C.isConstantStride() == false)
379  C.assign(*Ctmp);
380  }
381 
389  static void
390  MvAddMv (Scalar alpha,
392  Scalar beta,
395  {
396  mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
397  }
398 
399  static void
401  const Scalar& alpha)
402  {
403  mv.scale (alpha);
404  }
405 
406  static void
408  const std::vector<BaseScalar>& alphas)
409  {
410  std::vector<Scalar> alphas_mp(alphas.size());
411  const size_t sz = alphas.size();
412  for (size_t i=0; i<sz; ++i)
413  alphas_mp[i] = alphas[i];
414  mv.scale (alphas_mp);
415  }
416 
417  static void
419  const std::vector<Scalar>& alphas)
420  {
421  mv.scale (alphas);
422  }
423 
424  static void
428  Teuchos::SerialDenseMatrix<int,dot_type>& C)
429  {
430  using Teuchos::Comm;
431  using Teuchos::RCP;
432  using Teuchos::rcp;
433  using Teuchos::REDUCE_SUM;
434  using Teuchos::reduceAll;
435 
436  // Check if numRowsC == numColsC == 1, in which case we can call dot()
437  const int numRowsC = C.numRows ();
438  const int numColsC = C.numCols ();
439  const int strideC = C.stride ();
440  if (numRowsC == 1 && numColsC == 1) {
441  if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
442  // Short-circuit, as required by BLAS semantics.
443  C(0,0) = alpha;
444  return;
445  }
446  A.dot (B, Teuchos::ArrayView<dot_type> (C.values (), 1));
447  if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
448  C(0,0) *= alpha;
449  }
450  return;
451  }
452 
453  // Ensure A and B have constant stride
454  RCP<const MV> Atmp, Btmp;
455  if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
456  else Atmp = rcp(&A,false);
457 
458  if (B.isConstantStride() == false) Btmp = rcp (new MV (B, Teuchos::Copy));
459  else Btmp = rcp(&B,false);
460 
461  // Create flattened Kokkos::MultiVector's
463  typedef typename FMV::dual_view_type::t_dev flat_view_type;
465  flat_view_type flat_A_view = Atmp->template getLocalView<execution_space>();
466  flat_view_type flat_B_view = Btmp->template getLocalView<execution_space>();
467 
468  // Create a view for C on the host
469  typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> c_host_view_type;
470  c_host_view_type C_view_host( C.values(), strideC, numColsC);
471 
472  // Create view for C on the device -- need to be careful to get the
473  // right stride to match C (allow setting to 0 for first-touch)
474  typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> c_view_type;
475  typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> c_1d_view_type;
476  c_1d_view_type C_1d_view_dev("C", strideC*numColsC);
477  c_view_type C_view_dev( C_1d_view_dev.ptr_on_device(), strideC, numColsC);
478 
479  // Do local multiply
480  Kokkos::DeviceGEMM<dot_type,execution_space>::GEMM(
481  Teuchos::CONJ_TRANS, Teuchos::NO_TRANS,
482  alpha, flat_A_view, flat_B_view,
483  Kokkos::Details::ArithTraits<dot_type>::zero(),
484  C_view_dev);
485 
486  // reduce across processors -- could check for RDMA
487  RCP<const Comm<int> > pcomm = A.getMap()->getComm ();
488  if (pcomm->getSize () == 1)
489  Kokkos::deep_copy(C_view_host, C_view_dev);
490  else {
491  typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, Kokkos::HostSpace> c_1d_host_view_type;
492  c_1d_host_view_type C_1d_view_tmp(Kokkos::ViewAllocateWithoutInitializing("C_tmp"), strideC*numColsC);
493  c_host_view_type C_view_tmp( C_1d_view_tmp.ptr_on_device(),
494  strideC, numColsC);
495  Kokkos::deep_copy(C_view_tmp, C_view_dev);
496  reduceAll<int> (*pcomm, REDUCE_SUM, strideC*numColsC,
497  C_view_tmp.ptr_on_device(),
498  C_view_host.ptr_on_device());
499  }
500  }
501 
503  static void
506  std::vector<dot_type>& dots)
507  {
508  const size_t numVecs = A.getNumVectors ();
509 
510  TEUCHOS_TEST_FOR_EXCEPTION(
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
517  TEUCHOS_TEST_FOR_EXCEPTION(
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 
525  Teuchos::ArrayView<dot_type> av (dots);
526  A.dot (B, av (0, numVecs));
527  }
528 
530  static void
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;
538  TEUCHOS_TEST_FOR_EXCEPTION(
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.
562  TEUCHOS_TEST_FOR_EXCEPTION(
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
580  TEUCHOS_TEST_FOR_EXCEPTION(
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): ";
614  TEUCHOS_TEST_FOR_EXCEPTION(
615  maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
616  "of columns (size_t) in the input MultiVector 'A' overflows int.");
617  TEUCHOS_TEST_FOR_EXCEPTION(
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): ";
634  TEUCHOS_TEST_FOR_EXCEPTION(
635  index.lbound() < 0, std::invalid_argument,
636  os.str() << "Range lower bound must be nonnegative.");
637  TEUCHOS_TEST_FOR_EXCEPTION(
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.");
641  TEUCHOS_TEST_FOR_EXCEPTION(
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.");
645  TEUCHOS_TEST_FOR_EXCEPTION(
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) {
688  TEUCHOS_TEST_FOR_EXCEPTION(
689  maxInt < A.getNumVectors(), std::range_error,
690  errPrefix << "Number of columns in the input multivector 'A' "
691  "(a size_t) overflows int.");
692  TEUCHOS_TEST_FOR_EXCEPTION(
693  maxInt < mv.getNumVectors(), std::range_error,
694  errPrefix << "Number of columns in the output multivector 'mv' "
695  "(a size_t) overflows int.");
696  TEUCHOS_TEST_FOR_EXCEPTION(
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) {
703  TEUCHOS_TEST_FOR_EXCEPTION(
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,
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
Kokkos::DefaultExecutionSpace execution_space
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)
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< Scalar > &alphas)
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)
KOKKOS_INLINE_FUNCTION PCE< Storage > max(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
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)
expr expr expr expr j
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)
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.
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 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)