Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
BelosPCETpetraAdapter.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_PCE_TPETRA_ADAPTER_HPP
43 #define BELOS_PCE_TPETRA_ADAPTER_HPP
44 
49 // TODO: the assumption is made that the solver, multivector and operator are templated on the same scalar. this will need to be modified.
50 
51 #include <Tpetra_MultiVector.hpp>
52 #include <Tpetra_Operator.hpp>
53 #include <Teuchos_Assert.hpp>
54 #include <Teuchos_ScalarTraits.hpp>
55 #include <Teuchos_TypeNameTraits.hpp>
56 #include <Teuchos_Array.hpp>
57 #include <Teuchos_DefaultSerialComm.hpp>
58 
59 #include <BelosConfigDefs.hpp>
60 #include <BelosTypes.hpp>
61 #include <BelosMultiVecTraits.hpp>
62 #include <BelosOperatorTraits.hpp>
63 #include <Kokkos_NodeAPIConfigDefs.hpp>
64 
65 #ifdef HAVE_BELOS_TSQR
66 # include <Tpetra_TsqrAdaptor.hpp>
67 #endif // HAVE_BELOS_TSQR
68 
69 
70 namespace Belos {
71 
73  //
74  // Implementation of the Belos::MultiVecTraits for Tpetra::MultiVector.
75  //
77 
86  template<class BaseScalar, class Storage, class LO, class GO, class Node>
87  class MultiVecTraits<BaseScalar, Tpetra::MultiVector< Sacado::PCE::OrthogPoly<BaseScalar,Storage>,LO,GO,Node> >
88  {
89  public:
90 #ifdef HAVE_BELOS_TPETRA_TIMERS
91  static Teuchos::RCP<Teuchos::Time> mvTimesMatAddMvTimer_, mvTransMvTimer_;
92 #endif
93 
95 
96  static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > Clone( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const int numvecs )
97  {
98  return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>(mv.getMap(),numvecs));
99  }
100 
101  static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
102  {
103  return Teuchos::rcp( new Tpetra::MultiVector<Scalar,LO,GO,Node>( mv ) );
104  }
105 
106  static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneCopy( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
107  {
108  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
109  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneCopy(mv,index): numvecs must be greater than zero.");
110 #ifdef HAVE_TPETRA_DEBUG
111  TEUCHOS_TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::runtime_error,
112  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneCopy(mv,index): indices must be >= zero.");
113  TEUCHOS_TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::runtime_error,
114  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneCopy(mv,index): indices must be < mv.getNumVectors().");
115 #endif
116  for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
117  if (index[j] != index[j-1]+1) {
118  // not contiguous; short circuit
119  Teuchos::Array<size_t> stinds(index.begin(), index.end());
120  return mv.subCopy(stinds);
121  }
122  }
123  // contiguous
124  return mv.subCopy(Teuchos::Range1D(index.front(),index.back()));
125  }
126 
127  static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> >
129  const Teuchos::Range1D& index)
130  {
131  const bool validRange = index.size() > 0 &&
132  index.lbound() >= 0 &&
133  index.ubound() < GetNumberVecs(mv);
134  if (! validRange)
135  {
136  std::ostringstream os;
137  os << "Belos::MultiVecTraits<Scalar, Tpetra::MultiVector<...> >::"
138  "CloneCopy(mv,index=[" << index.lbound() << ", " << index.ubound()
139  << "]): ";
140  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
141  os.str() << "Empty index range is not allowed.");
142  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
143  os.str() << "Index range includes negative "
144  "index/ices, which is not allowed.");
145  // Range1D bounds are signed; size_t is unsigned.
146  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= GetNumberVecs(mv),
147  std::invalid_argument,
148  os.str() << "Index range exceeds number of vectors "
149  << mv.getNumVectors() << " in the input multivector.");
150  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
151  os.str() << "Should never get here!");
152  }
153  return mv.subCopy (index);
154  }
155 
156 
157  static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneViewNonConst( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
158  {
159  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
160  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
161 #ifdef HAVE_TPETRA_DEBUG
162  TEUCHOS_TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
163  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be >= zero.");
164  TEUCHOS_TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
165  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be < mv.getNumVectors().");
166 #endif
167  for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
168  if (index[j] != index[j-1]+1) {
169  // not contiguous; short circuit
170  Teuchos::Array<size_t> stinds(index.begin(), index.end());
171  return mv.subViewNonConst(stinds);
172  }
173  }
174  // contiguous
175  return mv.subViewNonConst(Teuchos::Range1D(index.front(),index.back()));
176  }
177 
178 
179  static Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> >
181  const Teuchos::Range1D& index)
182  {
183  // NOTE (mfh 11 Jan 2011) We really should check for possible
184  // overflow of int here. However, the number of columns in a
185  // multivector typically fits in an int.
186  const int numCols = static_cast<int> (mv.getNumVectors());
187  const bool validRange = index.size() > 0 &&
188  index.lbound() >= 0 && index.ubound() < numCols;
189  if (! validRange)
190  {
191  std::ostringstream os;
192  os << "Belos::MultiVecTraits<Scalar, Tpetra::MultiVector<...> >::"
193  "CloneViewNonConst(mv,index=[" << index.lbound() << ", "
194  << index.ubound() << "]): ";
195  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
196  os.str() << "Empty index range is not allowed.");
197  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
198  os.str() << "Index range includes negative "
199  "index/ices, which is not allowed.");
200  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numCols, std::invalid_argument,
201  os.str() << "Index range exceeds number of "
202  "vectors " << numCols << " in the input "
203  "multivector.");
204  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
205  os.str() << "Should never get here!");
206  }
207  return mv.subViewNonConst (index);
208  }
209 
210 
211  static Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > CloneView(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<int>& index )
212  {
213  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0,std::invalid_argument,
214  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): numvecs must be greater than zero.");
215 #ifdef HAVE_TPETRA_DEBUG
216  TEUCHOS_TEST_FOR_EXCEPTION( *std::min_element(index.begin(),index.end()) < 0, std::invalid_argument,
217  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be >= zero.");
218  TEUCHOS_TEST_FOR_EXCEPTION( (size_t)*std::max_element(index.begin(),index.end()) >= mv.getNumVectors(), std::invalid_argument,
219  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::CloneView(mv,index): indices must be < mv.getNumVectors().");
220 #endif
221  for (typename std::vector<int>::size_type j=1; j<index.size(); ++j) {
222  if (index[j] != index[j-1]+1) {
223  // not contiguous; short circuit
224  Teuchos::Array<size_t> stinds(index.begin(), index.end());
225  return mv.subView(stinds);
226  }
227  }
228  // contiguous
229  return mv.subView(Teuchos::Range1D(index.front(),index.back()));
230  }
231 
232  static Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> >
234  const Teuchos::Range1D& index)
235  {
236  // NOTE (mfh 11 Jan 2011) We really should check for possible
237  // overflow of int here. However, the number of columns in a
238  // multivector typically fits in an int.
239  const int numCols = static_cast<int> (mv.getNumVectors());
240  const bool validRange = index.size() > 0 &&
241  index.lbound() >= 0 && index.ubound() < numCols;
242  if (! validRange)
243  {
244  std::ostringstream os;
245  os << "Belos::MultiVecTraits<Scalar, Tpetra::MultiVector<...> >::"
246  "CloneView(mv, index=[" << index.lbound() << ", "
247  << index.ubound() << "]): ";
248  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
249  os.str() << "Empty index range is not allowed.");
250  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
251  os.str() << "Index range includes negative "
252  "index/ices, which is not allowed.");
253  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numCols, std::invalid_argument,
254  os.str() << "Index range exceeds number of "
255  "vectors " << numCols << " in the input "
256  "multivector.");
257  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
258  os.str() << "Should never get here!");
259  }
260  return mv.subView (index);
261  }
262 
264  { return static_cast<ptrdiff_t>(mv.getGlobalLength()); }
265 
267  { return mv.getNumVectors(); }
268 
270  { return mv.isConstantStride(); }
271 
273  const Teuchos::SerialDenseMatrix<int,BaseScalar>& B,
275  {
276  Teuchos::SerialDenseMatrix<int,Scalar> B_pce(B.numRows(), B.numCols());
277  for (int i=0; i<B.numRows(); i++)
278  for (int j=0; j<B.numCols(); j++)
279  B_pce(i,j) = B(i,j);
280  MvTimesMatAddMv(alpha, A, B_pce, beta, mv);
281  }
283  const Teuchos::SerialDenseMatrix<int,Scalar>& B,
285  {
286 #ifdef HAVE_BELOS_TPETRA_TIMERS
287  Teuchos::TimeMonitor lcltimer(*mvTimesMatAddMvTimer_);
288 #endif
289  // create local map
290  Teuchos::SerialComm<int> scomm;
291  Tpetra::Map<LO,GO,Node> LocalMap(B.numRows(), 0, Teuchos::rcpFromRef< const Teuchos::Comm<int> >(scomm), Tpetra::LocallyReplicated, A.getMap()->getNode());
292  // encapsulate Teuchos::SerialDenseMatrix data in ArrayView
293  Teuchos::ArrayView<const Scalar> Bvalues(B.values(),B.stride()*B.numCols());
294  // create locally replicated MultiVector with a copy of this data
295  Tpetra::MultiVector<Scalar,LO,GO,Node> B_mv(Teuchos::rcpFromRef(LocalMap),Bvalues,B.stride(),B.numCols());
296  // multiply
297  mv.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
298  }
299 
301  {
302  mv.update(alpha,A,beta,B,Teuchos::ScalarTraits<Scalar>::zero());
303  }
304 
306  { mv.scale(alpha); }
307 
308  static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<BaseScalar>& alphas )
309  {
310  std::vector<Scalar> alphas_pce(alphas.size());
311  for (int i=0; i<alphas.size(); i++) alphas_pce[i] = alphas[i];
312  mv.scale(alphas_pce);
313  }
314  static void MvScale ( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, const std::vector<Scalar>& alphas )
315  { mv.scale(alphas); }
316 
317  static void MvTransMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, Teuchos::SerialDenseMatrix<int,BaseScalar>& C)
318  {
319  Teuchos::SerialDenseMatrix<int,Scalar> C_pce(C.numRows(), C.numCols());
320  MvTransMv(alpha, A, B, C_pce);
321  for (int i=0; i<C.numRows(); i++)
322  for (int j=0; j<C.numCols(); j++)
323  C(i,j) = C_pce(i,j).coeff(0);
324  }
325  static void MvTransMv( Scalar alpha, const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, Teuchos::SerialDenseMatrix<int,Scalar>& C)
326  {
327 #ifdef HAVE_BELOS_TPETRA_TIMERS
328  Teuchos::TimeMonitor lcltimer(*mvTransMvTimer_);
329 #endif
330  // form alpha * A^H * B, then copy into SDM
331  // we will create a multivector C_mv from a a local map
332  // this map has a serial comm, the purpose being to short-circuit the MultiVector::reduce() call at the end of MultiVector::multiply()
333  // otherwise, the reduced multivector data would be copied back to the GPU, only to turn around and have to get it back here.
334  // this saves us a round trip for this data.
335  const int numRowsC = C.numRows(),
336  numColsC = C.numCols(),
337  strideC = C.stride();
338  Teuchos::SerialComm<int> scomm;
339  // create local map with serial comm
340  Tpetra::Map<LO,GO,Node> LocalMap(numRowsC, 0, Teuchos::rcpFromRef< const Teuchos::Comm<int> >(scomm), Tpetra::LocallyReplicated, A.getMap()->getNode());
341  // create local multivector to hold the result
342  const bool INIT_TO_ZERO = true;
343  Tpetra::MultiVector<Scalar,LO,GO,Node> C_mv(Teuchos::rcpFromRef(LocalMap),numColsC, INIT_TO_ZERO);
344  // multiply result into local multivector
345  C_mv.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,alpha,A,B,Teuchos::ScalarTraits<Scalar>::zero());
346  // get comm
347  Teuchos::RCP< const Teuchos::Comm<int> > pcomm = A.getMap()->getComm();
348  // create arrayview encapsulating the Teuchos::SerialDenseMatrix
349  Teuchos::ArrayView<Scalar> C_view(C.values(),strideC*numColsC);
350  if (pcomm->getSize() == 1) {
351  // no accumulation to do; simply extract the multivector data into C
352  // extract a copy of the result into the array view (and therefore, the SerialDenseMatrix)
353  C_mv.get1dCopy(C_view,strideC);
354  }
355  else {
356  // get a const host view of the data in C_mv
357  Teuchos::ArrayRCP<const Scalar> C_mv_view = C_mv.get1dView();
358  if (strideC == numRowsC) {
359  // sumall into C
360  Teuchos::reduceAll<int,Scalar>(*pcomm,Teuchos::REDUCE_SUM,numColsC*numRowsC,C_mv_view.getRawPtr(),C_view.getRawPtr());
361  }
362  else {
363  // sumall into temp, copy into C
364  Teuchos::Array<Scalar> destBuff(numColsC*numRowsC);
365  Teuchos::reduceAll<int,Scalar>(*pcomm,Teuchos::REDUCE_SUM,numColsC*numRowsC,C_mv_view.getRawPtr(),destBuff.getRawPtr());
366  for (int j=0; j < numColsC; ++j) {
367  for (int i=0; i < numRowsC; ++i) {
368  C_view[strideC*j+i] = destBuff[numRowsC*j+i];
369  }
370  }
371  }
372  }
373  }
374 
375  static void MvDot( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const Tpetra::MultiVector<Scalar,LO,GO,Node>& B, std::vector<BaseScalar> &dots)
376  {
377  TEUCHOS_TEST_FOR_EXCEPTION(A.getNumVectors() != B.getNumVectors(),std::invalid_argument,
378  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): A and B must have the same number of vectors.");
379 #ifdef HAVE_TPETRA_DEBUG
380  TEUCHOS_TEST_FOR_EXCEPTION(dots.size() < (typename std::vector<int>::size_type)A.getNumVectors(),std::invalid_argument,
381  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvDot(A,B,dots): dots must have room for all dot products.");
382 #endif
383  // Teuchos::ArrayView<Scalar> av(dots);
384  // A.dot(B,av(0,A.getNumVectors()));
385  Teuchos::Array<Scalar> pce_dots(A.getNumVectors());
386  A.dot(B, pce_dots);
387  for (unsigned int i=0; i<A.getNumVectors(); i++)
388  dots[i] = pce_dots[i].coeff(0);
389  }
390 
391  static void MvNorm(const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::vector<typename Teuchos::ScalarTraits<BaseScalar>::magnitudeType> &normvec, NormType type=TwoNorm)
392  {
393 #ifdef HAVE_TPETRA_DEBUG
394  TEUCHOS_TEST_FOR_EXCEPTION(normvec.size() < (typename std::vector<int>::size_type)mv.getNumVectors(),std::invalid_argument,
395  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::MvNorm(mv,normvec): normvec must have room for all norms.");
396 #endif
397  Teuchos::ArrayView<typename Teuchos::ScalarTraits<BaseScalar>::magnitudeType> av(normvec);
398  //Teuchos::Array<Scalar> pce_norms(mv.getNumVectors());
399  switch (type) {
400  case OneNorm:
401  mv.norm1(av(0,mv.getNumVectors()));
402  // mv.norm1(pce_norms);
403  // for (unsigned int i=0; i<mv.getNumVectors(); i++)
404  // normvec[i] = pce_norms[i].coeff(0);
405  break;
406  case TwoNorm:
407  mv.norm2(av(0,mv.getNumVectors()));
408  // mv.dot(mv, pce_norms);
409  // for (unsigned int i=0; i<mv.getNumVectors(); i++)
410  // normvec[i] = std::sqrt(pce_norms[i].coeff(0));
411  break;
412  case InfNorm:
413  mv.normInf(av(0,mv.getNumVectors()));
414  // mv.normInf(pce_norms);
415  // for (unsigned int i=0; i<mv.getNumVectors(); i++)
416  // normvec[i] = pce_norms[i].coeff(0);
417  break;
418  }
419 
420  }
421 
422  static void SetBlock( const Tpetra::MultiVector<Scalar,LO,GO,Node>& A, const std::vector<int>& index, Tpetra::MultiVector<Scalar,LO,GO,Node>& mv )
423  {
424 #ifdef HAVE_TPETRA_DEBUG
425  TEUCHOS_TEST_FOR_EXCEPTION((typename std::vector<int>::size_type)A.getNumVectors() < index.size(),std::invalid_argument,
426  "Belos::MultiVecTraits<Scalar,Tpetra::MultiVector>::SetBlock(A,index,mv): index must be the same size as A.");
427 #endif
428  Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > mvsub = CloneViewNonConst(mv,index);
429  if ((typename std::vector<int>::size_type)A.getNumVectors() > index.size()) {
430  Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > Asub = A.subView(Teuchos::Range1D(0,index.size()-1));
431  (*mvsub) = (*Asub);
432  }
433  else {
434  (*mvsub) = A;
435  }
436  mvsub = Teuchos::null;
437  }
438 
439  static void
441  const Teuchos::Range1D& index,
443  {
444  // Range1D bounds are signed; size_t is unsigned.
445  // Assignment of Tpetra::MultiVector is a deep copy.
446 
447  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
448  // fair to assume that the number of vectors won't overflow int,
449  // since the typical use case of multivectors involves few
450  // columns, but it's friendly to check just in case.
451  const size_t maxInt = static_cast<size_t> (Teuchos::OrdinalTraits<int>::max());
452  const bool overflow = maxInt < A.getNumVectors() && maxInt < mv.getNumVectors();
453  if (overflow)
454  {
455  std::ostringstream os;
456  os << "Belos::MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar, ..."
457  "> >::SetBlock(A, index=[" << index.lbound() << ", "
458  << index.ubound() << "], mv): ";
459  TEUCHOS_TEST_FOR_EXCEPTION(maxInt < A.getNumVectors(), std::range_error,
460  os.str() << "Number of columns in the input multi"
461  "vector 'A' (a size_t) overflows int.");
462  TEUCHOS_TEST_FOR_EXCEPTION(maxInt < mv.getNumVectors(), std::range_error,
463  os.str() << "Number of columns in the output multi"
464  "vector 'mv' (a size_t) overflows int.");
465  }
466  // We've already validated the static casts above.
467  const int numColsA = static_cast<int> (A.getNumVectors());
468  const int numColsMv = static_cast<int> (mv.getNumVectors());
469  // 'index' indexes into mv; it's the index set of the target.
470  const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
471  // We can't take more columns out of A than A has.
472  const bool validSource = index.size() <= numColsA;
473 
474  if (! validIndex || ! validSource)
475  {
476  std::ostringstream os;
477  os << "Belos::MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar, ..."
478  "> >::SetBlock(A, index=[" << index.lbound() << ", "
479  << index.ubound() << "], mv): ";
480  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
481  os.str() << "Range lower bound must be nonnegative.");
482  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
483  os.str() << "Range upper bound must be less than "
484  "the number of columns " << numColsA << " in the "
485  "'mv' output argument.");
486  TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
487  os.str() << "Range must have no more elements than"
488  " the number of columns " << numColsA << " in the "
489  "'A' input argument.");
490  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
491  }
492  typedef Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > MV_ptr;
493  typedef Teuchos::RCP<const Tpetra::MultiVector<Scalar,LO,GO,Node> > const_MV_ptr;
494 
495  // View of the relevant column(s) of the target multivector mv.
496  // We avoid view creation overhead by only creating a view if
497  // the index range is different than [0, (# columns in mv) - 1].
498  MV_ptr mv_view;
499  if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
500  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
501  else
502  mv_view = CloneViewNonConst (mv, index);
503 
504  // View of the relevant column(s) of the source multivector A.
505  // If A has fewer columns than mv_view, then create a view of
506  // the first index.size() columns of A.
507  const_MV_ptr A_view;
508  if (index.size() == numColsA)
509  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
510  else
511  A_view = CloneView (A, Teuchos::Range1D(0, index.size()-1));
512 
513  // Assignment of Tpetra::MultiVector objects via operator=()
514  // assumes that both arguments have compatible Maps. If
515  // HAVE_TPETRA_DEBUG is defined at compile time, operator=()
516  // will throw an std::runtime_error if the Maps are
517  // incompatible.
518  *mv_view = *A_view;
519  }
520 
521  static void
524  {
525  // Range1D bounds are signed; size_t is unsigned.
526  // Assignment of Tpetra::MultiVector is a deep copy.
527 
528  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
529  // fair to assume that the number of vectors won't overflow int,
530  // since the typical use case of multivectors involves few
531  // columns, but it's friendly to check just in case.
532  const size_t maxInt = static_cast<size_t> (Teuchos::OrdinalTraits<int>::max());
533  const bool overflow = maxInt < A.getNumVectors() && maxInt < mv.getNumVectors();
534  if (overflow)
535  {
536  std::ostringstream os;
537  os << "Belos::MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar, ..."
538  "> >::Assign(A, mv): ";
539  TEUCHOS_TEST_FOR_EXCEPTION(maxInt < A.getNumVectors(), std::range_error,
540  os.str() << "Number of columns in the input multi"
541  "vector 'A' (a size_t) overflows int.");
542  TEUCHOS_TEST_FOR_EXCEPTION(maxInt < mv.getNumVectors(), std::range_error,
543  os.str() << "Number of columns in the output multi"
544  "vector 'mv' (a size_t) overflows int.");
545  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
546  }
547  // We've already validated the static casts above.
548  const int numColsA = static_cast<int> (A.getNumVectors());
549  const int numColsMv = static_cast<int> (mv.getNumVectors());
550  if (numColsA > numColsMv)
551  {
552  std::ostringstream os;
553  os << "Belos::MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar, ..."
554  "> >::Assign(A, mv): ";
555  TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
556  os.str() << "Input multivector 'A' has "
557  << numColsA << " columns, but output multivector "
558  "'mv' has only " << numColsMv << " columns.");
559  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
560  }
561  // Assignment of Tpetra::MultiVector objects via operator=()
562  // assumes that both arguments have compatible Maps. If
563  // HAVE_TPETRA_DEBUG is defined at compile time, operator=()
564  // will throw an std::runtime_error if the Maps are
565  // incompatible.
566  if (numColsA == numColsMv)
567  mv = A;
568  else
569  {
570  Teuchos::RCP<Tpetra::MultiVector<Scalar,LO,GO,Node> > mv_view =
571  CloneViewNonConst (mv, Teuchos::Range1D(0, numColsA-1));
572  *mv_view = A;
573  }
574  }
575 
576 
578  {
579  mv.randomize();
580  }
581 
582  static void MvInit( Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero() )
583  { mv.putScalar(alpha); }
584 
585  static void MvPrint( const Tpetra::MultiVector<Scalar,LO,GO,Node>& mv, std::ostream& os )
586  {
587  Teuchos::FancyOStream fos(Teuchos::rcp(&os,false));
588  mv.describe(fos,Teuchos::VERB_EXTREME);
589  }
590 
591 #ifdef HAVE_BELOS_TSQR
592  typedef Tpetra::TsqrAdaptor< Tpetra::MultiVector< Scalar, LO, GO, Node > > tsqr_adaptor_type;
596 #endif // HAVE_BELOS_TSQR
597  };
598 
600  //
601  // Implementation of the Belos::OperatorTraits for Tpetra::Operator.
602  //
604 
606  template <class BaseScalar, class Storage, class LO, class GO, class Node>
607  class OperatorTraits <BaseScalar, Tpetra::MultiVector<Sacado::PCE::OrthogPoly<BaseScalar,Storage>,LO,GO,Node>, Tpetra::Operator<Sacado::PCE::OrthogPoly<BaseScalar,Storage>,LO,GO,Node> >
608  {
609  public:
611  static void
612  Apply (const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
615  ETrans trans=NOTRANS)
616  {
617  switch (trans) {
618  case NOTRANS:
619  Op.apply(X,Y,Teuchos::NO_TRANS);
620  break;
621  case TRANS:
622  Op.apply(X,Y,Teuchos::TRANS);
623  break;
624  case CONJTRANS:
625  Op.apply(X,Y,Teuchos::CONJ_TRANS);
626  break;
627  default:
628  const std::string scalarName = Teuchos::TypeNameTraits<Scalar>::name();
629  const std::string loName = Teuchos::TypeNameTraits<LO>::name();
630  const std::string goName = Teuchos::TypeNameTraits<GO>::name();
631  const std::string nodeName = Teuchos::TypeNameTraits<Node>::name();
632  const std::string otName = "Belos::OperatorTraits<" + scalarName
633  + "," + loName + "," + goName + "," + nodeName + ">";
634  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, otName << ": Should never "
635  "get here; fell through a switch statement. "
636  "Please report this bug to the Belos developers.");
637  }
638  }
639 
640  static bool
641  HasApplyTranspose (const Tpetra::Operator<Scalar,LO,GO,Node>& Op)
642  {
643  return Op.hasTransposeApply ();
644  }
645  };
646 
647 } // end of Belos namespace
648 
649 #endif
650 // end of file BELOS_PCE_TPETRA_ADAPTER_HPP
static void MvInit(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, Scalar alpha=Teuchos::ScalarTraits< Scalar >::zero())
static Teuchos::RCP< Tpetra::MultiVector< Scalar, LO, GO, Node > > CloneCopy(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< int > &index)
static void Assign(const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, Tpetra::MultiVector< Scalar, LO, GO, Node > &mv)
static void MvTransMv(Scalar alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, Teuchos::SerialDenseMatrix< int, BaseScalar > &C)
static void MvTimesMatAddMv(Scalar alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Teuchos::SerialDenseMatrix< int, Scalar > &B, Scalar beta, Tpetra::MultiVector< Scalar, LO, GO, Node > &mv)
static void SetBlock(const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const std::vector< int > &index, Tpetra::MultiVector< Scalar, LO, GO, Node > &mv)
static Teuchos::RCP< const Tpetra::MultiVector< Scalar, LO, GO, Node > > CloneView(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const Teuchos::Range1D &index)
static Teuchos::RCP< Tpetra::MultiVector< Scalar, LO, GO, Node > > CloneCopy(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv)
static Teuchos::RCP< Tpetra::MultiVector< Scalar, LO, GO, Node > > Clone(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const int numvecs)
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< BaseScalar > &alphas)
KOKKOS_INLINE_FUNCTION PCE< Storage > max(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
static void MvTransMv(Scalar alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, Teuchos::SerialDenseMatrix< int, Scalar > &C)
static void MvTimesMatAddMv(Scalar alpha, const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Teuchos::SerialDenseMatrix< int, BaseScalar > &B, Scalar beta, Tpetra::MultiVector< Scalar, LO, GO, Node > &mv)
static void SetBlock(const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Teuchos::Range1D &index, Tpetra::MultiVector< Scalar, LO, GO, Node > &mv)
expr expr expr expr j
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 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)
static Teuchos::RCP< Tpetra::MultiVector< Scalar, LO, GO, Node > > CloneViewNonConst(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< int > &index)
static Teuchos::RCP< Tpetra::MultiVector< Scalar, LO, GO, Node > > CloneCopy(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const Teuchos::Range1D &index)
static void MvNorm(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, std::vector< typename Teuchos::ScalarTraits< BaseScalar >::magnitudeType > &normvec, NormType type=TwoNorm)
static void MvDot(const Tpetra::MultiVector< Scalar, LO, GO, Node > &A, const Tpetra::MultiVector< Scalar, LO, GO, Node > &B, std::vector< BaseScalar > &dots)
static Teuchos::RCP< const Tpetra::MultiVector< Scalar, LO, GO, Node > > CloneView(const Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< int > &index)
static Teuchos::RCP< Tpetra::MultiVector< Scalar, LO, GO, Node > > CloneViewNonConst(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const Teuchos::Range1D &index)
static void MvScale(Tpetra::MultiVector< Scalar, LO, GO, Node > &mv, const std::vector< Scalar > &alphas)