Anasazi  Version of the Day
AnasaziThyraAdapter.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Anasazi: Block Eigensolvers Package
5 // Copyright (2004) 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 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
34 #ifndef ANASAZI_THYRA_ADAPTER_HPP
35 #define ANASAZI_THYRA_ADAPTER_HPP
36 
39 #include "AnasaziConfigDefs.hpp"
40 
41 #include <Thyra_DetachedMultiVectorView.hpp>
42 #include <Thyra_MultiVectorBase.hpp>
43 #include <Thyra_MultiVectorStdOps.hpp>
44 #include <Thyra_VectorStdOps.hpp>
45 
46 #include <Teuchos_Ptr.hpp>
47 #include <Teuchos_ArrayRCP.hpp>
48 #include <Teuchos_ArrayView.hpp>
49 
50 namespace Anasazi {
51 
53  //
54  // Implementation of the Anasazi::MultiVecTraits for Thyra::MultiVectorBase
55  //
57 
66  template<class ScalarType>
67  class MultiVecTraits< ScalarType, Thyra::MultiVectorBase<ScalarType> >
68  {
69  private:
70  typedef Thyra::MultiVectorBase<ScalarType> TMVB;
71  typedef Teuchos::ScalarTraits<ScalarType> ST;
72  typedef typename ST::magnitudeType magType;
73 
74  public:
75 
78 
83  static Teuchos::RCP<TMVB> Clone( const TMVB & mv, const int numvecs )
84  {
85  Teuchos::RCP<TMVB> c = Thyra::createMembers( mv.range(), numvecs );
86  return c;
87  }
88 
93  static Teuchos::RCP<TMVB>
94  CloneCopy (const TMVB& mv)
95  {
96  const int numvecs = mv.domain()->dim();
97  // create the new multivector
98  Teuchos::RCP< TMVB > cc = Thyra::createMembers (mv.range(), numvecs);
99  // copy the data from the source multivector to the new multivector
100  Thyra::assign (Teuchos::outArg (*cc), mv);
101  return cc;
102  }
103 
109  static Teuchos::RCP< TMVB > CloneCopy( const TMVB & mv, const std::vector<int>& index )
110  {
111  const int numvecs = index.size();
112  // create the new multivector
113  Teuchos::RCP<TMVB > cc = Thyra::createMembers (mv.range(), numvecs);
114  // create a view to the relevant part of the source multivector
115  Teuchos::RCP<const TMVB > view = mv.subView ( Teuchos::arrayViewFromVector( index ) );
116  // copy the data from the relevant view to the new multivector
117  Thyra::assign (Teuchos::outArg (*cc), *view);
118  return cc;
119  }
120 
121  static Teuchos::RCP<TMVB>
122  CloneCopy (const TMVB& mv, const Teuchos::Range1D& index)
123  {
124  const 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 view to the new multivector.
130  Thyra::assign (Teuchos::outArg (*cc), *view);
131  return cc;
132  }
133 
139  static Teuchos::RCP< TMVB > CloneViewNonConst( TMVB & mv, const std::vector<int>& index )
140  {
141  int numvecs = index.size();
142 
143  // We do not assume that the indices are sorted, nor do we check that
144  // index.size() > 0. This code is fail-safe, in the sense that a zero
145  // length index vector will pass the error on the Thyra.
146 
147  // Thyra has two ways to create an indexed View:
148  // * contiguous (via a range of columns)
149  // * indexed (via a vector of column indices)
150  // The former is significantly more efficient than the latter, in terms of
151  // computations performed with/against the created view.
152  // We will therefore check to see if the given indices are contiguous, and
153  // if so, we will use the contiguous view creation method.
154 
155  int lb = index[0];
156  bool contig = true;
157  for (int i=0; i<numvecs; i++) {
158  if (lb+i != index[i]) contig = false;
159  }
160 
161  Teuchos::RCP< TMVB > cc;
162  if (contig) {
163  const Thyra::Range1D rng(lb,lb+numvecs-1);
164  // create a contiguous view to the relevant part of the source multivector
165  cc = mv.subView(rng);
166  }
167  else {
168  // create an indexed view to the relevant part of the source multivector
169  cc = mv.subView( Teuchos::arrayViewFromVector( index ) );
170  }
171  return cc;
172  }
173 
174  static Teuchos::RCP<TMVB>
175  CloneViewNonConst (TMVB& mv, const Teuchos::Range1D& index)
176  {
177  // We let Thyra be responsible for checking that the index range
178  // is nonempty.
179  //
180  // Create and return a contiguous view to the relevant part of
181  // the source multivector.
182  return mv.subView (index);
183  }
184 
190  static Teuchos::RCP<const TMVB > CloneView( const TMVB & mv, const std::vector<int>& index )
191  {
192  int numvecs = index.size();
193 
194  // We do not assume that the indices are sorted, nor do we check that
195  // index.size() > 0. This code is fail-safe, in the sense that a zero
196  // length index vector will pass the error on the Thyra.
197 
198  // Thyra has two ways to create an indexed View:
199  // * contiguous (via a range of columns)
200  // * indexed (via a vector of column indices)
201  // The former is significantly more efficient than the latter, in terms of
202  // computations performed with/against the created view.
203  // We will therefore check to see if the given indices are contiguous, and
204  // if so, we will use the contiguous view creation method.
205 
206  int lb = index[0];
207  bool contig = true;
208  for (int i=0; i<numvecs; i++) {
209  if (lb+i != index[i]) contig = false;
210  }
211 
212  Teuchos::RCP< const TMVB > cc;
213  if (contig) {
214  const Thyra::Range1D rng(lb,lb+numvecs-1);
215  // create a contiguous view to the relevant part of the source multivector
216  cc = mv.subView(rng);
217  }
218  else {
219  // create an indexed view to the relevant part of the source multivector
220  cc = mv.subView(Teuchos::arrayViewFromVector( index ) );
221  }
222  return cc;
223  }
224 
225  static Teuchos::RCP<const TMVB>
226  CloneView (const TMVB& mv, const Teuchos::Range1D& index)
227  {
228  // We let Thyra be responsible for checking that the index range
229  // is nonempty.
230  //
231  // Create and return a contiguous view to the relevant part of
232  // the source multivector.
233  return mv.subView (index);
234  }
235 
237 
240 
242  static ptrdiff_t GetGlobalLength( const TMVB & mv )
243  { return mv.range()->dim(); }
244 
246  static int GetNumberVecs( const TMVB & mv )
247  { return mv.domain()->dim(); }
248 
250 
253 
256  static void MvTimesMatAddMv( const ScalarType alpha, const TMVB & A,
257  const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
258  const ScalarType beta, TMVB & mv )
259  {
260  int m = B.numRows();
261  int n = B.numCols();
262  // Create a view of the B object!
263  Teuchos::RCP< const TMVB >
264  B_thyra = Thyra::createMembersView(
265  A.domain(),
266  RTOpPack::ConstSubMultiVectorView<ScalarType>(
267  0, m, 0, n,
268  Teuchos::arcpFromArrayView(Teuchos::arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
269  )
270  );
271  // perform the operation via A: mv <- alpha*A*B_thyra + beta*mv
272  A.apply(Thyra::NOTRANS,*B_thyra,Teuchos::outArg (mv),alpha,beta);
273  }
274 
277  static void MvAddMv( const ScalarType alpha, const TMVB & A,
278  const ScalarType beta, const TMVB & B, TMVB & mv )
279  {
280  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::inoutArg;
281 
282  Thyra::linear_combination<ScalarType>(
283  tuple(alpha, beta)(), tuple(ptrInArg(A), ptrInArg(B))(), ST::zero(), inoutArg(mv));
284  }
285 
288  static void MvTransMv( const ScalarType alpha, const TMVB & A, const TMVB & mv,
289  Teuchos::SerialDenseMatrix<int,ScalarType>& B )
290  {
291  // Create a multivector to hold the result (m by n)
292  int m = A.domain()->dim();
293  int n = mv.domain()->dim();
294  // Create a view of the B object!
295  Teuchos::RCP< TMVB >
296  B_thyra = Thyra::createMembersView(
297  A.domain(),
298  RTOpPack::SubMultiVectorView<ScalarType>(
299  0, m, 0, n,
300  Teuchos::arcpFromArrayView(Teuchos::arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
301  )
302  );
303  A.apply(Thyra::CONJTRANS,mv,B_thyra.ptr(),alpha,Teuchos::ScalarTraits<ScalarType>::zero());
304  }
305 
309  static void MvDot( const TMVB & mv, const TMVB & A, std::vector<ScalarType> &b )
310  { Thyra::dots(mv,A,Teuchos::arrayViewFromVector( b )); }
311 
314  static void
315  MvScale (TMVB& mv,
316  const ScalarType alpha)
317  {
318  Thyra::scale (alpha, Teuchos::inOutArg (mv));
319  }
320 
323  static void
324  MvScale (TMVB& mv,
325  const std::vector<ScalarType>& alpha)
326  {
327  for (unsigned int i=0; i<alpha.size(); i++) {
328  Thyra::scale (alpha[i], mv.col(i).ptr());
329  }
330  }
331 
333 
336 
340  static void MvNorm( const TMVB & mv, std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> &normvec )
341  { Thyra::norms_2(mv,Teuchos::arrayViewFromVector( normvec )); }
342 
344 
347 
350  static void SetBlock( const TMVB & A, const std::vector<int>& index, TMVB & mv )
351  {
352  // Extract the "numvecs" columns of mv indicated by the index vector.
353  int numvecs = index.size();
354  std::vector<int> indexA(numvecs);
355  int numAcols = A.domain()->dim();
356  for (int i=0; i<numvecs; i++) {
357  indexA[i] = i;
358  }
359  // Thyra::assign requires that both arguments have the same number of
360  // vectors. Enforce this, by shrinking one to match the other.
361  if ( numAcols < numvecs ) {
362  // A does not have enough columns to satisfy index_plus. Shrink
363  // index_plus.
364  numvecs = numAcols;
365  }
366  else if ( numAcols > numvecs ) {
367  numAcols = numvecs;
368  indexA.resize( numAcols );
369  }
370  // create a view to the relevant part of the source multivector
371  Teuchos::RCP< const TMVB > relsource = A.subView( Teuchos::arrayViewFromVector( indexA ) );
372  // create a view to the relevant part of the destination multivector
373  Teuchos::RCP< TMVB > reldest = mv.subView( Teuchos::arrayViewFromVector( index ) );
374  // copy the data to the destination multivector subview
375  Thyra::assign (Teuchos::outArg (*reldest), *relsource);
376  }
377 
378  static void
379  SetBlock (const TMVB& A, const Teuchos::Range1D& index, TMVB& mv)
380  {
381  const int numColsA = A.domain()->dim();
382  const int numColsMv = mv.domain()->dim();
383  // 'index' indexes into mv; it's the index set of the target.
384  const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
385  // We can't take more columns out of A than A has.
386  const bool validSource = index.size() <= numColsA;
387 
388  if (! validIndex || ! validSource)
389  {
390  std::ostringstream os;
391  os << "Anasazi::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar> "
392  ">::SetBlock(A, [" << index.lbound() << ", " << index.ubound()
393  << "], mv): ";
394  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
395  os.str() << "Range lower bound must be nonnegative.");
396  TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
397  os.str() << "Range upper bound must be less than "
398  "the number of columns " << numColsA << " in the "
399  "'mv' output argument.");
400  TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
401  os.str() << "Range must have no more elements than"
402  " the number of columns " << numColsA << " in the "
403  "'A' input argument.");
404  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
405  }
406 
407  // View of the relevant column(s) of the target multivector mv.
408  // We avoid view creation overhead by only creating a view if
409  // the index range is different than [0, (# columns in mv) - 1].
410  Teuchos::RCP<TMVB> mv_view;
411  if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
412  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
413  else
414  mv_view = mv.subView (index);
415 
416  // View of the relevant column(s) of the source multivector A.
417  // If A has fewer columns than mv_view, then create a view of
418  // the first index.size() columns of A.
419  Teuchos::RCP<const TMVB> A_view;
420  if (index.size() == numColsA)
421  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
422  else
423  A_view = A.subView (Teuchos::Range1D(0, index.size()-1));
424 
425  // Copy the data to the destination multivector.
426  Thyra::assign (Teuchos::outArg (*mv_view), *A_view);
427  }
428 
429  static void
430  Assign (const TMVB& A, TMVB& mv)
431  {
432  using Teuchos::ptr;
433  using Teuchos::Range1D;
434  using Teuchos::RCP;
435 
436  const int numColsA = A.domain()->dim();
437  const int numColsMv = mv.domain()->dim();
438  if (numColsA > numColsMv) {
439  const std::string prefix ("Anasazi::MultiVecTraits<Scalar, "
440  "Thyra::MultiVectorBase<Scalar>"
441  " >::Assign(A, mv): ");
442  TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
443  prefix << "Input multivector 'A' has "
444  << numColsA << " columns, but output multivector "
445  "'mv' has only " << numColsMv << " columns.");
446  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
447  }
448  // Copy the data to the destination multivector.
449  if (numColsA == numColsMv) {
450  Thyra::assign (outArg (mv), A);
451  }
452  else {
453  RCP<TMVB> mv_view = CloneViewNonConst (mv, Range1D(0, numColsA-1));
454  Thyra::assign (outArg (*mv_view), A);
455  }
456  }
457 
460  static void MvRandom( TMVB & mv )
461  {
462  // Thyra::randomize generates via a uniform distribution on [l,u]
463  // We will use this to generate on [-1,1]
464  Thyra::randomize(-Teuchos::ScalarTraits<ScalarType>::one(),
465  Teuchos::ScalarTraits<ScalarType>::one(),
466  Teuchos::outArg (mv));
467  }
468 
471  static void
472  MvInit (TMVB& mv,
473  ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero())
474  {
475  Thyra::assign (Teuchos::outArg (mv), alpha);
476  }
477 
479 
482 
485  static void MvPrint( const TMVB & mv, std::ostream& os )
486  {
487  Teuchos::RCP<Teuchos::FancyOStream> out = Teuchos::getFancyOStream(Teuchos::rcp(&os,false));
488  out->setf(std::ios_base::scientific);
489  out->precision(16);
490  mv.describe(*out,Teuchos::VERB_EXTREME);
491  }
492 
494 
495  };
496 
497 
499  //
500  // Implementation of the Anasazi::OperatorTraits for Thyra::LinearOpBase
501  //
503 
513  template <class ScalarType>
514  class OperatorTraits < ScalarType, Thyra::MultiVectorBase<ScalarType>, Thyra::LinearOpBase<ScalarType> >
515  {
516  public:
517 
521  static void Apply ( const Thyra::LinearOpBase< ScalarType >& Op, const Thyra::MultiVectorBase< ScalarType > & x, Thyra::MultiVectorBase< ScalarType > & y )
522  {
523  Op.apply(Thyra::NOTRANS,x,Teuchos::outArg (y), Teuchos::ScalarTraits<ScalarType>::one(),Teuchos::ScalarTraits<ScalarType>::zero());
524  }
525 
526  };
527 
528 } // end of Anasazi namespace
529 
530 #endif
531 // end of file ANASAZI_THYRA_ADAPTER_HPP
static void MvNorm(const TMVB &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec)
Compute the 2-norm of each individual vector of mv. Upon return, normvec[i] holds the value of ...
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 vector (deep copy)...
static void MvDot(const TMVB &mv, const TMVB &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
static void MvTimesMatAddMv(const ScalarType alpha, const TMVB &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, TMVB &mv)
Update mv with .
static void Assign(const MV &A, MV &mv)
mv := A
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
static void MvPrint(const TMVB &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
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 .
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Virtual base class which defines basic traits for the operator type.
static void MvAddMv(const ScalarType alpha, const TMVB &A, const ScalarType beta, const TMVB &B, TMVB &mv)
Replace mv with .
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
static void MvRandom(TMVB &mv)
Replace the vectors in mv with random vectors.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static int GetNumberVecs(const TMVB &mv)
Obtain the number of vectors in mv.
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 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 void MvInit(TMVB &mv, ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv)
Creates a new MultiVectorBase and copies contents of mv into the new vector (deep copy)...
static void MvScale(TMVB &mv, const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
static ptrdiff_t GetGlobalLength(const TMVB &mv)
Obtain the vector length of mv.
static void Apply(const Thyra::LinearOpBase< ScalarType > &Op, const Thyra::MultiVectorBase< ScalarType > &x, Thyra::MultiVectorBase< ScalarType > &y)
This method takes the MultiVectorBase x and applies the LinearOpBase Op to it resulting in the MultiV...
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 void MvScale(TMVB &mv, const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
static Teuchos::RCP< TMVB > Clone(const TMVB &mv, const int numvecs)
Creates a new empty MultiVectorBase containing numvecs columns.