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