Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
BelosThyraAdapter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Stratimikos: Thyra-based strategies for linear solvers
5 // Copyright (2006) 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 Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
52 #ifndef BELOS_THYRA_ADAPTER_HPP
53 #define BELOS_THYRA_ADAPTER_HPP
54 
55 #include "BelosConfigDefs.hpp"
56 #include "BelosMultiVecTraits.hpp"
57 #include "BelosOperatorTraits.hpp"
58 
59 #include <Thyra_DetachedMultiVectorView.hpp>
60 #include <Thyra_MultiVectorBase.hpp>
61 #include <Thyra_MultiVectorStdOps.hpp>
62 #ifdef HAVE_BELOS_TSQR
63 # include <Thyra_TsqrAdaptor.hpp>
64 #endif // HAVE_BELOS_TSQR
65 
66 namespace Belos {
67 
69  //
70  // Implementation of the Belos::MultiVecTraits for Thyra::MultiVectorBase
71  //
73 
80  template<class ScalarType>
81  class MultiVecTraits< ScalarType, Thyra::MultiVectorBase<ScalarType> >
82  {
83  private:
84  typedef Thyra::MultiVectorBase<ScalarType> TMVB;
85  typedef Teuchos::ScalarTraits<ScalarType> ST;
86  typedef typename ST::magnitudeType magType;
87 
88  public:
89 
92 
97  static Teuchos::RCP<TMVB> Clone( const TMVB& mv, const int numvecs )
98  {
99  Teuchos::RCP<TMVB> c = Thyra::createMembers( mv.range(), numvecs );
100  return c;
101  }
102 
107  static Teuchos::RCP<TMVB> CloneCopy( const TMVB& mv )
108  {
109  int numvecs = mv.domain()->dim();
110  // create the new multivector
111  Teuchos::RCP< TMVB > cc = Thyra::createMembers( mv.range(), numvecs );
112  // copy the data from the source multivector to the new multivector
113  Thyra::assign(cc.ptr(), mv);
114  return cc;
115  }
116 
122  static Teuchos::RCP<TMVB> CloneCopy( const TMVB& mv, const std::vector<int>& index )
123  {
124  int numvecs = index.size();
125  // create the new multivector
126  Teuchos::RCP<TMVB> cc = Thyra::createMembers( mv.range(), numvecs );
127  // create a view to the relevant part of the source multivector
128  Teuchos::RCP<const TMVB> view = mv.subView(index);
129  // copy the data from the relevant view to the new multivector
130  Thyra::assign(cc.ptr(), *view);
131  return cc;
132  }
133 
134  static Teuchos::RCP<TMVB>
135  CloneCopy (const TMVB& mv, const Teuchos::Range1D& index)
136  {
137  const int numVecs = index.size();
138  // Create the new multivector
139  Teuchos::RCP<TMVB> cc = Thyra::createMembers (mv.range(), numVecs);
140  // Create a view to the relevant part of the source multivector
141  Teuchos::RCP<const TMVB> view = mv.subView (index);
142  // Copy the data from the view to the new multivector.
143  Thyra::assign (cc.ptr(), *view);
144  return cc;
145  }
146 
152  static Teuchos::RCP<TMVB> CloneViewNonConst( TMVB& mv, const std::vector<int>& index )
153  {
154  int numvecs = index.size();
155 
156  // We do not assume that the indices are sorted, nor do we check that
157  // index.size() > 0. This code is fail-safe, in the sense that a zero
158  // length index std::vector will pass the error on the Thyra.
159 
160  // Thyra has two ways to create an indexed View:
161  // * contiguous (via a range of columns)
162  // * indexed (via a std::vector of column indices)
163  // The former is significantly more efficient than the latter, in terms of
164  // computations performed with/against the created view.
165  // We will therefore check to see if the given indices are contiguous, and
166  // if so, we will use the contiguous view creation method.
167 
168  int lb = index[0];
169  bool contig = true;
170  for (int i=0; i<numvecs; i++) {
171  if (lb+i != index[i]) contig = false;
172  }
173 
174  Teuchos::RCP< TMVB > cc;
175  if (contig) {
176  const Thyra::Range1D rng(lb,lb+numvecs-1);
177  // create a contiguous view to the relevant part of the source multivector
178  cc = mv.subView(rng);
179  }
180  else {
181  // create an indexed view to the relevant part of the source multivector
182  cc = mv.subView(index);
183  }
184  return cc;
185  }
186 
187  static Teuchos::RCP<TMVB>
188  CloneViewNonConst (TMVB& mv, const Teuchos::Range1D& index)
189  {
190  // We let Thyra be responsible for checking that the index range
191  // is nonempty.
192  //
193  // Create and return a contiguous view to the relevant part of
194  // the source multivector.
195  return mv.subView (index);
196  }
197 
198 
204  static Teuchos::RCP<const TMVB> CloneView( const TMVB& mv, const std::vector<int>& index )
205  {
206  int numvecs = index.size();
207 
208  // We do not assume that the indices are sorted, nor do we check that
209  // index.size() > 0. This code is fail-safe, in the sense that a zero
210  // length index std::vector will pass the error on the Thyra.
211 
212  // Thyra has two ways to create an indexed View:
213  // * contiguous (via a range of columns)
214  // * indexed (via a std::vector of column indices)
215  // The former is significantly more efficient than the latter, in terms of
216  // computations performed with/against the created view.
217  // We will therefore check to see if the given indices are contiguous, and
218  // if so, we will use the contiguous view creation method.
219 
220  int lb = index[0];
221  bool contig = true;
222  for (int i=0; i<numvecs; i++) {
223  if (lb+i != index[i]) contig = false;
224  }
225 
226  Teuchos::RCP< const TMVB > cc;
227  if (contig) {
228  const Thyra::Range1D rng(lb,lb+numvecs-1);
229  // create a contiguous view to the relevant part of the source multivector
230  cc = mv.subView(rng);
231  }
232  else {
233  // create an indexed view to the relevant part of the source multivector
234  cc = mv.subView(index);
235  }
236  return cc;
237  }
238 
239  static Teuchos::RCP<const TMVB>
240  CloneView (const TMVB& mv, const Teuchos::Range1D& index)
241  {
242  // We let Thyra be responsible for checking that the index range
243  // is nonempty.
244  //
245  // Create and return a contiguous view to the relevant part of
246  // the source multivector.
247  return mv.subView (index);
248  }
249 
251 
254 
256  static ptrdiff_t GetGlobalLength( const TMVB& mv ) {
257  return Teuchos::as<ptrdiff_t>(mv.range()->dim());
258  }
259 
261  static int GetNumberVecs( const TMVB& mv )
262  { return mv.domain()->dim(); }
263 
265 
268 
271  static void MvTimesMatAddMv( const ScalarType alpha, const TMVB& A,
272  const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
273  const ScalarType beta, TMVB& mv )
274  {
275  using Teuchos::arrayView; using Teuchos::arcpFromArrayView;
276  const int m = B.numRows();
277  const int n = B.numCols();
278  // Create a view of the B object!
279  Teuchos::RCP< const TMVB >
280  B_thyra = Thyra::createMembersView(
281  A.domain(),
282  RTOpPack::ConstSubMultiVectorView<ScalarType>(
283  0, m, 0, n,
284  arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
285  )
286  );
287  // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
288  Thyra::apply<ScalarType>(A, Thyra::NOTRANS, *B_thyra, Teuchos::outArg(mv), alpha, beta);
289  }
290 
293  static void MvAddMv( const ScalarType alpha, const TMVB& A,
294  const ScalarType beta, const TMVB& B, TMVB& mv )
295  {
296  typedef Teuchos::ScalarTraits<ScalarType> ST;
297  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::inoutArg;
298 
299  Thyra::linear_combination<ScalarType>(
300  tuple(alpha, beta)(), tuple(ptrInArg(A), ptrInArg(B))(), ST::zero(), inoutArg(mv));
301  }
302 
305  static void MvScale ( TMVB& mv, const ScalarType alpha )
306  { Thyra::scale(alpha, Teuchos::inoutArg(mv)); }
307 
310  static void MvScale (TMVB& mv, const std::vector<ScalarType>& alpha)
311  {
312  for (unsigned int i=0; i<alpha.size(); i++) {
313  Thyra::scale<ScalarType> (alpha[i], mv.col(i).ptr());
314  }
315  }
316 
319  static void MvTransMv( const ScalarType alpha, const TMVB& A, const TMVB& mv,
320  Teuchos::SerialDenseMatrix<int,ScalarType>& B )
321  {
322  using Teuchos::arrayView; using Teuchos::arcpFromArrayView;
323  // Create a multivector to hold the result (m by n)
324  int m = A.domain()->dim();
325  int n = mv.domain()->dim();
326  // Create a view of the B object!
327  Teuchos::RCP< TMVB >
328  B_thyra = Thyra::createMembersView(
329  A.domain(),
330  RTOpPack::SubMultiVectorView<ScalarType>(
331  0, m, 0, n,
332  arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
333  )
334  );
335  Thyra::apply<ScalarType>(A, Thyra::CONJTRANS, mv, B_thyra.ptr(), alpha);
336  }
337 
341  static void MvDot( const TMVB& mv, const TMVB& A, std::vector<ScalarType>& b )
342  { Thyra::dots(mv, A, Teuchos::arrayViewFromVector(b)); }
343 
345 
348 
352  static void MvNorm( const TMVB& mv, std::vector<magType>& normvec,
353  NormType type = TwoNorm ) {
354  if(type == TwoNorm)
355  Thyra::norms_2(mv, Teuchos::arrayViewFromVector(normvec));
356  else if(type == OneNorm)
357  Thyra::norms_1(mv, Teuchos::arrayViewFromVector(normvec));
358  else if(type == InfNorm)
359  Thyra::norms_inf(mv, Teuchos::arrayViewFromVector(normvec));
360  else
361  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
362  "Belos::MultiVecTraits::MvNorm (Thyra specialization): "
363  "invalid norm type. Must be either TwoNorm, OneNorm or InfNorm");
364  }
365 
367 
370 
373  static void SetBlock( const TMVB& A, const std::vector<int>& index, TMVB& mv )
374  {
375  // Extract the "numvecs" columns of mv indicated by the index std::vector.
376  int numvecs = index.size();
377  std::vector<int> indexA(numvecs);
378  int numAcols = A.domain()->dim();
379  for (int i=0; i<numvecs; i++) {
380  indexA[i] = i;
381  }
382  // Thyra::assign requires that both arguments have the same number of
383  // vectors. Enforce this, by shrinking one to match the other.
384  if ( numAcols < numvecs ) {
385  // A does not have enough columns to satisfy index_plus. Shrink
386  // index_plus.
387  numvecs = numAcols;
388  }
389  else if ( numAcols > numvecs ) {
390  numAcols = numvecs;
391  indexA.resize( numAcols );
392  }
393  // create a view to the relevant part of the source multivector
394  Teuchos::RCP< const TMVB > relsource = A.subView(indexA);
395  // create a view to the relevant part of the destination multivector
396  Teuchos::RCP< TMVB > reldest = mv.subView(index);
397  // copy the data to the destination multivector subview
398  Thyra::assign(reldest.ptr(), *relsource);
399  }
400 
401  static void
402  SetBlock (const TMVB& A, const Teuchos::Range1D& index, TMVB& mv)
403  {
404  const int numColsA = A.domain()->dim();
405  const int numColsMv = mv.domain()->dim();
406  // 'index' indexes into mv; it's the index set of the target.
407  const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
408  // We can't take more columns out of A than A has.
409  const bool validSource = index.size() <= numColsA;
410 
411  if (! validIndex || ! validSource)
412  {
413  std::ostringstream os;
414  os << "Belos::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar> "
415  ">::SetBlock(A, [" << index.lbound() << ", " << index.ubound()
416  << "], mv): ";
417  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
418  os.str() << "Range lower bound must be nonnegative.");
419  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
420  os.str() << "Range upper bound must be less than "
421  "the number of columns " << numColsA << " in the "
422  "'mv' output argument.");
423  TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
424  os.str() << "Range must have no more elements than"
425  " the number of columns " << numColsA << " in the "
426  "'A' input argument.");
427  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
428  }
429 
430  // View of the relevant column(s) of the target multivector mv.
431  // We avoid view creation overhead by only creating a view if
432  // the index range is different than [0, (# columns in mv) - 1].
433  Teuchos::RCP<TMVB> mv_view;
434  if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
435  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
436  else
437  mv_view = mv.subView (index);
438 
439  // View of the relevant column(s) of the source multivector A.
440  // If A has fewer columns than mv_view, then create a view of
441  // the first index.size() columns of A.
442  Teuchos::RCP<const TMVB> A_view;
443  if (index.size() == numColsA)
444  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
445  else
446  A_view = A.subView (Teuchos::Range1D(0, index.size()-1));
447 
448  // Copy the data to the destination multivector.
449  Thyra::assign(mv_view.ptr(), *A_view);
450  }
451 
452  static void
453  Assign (const TMVB& A, TMVB& mv)
454  {
455  const int numColsA = A.domain()->dim();
456  const int numColsMv = mv.domain()->dim();
457  if (numColsA > numColsMv)
458  {
459  std::ostringstream os;
460  os << "Belos::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar>"
461  " >::Assign(A, mv): ";
462  TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
463  os.str() << "Input multivector 'A' has "
464  << numColsA << " columns, but output multivector "
465  "'mv' has only " << numColsMv << " columns.");
466  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
467  }
468  // Copy the data to the destination multivector.
469  if (numColsA == numColsMv) {
470  Thyra::assign (Teuchos::outArg (mv), A);
471  } else {
472  Teuchos::RCP<TMVB> mv_view =
473  CloneViewNonConst (mv, Teuchos::Range1D(0, numColsA-1));
474  Thyra::assign (mv_view.ptr(), A);
475  }
476  }
477 
480  static void MvRandom( TMVB& mv )
481  {
482  // Thyra::randomize generates via a uniform distribution on [l,u]
483  // We will use this to generate on [-1,1]
484  Thyra::randomize<ScalarType>(
485  -Teuchos::ScalarTraits<ScalarType>::one(),
486  Teuchos::ScalarTraits<ScalarType>::one(),
487  Teuchos::outArg(mv));
488  }
489 
491  static void
492  MvInit (TMVB& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero())
493  {
494  Thyra::assign (Teuchos::outArg (mv), alpha);
495  }
496 
498 
501 
504  static void MvPrint( const TMVB& mv, std::ostream& os )
505  { os << describe(mv,Teuchos::VERB_EXTREME); }
506 
508 
509 #ifdef HAVE_BELOS_TSQR
510  typedef Thyra::TsqrAdaptor< ScalarType > tsqr_adaptor_type;
516 #endif // HAVE_BELOS_TSQR
517  };
518 
520  //
521  // Implementation of the Belos::OperatorTraits for Thyra::LinearOpBase
522  //
524 
532  template<class ScalarType>
533  class OperatorTraits <ScalarType,
534  Thyra::MultiVectorBase<ScalarType>,
535  Thyra::LinearOpBase<ScalarType> >
536  {
537  private:
538  typedef Thyra::MultiVectorBase<ScalarType> TMVB;
539  typedef Thyra::LinearOpBase<ScalarType> TLOB;
540 
541  public:
557  static void
558  Apply (const TLOB& Op,
559  const TMVB& x,
560  TMVB& y,
561  ETrans trans = NOTRANS)
562  {
563  Thyra::EOpTransp whichOp;
564 
565  // We don't check here whether the operator implements the
566  // requested operation. Call HasApplyTranspose() to check.
567  // Thyra::LinearOpBase implementations are not required to
568  // implement NOTRANS. However, Belos needs NOTRANS
569  // (obviously!), so we assume that Op implements NOTRANS.
570  if (trans == NOTRANS)
571  whichOp = Thyra::NOTRANS;
572  else if (trans == TRANS)
573  whichOp = Thyra::TRANS;
574  else if (trans == CONJTRANS)
575  whichOp = Thyra::CONJTRANS;
576  else
577  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
578  "Belos::OperatorTraits::Apply (Thyra specialization): "
579  "'trans' argument must be neither NOTRANS=" << NOTRANS
580  << ", TRANS=" << TRANS << ", or CONJTRANS=" << CONJTRANS
581  << ", but instead has an invalid value of " << trans << ".");
582  Thyra::apply<ScalarType>(Op, whichOp, x, Teuchos::outArg(y));
583  }
584 
586  static bool HasApplyTranspose (const TLOB& Op)
587  {
588  typedef Teuchos::ScalarTraits<ScalarType> STS;
589 
590  // Thyra::LinearOpBase's interface lets you check whether the
591  // operator implements any of all four possible combinations of
592  // conjugation and transpose. Belos only needs transpose
593  // (TRANS) if the operator is real; in that case, Apply() does
594  // the same thing with trans = CONJTRANS or TRANS. If the
595  // operator is complex, Belos needs both transpose and conjugate
596  // transpose (CONJTRANS) if the operator is complex.
597  return Op.opSupported (Thyra::TRANS) &&
598  (! STS::isComplex || Op.opSupported (Thyra::CONJTRANS));
599  }
600  };
601 
602 } // end of Belos namespace
603 
604 #endif
605 // end of file BELOS_THYRA_ADAPTER_HPP
static int GetNumberVecs(const TMVB &mv)
Obtain the number of vectors in mv.
static ptrdiff_t GetGlobalLength(const TMVB &mv)
Obtain the std::vector length of mv.
static void MvPrint(const TMVB &mv, std::ostream &os)
Print the mv multi-std::vector to the os output stream.
Stub adaptor from Thyra::MultiVectorBase to TSQR.
double ST
static Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv, const std::vector< int > &index)
Creates a new MultiVectorBase and copies the selected contents of mv into the new std::vector (deep c...
static void MvTransMv(const ScalarType alpha, const TMVB &A, const TMVB &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
static void MvScale(TMVB &mv, const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
static void MvDot(const TMVB &mv, const TMVB &A, std::vector< ScalarType > &b)
Compute a std::vector b where the components are the individual dot-products of the i-th columns of A...
static Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv)
Creates a new MultiVectorBase and copies contents of mv into the new std::vector (deep copy)...
static void MvScale(TMVB &mv, const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
static void MvTimesMatAddMv(const ScalarType alpha, const TMVB &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, TMVB &mv)
Update mv with .
static Teuchos::RCP< const TMVB > CloneView(const TMVB &mv, const std::vector< int > &index)
Creates a new const MultiVectorBase that shares the selected contents of mv (shallow copy)...
static Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv, const Teuchos::Range1D &index)
static void Apply(const TLOB &Op, const TMVB &x, TMVB &y, ETrans trans=NOTRANS)
Apply Op to x, storing the result in y.
static void MvRandom(TMVB &mv)
Replace the vectors in mv with random vectors.
static Teuchos::RCP< TMVB > CloneViewNonConst(TMVB &mv, const Teuchos::Range1D &index)
static void SetBlock(const TMVB &A, const std::vector< int > &index, TMVB &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static void SetBlock(const TMVB &A, const Teuchos::Range1D &index, TMVB &mv)
static void MvNorm(const TMVB &mv, std::vector< magType > &normvec, NormType type=TwoNorm)
Compute the 2-norm of each individual std::vector of mv. Upon return, normvec[i] holds the value of ...
static void MvInit(TMVB &mv, ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static bool HasApplyTranspose(const TLOB &Op)
Whether the operator implements applying the transpose.
static Teuchos::RCP< TMVB > CloneViewNonConst(TMVB &mv, const std::vector< int > &index)
Creates a new MultiVectorBase that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< const TMVB > CloneView(const TMVB &mv, const Teuchos::Range1D &index)
static Teuchos::RCP< TMVB > Clone(const TMVB &mv, const int numvecs)
Creates a new empty MultiVectorBase containing numvecs columns.
static void MvAddMv(const ScalarType alpha, const TMVB &A, const ScalarType beta, const TMVB &B, TMVB &mv)
Replace mv with .