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_MP_Vector.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_MP_VECTOR_HPP
43 #define BELOS_TPETRA_ADAPTER_MP_VECTOR_HPP
44 
45 #include "BelosTpetraAdapter.hpp"
47 #include "KokkosBlas.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::MP::Vector<Storage>,
74  LO, GO, Node > > {
75  public:
76  typedef typename Storage::ordinal_type s_ordinal;
77  typedef typename Storage::value_type BaseScalar;
79  typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
80  public:
81  typedef typename Tpetra::MultiVector<Scalar,LO,GO,Node>::dot_type dot_type;
82  typedef typename Tpetra::MultiVector<Scalar,LO,GO,Node>::mag_type mag_type;
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 ();
130  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
131  std::runtime_error, fnName << ": All indices must be nonnegative.");
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() << "]): ";
170  index.size() == 0, std::invalid_argument,
171  os.str() << "Empty index range is not allowed.");
173  index.lbound() < 0, std::invalid_argument,
174  os.str() << "Index range includes negative index/ices, which is not "
175  "allowed.");
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 ();
196  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
197  std::invalid_argument,
198  fnName << ": All indices must be nonnegative.");
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>
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() << "]): ";
234  index.size() == 0, std::invalid_argument,
235  os.str() << "Empty index range is not allowed.");
237  index.lbound() < 0, std::invalid_argument,
238  os.str() << "Index range includes negative inde{x,ices}, which is "
239  "not allowed.");
241  index.ubound() >= numCols, std::invalid_argument,
242  os.str() << "Index range exceeds number of vectors " << numCols
243  << " in the input multivector.");
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 
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 ();
261  *std::min_element (index.begin (), index.end ()) < 0,
262  std::invalid_argument,
263  fnName << ": All indices must be nonnegative.");
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 
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.");
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,
328  const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
330  const dot_type& beta,
331  Tpetra::MultiVector<Scalar,LO,GO,Node>& C)
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
354  typedef Tpetra::MultiVector<dot_type,LO,GO,Node> FMV;
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_input( B.values(), strideB, numColsB);
363  auto B_view_host = Kokkos::subview(B_view_host_input,
364  Kokkos::pair<int,int>(0,numRowsB),
365  Kokkos::pair<int,int>(0,numColsB));
366 
367  // Create view for B on the device -- need to be careful to get the
368  // right stride to match B
369  typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> b_view_type;
370  typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> b_1d_view_type;
371  b_1d_view_type B_1d_view_dev(Kokkos::ViewAllocateWithoutInitializing("B"), numRowsB*numColsB);
372  b_view_type B_view_dev( B_1d_view_dev.data(), numRowsB, numColsB);
373  Kokkos::deep_copy(B_view_dev, B_view_host);
374 
375  // Do local multiply
376  {
377  const char ctransA = 'N', ctransB = 'N';
378  KokkosBlas::gemm (
379  &ctransA, &ctransB,
380  alpha, flat_A_view, B_view_dev, beta, flat_C_view);
381  }
382  // Copy back to C if we made a copy
383  if (C.isConstantStride() == false)
384  C.assign(*Ctmp);
385  }
386 
394  static void
395  MvAddMv (Scalar alpha,
396  const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
397  Scalar beta,
398  const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
399  Tpetra::MultiVector<Scalar,LO,GO,Node>& mv)
400  {
401  mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
402  }
403 
404  static void
405  MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
406  const Scalar& alpha)
407  {
408  mv.scale (alpha);
409  }
410 
411  static void
412  MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
413  const std::vector<BaseScalar>& alphas)
414  {
415  std::vector<Scalar> alphas_mp(alphas.size());
416  const size_t sz = alphas.size();
417  for (size_t i=0; i<sz; ++i)
418  alphas_mp[i] = alphas[i];
419  mv.scale (alphas_mp);
420  }
421 
422  static void
423  MvScale (Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
424  const std::vector<Scalar>& alphas)
425  {
426  mv.scale (alphas);
427  }
428 
429  static void
431  const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
432  const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
434  {
435  using Teuchos::Comm;
436  using Teuchos::RCP;
437  using Teuchos::rcp;
438  using Teuchos::REDUCE_SUM;
439  using Teuchos::reduceAll;
440 
441  // Check if numRowsC == numColsC == 1, in which case we can call dot()
442  const int numRowsC = C.numRows ();
443  const int numColsC = C.numCols ();
444  const int strideC = C.stride ();
445  if (numRowsC == 1 && numColsC == 1) {
446  if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
447  // Short-circuit, as required by BLAS semantics.
448  C(0,0) = alpha;
449  return;
450  }
451  A.dot (B, Teuchos::ArrayView<dot_type> (C.values (), 1));
452  if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
453  C(0,0) *= alpha;
454  }
455  return;
456  }
457 
458  // Ensure A and B have constant stride
459  RCP<const MV> Atmp, Btmp;
460  if (A.isConstantStride() == false) Atmp = rcp (new MV (A, Teuchos::Copy));
461  else Atmp = rcp(&A,false);
462 
463  if (B.isConstantStride() == false) Btmp = rcp (new MV (B, Teuchos::Copy));
464  else Btmp = rcp(&B,false);
465 
466  // Create flattened Kokkos::MultiVector's
467  typedef Tpetra::MultiVector<dot_type,LO,GO,Node> FMV;
468  typedef typename FMV::dual_view_type::t_dev flat_view_type;
470  flat_view_type flat_A_view = Atmp->template getLocalView<execution_space>();
471  flat_view_type flat_B_view = Btmp->template getLocalView<execution_space>();
472 
473  // Create a view for C on the host
474  typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, Kokkos::HostSpace> c_host_view_type;
475  c_host_view_type C_view_host_input( C.values(), strideC, numColsC);
476  auto C_view_host = Kokkos::subview(C_view_host_input,
477  Kokkos::pair<int,int>(0,numRowsC),
478  Kokkos::pair<int,int>(0,numColsC));
479 
480  // Create view for C on the device -- need to be careful to get the
481  // right stride to match C (allow setting to 0 for first-touch)
482  typedef Kokkos::View<dot_type**, Kokkos::LayoutLeft, execution_space> c_view_type;
483  typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, execution_space> c_1d_view_type;
484  c_1d_view_type C_1d_view_dev("C", numRowsC*numColsC);
485  c_view_type C_view_dev( C_1d_view_dev.data(), numRowsC, numColsC);
486 
487  // Do local multiply
488  {
489  const char ctransA = 'C', ctransB = 'N';
490  KokkosBlas::gemm (
491  &ctransA, &ctransB,
492  alpha, flat_A_view, flat_B_view,
493  Kokkos::Details::ArithTraits<dot_type>::zero(),
494  C_view_dev);
495  }
496 
497  // reduce across processors -- could check for RDMA
498  RCP<const Comm<int> > pcomm = A.getMap()->getComm ();
499  if (pcomm->getSize () == 1)
500  Kokkos::deep_copy(C_view_host, C_view_dev);
501  else {
502  typedef Kokkos::View<dot_type*, Kokkos::LayoutLeft, Kokkos::HostSpace> c_1d_host_view_type;
503  c_1d_host_view_type C_1d_view_tmp(Kokkos::ViewAllocateWithoutInitializing("C_tmp"), strideC*numColsC);
504  c_host_view_type C_view_tmp( C_1d_view_tmp.data(),
505  strideC, numColsC);
506  Kokkos::deep_copy(Kokkos::subview(C_view_tmp,
507  Kokkos::pair<int,int>(0,numRowsC),
508  Kokkos::pair<int,int>(0,numColsC)) ,
509  C_view_dev);
510  reduceAll<int> (*pcomm, REDUCE_SUM, strideC*numColsC,
511  C_view_tmp.data(),
512  C_view_host.data());
513  }
514  }
515 
517  static void
518  MvDot (const Tpetra::MultiVector<Scalar,LO,GO,Node>& A,
519  const Tpetra::MultiVector<Scalar,LO,GO,Node>& B,
520  std::vector<dot_type>& dots)
521  {
522  const size_t numVecs = A.getNumVectors ();
523 
525  numVecs != B.getNumVectors (), std::invalid_argument,
526  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): "
527  "A and B must have the same number of columns. "
528  "A has " << numVecs << " column(s), "
529  "but B has " << B.getNumVectors () << " column(s).");
530 #ifdef HAVE_TPETRA_DEBUG
532  dots.size() < numVecs, std::invalid_argument,
533  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): "
534  "The output array 'dots' must have room for all dot products. "
535  "A and B each have " << numVecs << " column(s), "
536  "but 'dots' only has " << dots.size() << " entry(/ies).");
537 #endif // HAVE_TPETRA_DEBUG
538 
540  A.dot (B, av (0, numVecs));
541  }
542 
544  static void
545  MvNorm (const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv,
546  std::vector<mag_type> &normvec,
547  NormType type=TwoNorm)
548  {
549 
550 #ifdef HAVE_TPETRA_DEBUG
551  typedef std::vector<int>::size_type size_type;
553  normvec.size () < static_cast<size_type> (mv.getNumVectors ()),
554  std::invalid_argument,
555  "Belos::MultiVecTraits::MvNorm(mv,normvec): The normvec output "
556  "argument must have at least as many entries as the number of vectors "
557  "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
558  << " < mv.getNumVectors() = " << mv.getNumVectors () << ".");
559 #endif
560 
561  Teuchos::ArrayView<mag_type> av(normvec);
562  switch (type) {
563  case OneNorm:
564  mv.norm1(av(0,mv.getNumVectors()));
565  break;
566  case TwoNorm:
567  mv.norm2(av(0,mv.getNumVectors()));
568  break;
569  case InfNorm:
570  mv.normInf(av(0,mv.getNumVectors()));
571  break;
572  default:
573  // Throw logic_error rather than invalid_argument, because if
574  // we get here, it's probably the fault of a Belos solver,
575  // rather than a user giving Belos an invalid input.
577  true, std::logic_error,
578  "Belos::MultiVecTraits::MvNorm: Invalid NormType value " << type
579  << ". Valid values are OneNorm=" << OneNorm << ", TwoNorm="
580  << TwoNorm <<", and InfNorm=" << InfNorm << ". If you are a Belos "
581  "user and have not modified Belos in any way, and you get this "
582  "message, then this is probably a bug in the Belos solver you were "
583  "using. Please report this to the Belos developers.");
584  }
585  }
586 
587  static void
588  SetBlock (const MV& A, const std::vector<int>& index, MV& mv)
589  {
590  using Teuchos::Range1D;
591  using Teuchos::RCP;
592  const size_t inNumVecs = A.getNumVectors ();
593 #ifdef HAVE_TPETRA_DEBUG
595  inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
596  "Belos::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must "
597  "have no more entries as the number of columns in the input MultiVector"
598  " A. A.getNumVectors() = " << inNumVecs << " < index.size () = "
599  << index.size () << ".");
600 #endif // HAVE_TPETRA_DEBUG
601  RCP<MV> mvsub = CloneViewNonConst (mv, index);
602  if (inNumVecs > static_cast<size_t> (index.size ())) {
603  RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
604  ::Tpetra::deep_copy (*mvsub, *Asub);
605  } else {
606  ::Tpetra::deep_copy (*mvsub, A);
607  }
608  }
609 
610  static void
611  SetBlock (const MV& A, const Teuchos::Range1D& index, MV& mv)
612  {
613  // Range1D bounds are signed; size_t is unsigned.
614  // Assignment of Tpetra::MultiVector is a deep copy.
615 
616  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
617  // fair to assume that the number of vectors won't overflow int,
618  // since the typical use case of multivectors involves few
619  // columns, but it's friendly to check just in case.
620  const size_t maxInt =
621  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
622  const bool overflow =
623  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
624  if (overflow) {
625  std::ostringstream os;
626  os << "Belos::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
627  << ", " << index.ubound () << "], mv): ";
629  maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
630  "of columns (size_t) in the input MultiVector 'A' overflows int.");
632  maxInt < mv.getNumVectors (), std::range_error, os.str () << "Number "
633  "of columns (size_t) in the output MultiVector 'mv' overflows int.");
634  }
635  // We've already validated the static casts above.
636  const int numColsA = static_cast<int> (A.getNumVectors ());
637  const int numColsMv = static_cast<int> (mv.getNumVectors ());
638  // 'index' indexes into mv; it's the index set of the target.
639  const bool validIndex =
640  index.lbound () >= 0 && index.ubound () < numColsMv;
641  // We can't take more columns out of A than A has.
642  const bool validSource = index.size () <= numColsA;
643 
644  if (! validIndex || ! validSource) {
645  std::ostringstream os;
646  os << "Belos::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
647  << ", " << index.ubound () << "], mv): ";
649  index.lbound() < 0, std::invalid_argument,
650  os.str() << "Range lower bound must be nonnegative.");
652  index.ubound() >= numColsMv, std::invalid_argument,
653  os.str() << "Range upper bound must be less than the number of "
654  "columns " << numColsA << " in the 'mv' output argument.");
656  index.size() > numColsA, std::invalid_argument,
657  os.str() << "Range must have no more elements than the number of "
658  "columns " << numColsA << " in the 'A' input argument.");
660  true, std::logic_error, "Should never get here!");
661  }
662 
663  // View of the relevant column(s) of the target multivector mv.
664  // We avoid view creation overhead by only creating a view if
665  // the index range is different than [0, (# columns in mv) - 1].
666  Teuchos::RCP<MV> mv_view;
667  if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
668  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
669  } else {
670  mv_view = CloneViewNonConst (mv, index);
671  }
672 
673  // View of the relevant column(s) of the source multivector A.
674  // If A has fewer columns than mv_view, then create a view of
675  // the first index.size() columns of A.
676  Teuchos::RCP<const MV> A_view;
677  if (index.size () == numColsA) {
678  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
679  } else {
680  A_view = CloneView (A, Teuchos::Range1D (0, index.size () - 1));
681  }
682 
683  ::Tpetra::deep_copy (*mv_view, *A_view);
684  }
685 
686  static void Assign (const MV& A, MV& mv)
687  {
688  const char errPrefix[] = "Belos::MultiVecTraits::Assign(A, mv): ";
689 
690  // Range1D bounds are signed; size_t is unsigned.
691  // Assignment of Tpetra::MultiVector is a deep copy.
692 
693  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
694  // fair to assume that the number of vectors won't overflow int,
695  // since the typical use case of multivectors involves few
696  // columns, but it's friendly to check just in case.
697  const size_t maxInt =
698  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
699  const bool overflow =
700  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
701  if (overflow) {
703  maxInt < A.getNumVectors(), std::range_error,
704  errPrefix << "Number of columns in the input multivector 'A' "
705  "(a size_t) overflows int.");
707  maxInt < mv.getNumVectors(), std::range_error,
708  errPrefix << "Number of columns in the output multivector 'mv' "
709  "(a size_t) overflows int.");
711  true, std::logic_error, "Should never get here!");
712  }
713  // We've already validated the static casts above.
714  const int numColsA = static_cast<int> (A.getNumVectors ());
715  const int numColsMv = static_cast<int> (mv.getNumVectors ());
716  if (numColsA > numColsMv) {
718  numColsA > numColsMv, std::invalid_argument,
719  errPrefix << "Input multivector 'A' has " << numColsA << " columns, "
720  "but output multivector 'mv' has only " << numColsMv << " columns.");
721  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
722  }
723  if (numColsA == numColsMv) {
724  ::Tpetra::deep_copy (mv, A);
725  } else {
726  Teuchos::RCP<MV> mv_view =
727  CloneViewNonConst (mv, Teuchos::Range1D (0, numColsA - 1));
728  ::Tpetra::deep_copy (*mv_view, A);
729  }
730  }
731 
732 
733  static void MvRandom (MV& mv) {
734  mv.randomize ();
735  }
736 
737  static void
738  MvInit (MV& mv, const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
739  {
740  mv.putScalar (alpha);
741  }
742 
743  static void MvPrint (const MV& mv, std::ostream& os) {
744  Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
745  mv.describe (fos, Teuchos::VERB_EXTREME);
746  }
747 
748 #ifdef HAVE_BELOS_TSQR
749  typedef Tpetra::TsqrAdaptor< Tpetra::MultiVector< Scalar, LO, GO, Node > > tsqr_adaptor_type;
753 #endif // HAVE_BELOS_TSQR
754  };
755 
757  //
758  // Implementation of the Belos::OperatorTraits for Tpetra::Operator.
759  //
761 
763  template <class Storage, class LO, class GO, class Node>
764  class OperatorTraits <typename Storage::value_type,
765  Tpetra::MultiVector<Sacado::MP::Vector<Storage>,
766  LO,GO,Node>,
767  Tpetra::Operator<Sacado::MP::Vector<Storage>,
768  LO,GO,Node> >
769  {
770  public:
772  static void
773  Apply (const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
774  const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
775  Tpetra::MultiVector<Scalar,LO,GO,Node>& Y,
776  ETrans trans=NOTRANS)
777  {
778  switch (trans) {
779  case NOTRANS:
780  Op.apply(X,Y,Teuchos::NO_TRANS);
781  break;
782  case TRANS:
783  Op.apply(X,Y,Teuchos::TRANS);
784  break;
785  case CONJTRANS:
786  Op.apply(X,Y,Teuchos::CONJ_TRANS);
787  break;
788  default:
789  const std::string scalarName = Teuchos::TypeNameTraits<Scalar>::name();
790  const std::string loName = Teuchos::TypeNameTraits<LO>::name();
791  const std::string goName = Teuchos::TypeNameTraits<GO>::name();
792  const std::string nodeName = Teuchos::TypeNameTraits<Node>::name();
793  const std::string otName = "Belos::OperatorTraits<" + scalarName
794  + "," + loName + "," + goName + "," + nodeName + ">";
795  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, otName << ": Should never "
796  "get here; fell through a switch statement. "
797  "Please report this bug to the Belos developers.");
798  }
799  }
800 
801  static bool
802  HasApplyTranspose (const Tpetra::Operator<Scalar,LO,GO,Node>& Op)
803  {
804  return Op.hasTransposeApply ();
805  }
806  };
807 
808 } // end of Belos namespace
809 
810 #endif
811 // end of file BELOS_TPETRA_ADAPTER_HPP
ScalarType * values() 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.
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< BaseScalar > &alphas)
Kokkos::DefaultExecutionSpace execution_space
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Ordinal ubound() 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)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< Scalar > &alphas)
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)
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
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 Teuchos::RCP< MV > Clone(const MV &X, const int numVecs)
Create a new MultiVector with numVecs columns.
Ordinal lbound() const
OrdinalType numCols() const
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 Teuchos::RCP< MV > CloneCopy(const MV &mv, const Teuchos::Range1D &index)
Create and return a deep copy of the given columns of mv.
Ordinal size() const
OrdinalType stride() const
static std::string name()
OrdinalType numRows() const
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 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