Anasazi  Version of the Day
AnasaziTpetraAdapter.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 
29 #ifndef ANASAZI_TPETRA_ADAPTER_HPP
30 #define ANASAZI_TPETRA_ADAPTER_HPP
31 
68 
69 #include <Tpetra_MultiVector.hpp>
70 #include <Tpetra_Operator.hpp>
71 
72 #include <Teuchos_Array.hpp>
73 #include <Teuchos_Assert.hpp>
74 #include <Teuchos_DefaultSerialComm.hpp>
75 #include <Teuchos_ScalarTraits.hpp>
76 
77 #include <AnasaziConfigDefs.hpp>
78 #include <AnasaziTypes.hpp>
81 
82 #ifdef HAVE_ANASAZI_TSQR
83 # include <Tpetra_TsqrAdaptor.hpp>
84 #endif // HAVE_ANASAZI_TSQR
85 
86 
87 namespace Anasazi {
88 
100  template<class Scalar, class LO, class GO, class Node>
101  class MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node> > {
102  typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
103  public:
109  static Teuchos::RCP<MV> Clone (const MV& X, const int numVecs) {
110  Teuchos::RCP<MV> Y (new MV (X.getMap (), numVecs, false));
111  Y->setCopyOrView (Teuchos::View);
112  return Y;
113  }
114 
116  static Teuchos::RCP<MV> CloneCopy (const MV& X)
117  {
118  // Make a deep copy of X. The one-argument copy constructor
119  // does a shallow copy by default; the second argument tells it
120  // to do a deep copy.
121  Teuchos::RCP<MV> X_copy (new MV (X, Teuchos::Copy));
122  // Make Tpetra::MultiVector use the new view semantics. This is
123  // a no-op for the Kokkos refactor version of Tpetra; it only
124  // does something for the "classic" version of Tpetra. This
125  // shouldn't matter because Belos only handles MV through RCP
126  // and through this interface anyway, but it doesn't hurt to set
127  // it and make sure that it works.
128  X_copy->setCopyOrView (Teuchos::View);
129  return X_copy;
130  }
131 
143  static Teuchos::RCP<MV>
144  CloneCopy (const MV& mv, const std::vector<int>& index)
145  {
146 #ifdef HAVE_TPETRA_DEBUG
147  const char fnName[] = "Anasazi::MultiVecTraits::CloneCopy(mv,index)";
148  const size_t inNumVecs = mv.getNumVectors ();
150  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
151  std::runtime_error, fnName << ": All indices must be nonnegative.");
153  index.size () > 0 &&
154  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
155  std::runtime_error,
156  fnName << ": All indices must be strictly less than the number of "
157  "columns " << inNumVecs << " of the input multivector mv.");
158 #endif // HAVE_TPETRA_DEBUG
159 
160  // Tpetra wants an array of size_t, not of int.
161  Teuchos::Array<size_t> columns (index.size ());
162  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
163  columns[j] = index[j];
164  }
165  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
166  // continuous column index range in MultiVector::subCopy, so we
167  // don't have to check here.
168  Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
169  X_copy->setCopyOrView (Teuchos::View);
170  return X_copy;
171  }
172 
179  static Teuchos::RCP<MV>
180  CloneCopy (const MV& mv, const Teuchos::Range1D& index)
181  {
182  const bool validRange = index.size() > 0 &&
183  index.lbound() >= 0 &&
184  index.ubound() < GetNumberVecs(mv);
185  if (! validRange) { // invalid range; generate error message
186  std::ostringstream os;
187  os << "Anasazi::MultiVecTraits::CloneCopy(mv,index=["
188  << index.lbound() << "," << index.ubound() << "]): ";
190  index.size() == 0, std::invalid_argument,
191  os.str() << "Empty index range is not allowed.");
193  index.lbound() < 0, std::invalid_argument,
194  os.str() << "Index range includes negative index/ices, which is not "
195  "allowed.");
197  index.ubound() >= GetNumberVecs(mv), std::invalid_argument,
198  os.str() << "Index range exceeds number of vectors "
199  << mv.getNumVectors() << " in the input multivector.");
200  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
201  os.str() << "Should never get here!");
202  }
203  Teuchos::RCP<MV> X_copy = mv.subCopy (index);
204  X_copy->setCopyOrView (Teuchos::View);
205  return X_copy;
206  }
207 
208  static Teuchos::RCP<MV>
209  CloneViewNonConst (MV& mv, const std::vector<int>& index)
210  {
211 #ifdef HAVE_TPETRA_DEBUG
212  const char fnName[] = "Anasazi::MultiVecTraits::CloneViewNonConst(mv,index)";
213  const size_t numVecs = mv.getNumVectors ();
215  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
216  std::invalid_argument,
217  fnName << ": All indices must be nonnegative.");
219  index.size () > 0 &&
220  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
221  std::invalid_argument,
222  fnName << ": All indices must be strictly less than the number of "
223  "columns " << numVecs << " in the input MultiVector mv.");
224 #endif // HAVE_TPETRA_DEBUG
225 
226  // Tpetra wants an array of size_t, not of int.
227  Teuchos::Array<size_t> columns (index.size ());
228  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
229  columns[j] = index[j];
230  }
231  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
232  // continuous column index range in
233  // MultiVector::subViewNonConst, so we don't have to check here.
234  Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
235  X_view->setCopyOrView (Teuchos::View);
236  return X_view;
237  }
238 
239  static Teuchos::RCP<MV>
240  CloneViewNonConst (MV& mv, const Teuchos::Range1D& index)
241  {
242  // NOTE (mfh 11 Jan 2011) We really should check for possible
243  // overflow of int here. However, the number of columns in a
244  // multivector typically fits in an int.
245  const int numCols = static_cast<int> (mv.getNumVectors());
246  const bool validRange = index.size() > 0 &&
247  index.lbound() >= 0 && index.ubound() < numCols;
248  if (! validRange) {
249  std::ostringstream os;
250  os << "Belos::MultiVecTraits::CloneViewNonConst(mv,index=["
251  << index.lbound() << ", " << index.ubound() << "]): ";
253  index.size() == 0, std::invalid_argument,
254  os.str() << "Empty index range is not allowed.");
256  index.lbound() < 0, std::invalid_argument,
257  os.str() << "Index range includes negative inde{x,ices}, which is "
258  "not allowed.");
260  index.ubound() >= numCols, std::invalid_argument,
261  os.str() << "Index range exceeds number of vectors " << numCols
262  << " in the input multivector.");
264  true, std::logic_error,
265  os.str() << "Should never get here!");
266  }
267  Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
268  X_view->setCopyOrView (Teuchos::View);
269  return X_view;
270  }
271 
273  CloneView (const MV& mv, const std::vector<int>& index)
274  {
275 #ifdef HAVE_TPETRA_DEBUG
276  const char fnName[] = "Belos::MultiVecTraits<Scalar, "
277  "Tpetra::MultiVector<...> >::CloneView(mv,index)";
278  const size_t numVecs = mv.getNumVectors ();
280  *std::min_element (index.begin (), index.end ()) < 0,
281  std::invalid_argument,
282  fnName << ": All indices must be nonnegative.");
284  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
285  std::invalid_argument,
286  fnName << ": All indices must be strictly less than the number of "
287  "columns " << numVecs << " in the input MultiVector mv.");
288 #endif // HAVE_TPETRA_DEBUG
289 
290  // Tpetra wants an array of size_t, not of int.
291  Teuchos::Array<size_t> columns (index.size ());
292  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
293  columns[j] = index[j];
294  }
295  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
296  // continuous column index range in MultiVector::subView, so we
297  // don't have to check here.
298  Teuchos::RCP<const MV> X_view = mv.subView (columns);
299  Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
300  return X_view;
301  }
302 
304  CloneView (const MV& mv, const Teuchos::Range1D& index)
305  {
306  // NOTE (mfh 11 Jan 2011) We really should check for possible
307  // overflow of int here. However, the number of columns in a
308  // multivector typically fits in an int.
309  const int numCols = static_cast<int> (mv.getNumVectors());
310  const bool validRange = index.size() > 0 &&
311  index.lbound() >= 0 && index.ubound() < numCols;
312  if (! validRange) {
313  std::ostringstream os;
314  os << "Anasazi::MultiVecTraits::CloneView(mv, index=["
315  << index.lbound () << ", " << index.ubound() << "]): ";
316  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
317  os.str() << "Empty index range is not allowed.");
318  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
319  os.str() << "Index range includes negative index/ices, which is not "
320  "allowed.");
322  index.ubound() >= numCols, std::invalid_argument,
323  os.str() << "Index range exceeds number of vectors " << numCols
324  << " in the input multivector.");
325  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
326  os.str() << "Should never get here!");
327  }
328  Teuchos::RCP<const MV> X_view = mv.subView (index);
329  Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
330  return X_view;
331  }
332 
333  static ptrdiff_t GetGlobalLength( const MV& mv ) {
334  return static_cast<ptrdiff_t> (mv.getGlobalLength ());
335  }
336 
337  static int GetNumberVecs (const MV& mv) {
338  return static_cast<int> (mv.getNumVectors ());
339  }
340 
341  static void
342  MvTimesMatAddMv (Scalar alpha,
343  const MV& A,
345  Scalar beta,
346  MV& mv)
347  {
348  using Teuchos::ArrayView;
349  using Teuchos::Comm;
350  using Teuchos::rcpFromRef;
351  typedef Tpetra::Map<LO, GO, Node> map_type;
352 
353 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
354  const std::string timerName ("Anasazi::MVT::MvTimesMatAddMv");
357  if (timer.is_null ()) {
358  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
359  }
361  timer.is_null (), std::logic_error,
362  "Anasazi::MultiVecTraits::MvTimesMatAddMv: "
363  "Failed to look up timer \"" << timerName << "\". "
364  "Please report this bug to the Belos developers.");
365 
366  // This starts the timer. It will be stopped on scope exit.
367  Teuchos::TimeMonitor timeMon (*timer);
368 #endif
369 
370  // Check if B is 1-by-1, in which case we can just call update()
371  if (B.numRows () == 1 && B.numCols () == 1) {
372  mv.update (alpha*B(0,0), A, beta);
373  return;
374  }
375 
376  // Create local map
377  Teuchos::SerialComm<int> serialComm;
378  map_type LocalMap (B.numRows (), A.getMap ()->getIndexBase (),
379  rcpFromRef<const Comm<int> > (serialComm),
380  Tpetra::LocallyReplicated, A.getMap ()->getNode ());
381  // encapsulate Teuchos::SerialDenseMatrix data in ArrayView
382  ArrayView<const Scalar> Bvalues (B.values (), B.stride () * B.numCols ());
383  // create locally replicated MultiVector with a copy of this data
384  MV B_mv (rcpFromRef (LocalMap), Bvalues, B.stride (), B.numCols ());
385  mv.multiply (Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
386  }
387 
395  static void
396  MvAddMv (Scalar alpha,
397  const MV& A,
398  Scalar beta,
399  const MV& B,
400  MV& mv)
401  {
402  mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
403  }
404 
405  static void MvScale (MV& mv, Scalar alpha) {
406  mv.scale (alpha);
407  }
408 
409  static void MvScale (MV& mv, const std::vector<Scalar>& alphas) {
410  mv.scale (alphas);
411  }
412 
413  static void
414  MvTransMv (const Scalar alpha,
415  const MV& A,
416  const MV& B,
418  {
419  using Tpetra::LocallyReplicated;
420  using Teuchos::Comm;
421  using Teuchos::RCP;
422  using Teuchos::rcp;
423  using Teuchos::TimeMonitor;
424  typedef Tpetra::Map<LO,GO,Node> map_type;
425 
426 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
427  const std::string timerName ("Anasazi::MVT::MvTransMv");
428  RCP<Teuchos::Time> timer = TimeMonitor::lookupCounter (timerName);
429  if (timer.is_null ()) {
430  timer = TimeMonitor::getNewCounter (timerName);
431  }
433  timer.is_null (), std::logic_error, "Anasazi::MvTransMv: "
434  "Failed to look up timer \"" << timerName << "\". "
435  "Please report this bug to the Belos developers.");
436 
437  // This starts the timer. It will be stopped on scope exit.
438  TimeMonitor timeMon (*timer);
439 #endif // HAVE_ANASAZI_TPETRA_TIMERS
440 
441  // Form alpha * A^H * B, then copy into the SerialDenseMatrix.
442  // We will create a multivector C_mv from a a local map. This
443  // map has a serial comm, the purpose being to short-circuit the
444  // MultiVector::reduce() call at the end of
445  // MultiVector::multiply(). Otherwise, the reduced multivector
446  // data would be copied back to the GPU, only to turn around and
447  // have to get it back here. This saves us a round trip for
448  // this data.
449  const int numRowsC = C.numRows ();
450  const int numColsC = C.numCols ();
451  const int strideC = C.stride ();
452 
453  // Check if numRowsC == numColsC == 1, in which case we can call dot()
454  if (numRowsC == 1 && numColsC == 1) {
455  if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
456  // Short-circuit, as required by BLAS semantics.
457  C(0,0) = alpha;
458  return;
459  }
460  A.dot (B, Teuchos::ArrayView<Scalar> (C.values (), 1));
461  if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
462  C(0,0) *= alpha;
463  }
464  return;
465  }
466 
467  // get comm
468  RCP<const Comm<int> > pcomm = A.getMap ()->getComm ();
469 
470  // create local map with comm
471  RCP<const map_type> LocalMap =
472  rcp (new map_type (numRowsC, 0, pcomm, LocallyReplicated,
473  A.getMap ()->getNode ()));
474  // create local multivector to hold the result
475  const bool INIT_TO_ZERO = true;
476  MV C_mv (LocalMap, numColsC, INIT_TO_ZERO);
477 
478  // multiply result into local multivector
479  C_mv.multiply (Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, alpha, A, B,
481 
482  // create arrayview encapsulating the Teuchos::SerialDenseMatrix
483  Teuchos::ArrayView<Scalar> C_view (C.values (), strideC * numColsC);
484 
485  // No accumulation to do; simply extract the multivector data
486  // into C. Extract a copy of the result into the array view
487  // (and therefore, the SerialDenseMatrix).
488  C_mv.get1dCopy (C_view, strideC);
489  }
490 
492  static void
493  MvDot (const MV& A, const MV& B, std::vector<Scalar> &dots)
494  {
495  const size_t numVecs = A.getNumVectors ();
497  numVecs != B.getNumVectors (), std::invalid_argument,
498  "Anasazi::MultiVecTraits::MvDot(A,B,dots): "
499  "A and B must have the same number of columns. "
500  "A has " << numVecs << " column(s), "
501  "but B has " << B.getNumVectors () << " column(s).");
502 #ifdef HAVE_TPETRA_DEBUG
504  dots.size() < numVecs, std::invalid_argument,
505  "Anasazi::MultiVecTraits::MvDot(A,B,dots): "
506  "The output array 'dots' must have room for all dot products. "
507  "A and B each have " << numVecs << " column(s), "
508  "but 'dots' only has " << dots.size() << " entry(/ies).");
509 #endif // HAVE_TPETRA_DEBUG
510 
511  Teuchos::ArrayView<Scalar> av (dots);
512  A.dot (B, av (0, numVecs));
513  }
514 
516  static void
517  MvNorm (const MV& mv,
518  std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
519  {
520  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
521 #ifdef HAVE_TPETRA_DEBUG
523  normvec.size () < static_cast<std::vector<int>::size_type> (mv.getNumVectors ()),
524  std::invalid_argument,
525  "Anasazi::MultiVecTraits::MvNorm(mv,normvec): The normvec output "
526  "argument must have at least as many entries as the number of vectors "
527  "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
528  << " < mv.getNumVectors() = " << mv.getNumVectors () << ".");
529 #endif // HAVE_TPETRA_DEBUG
531  mv.norm2 (av (0, mv.getNumVectors ()));
532  }
533 
534  static void
535  SetBlock (const MV& A, const std::vector<int>& index, MV& mv)
536  {
537  using Teuchos::Range1D;
538  using Teuchos::RCP;
539  const size_t inNumVecs = A.getNumVectors ();
540 #ifdef HAVE_TPETRA_DEBUG
542  inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
543  "Anasazi::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must "
544  "have no more entries as the number of columns in the input MultiVector"
545  " A. A.getNumVectors() = " << inNumVecs << " < index.size () = "
546  << index.size () << ".");
547 #endif // HAVE_TPETRA_DEBUG
548  RCP<MV> mvsub = CloneViewNonConst (mv, index);
549  if (inNumVecs > static_cast<size_t> (index.size ())) {
550  RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
551  Tpetra::deep_copy (*mvsub, *Asub);
552  } else {
553  Tpetra::deep_copy (*mvsub, A);
554  }
555  }
556 
557  static void
558  SetBlock (const MV& A, const Teuchos::Range1D& index, MV& mv)
559  {
560  // Range1D bounds are signed; size_t is unsigned.
561  // Assignment of Tpetra::MultiVector is a deep copy.
562 
563  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
564  // fair to assume that the number of vectors won't overflow int,
565  // since the typical use case of multivectors involves few
566  // columns, but it's friendly to check just in case.
567  const size_t maxInt =
568  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
569  const bool overflow =
570  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
571  if (overflow) {
572  std::ostringstream os;
573  os << "Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
574  << ", " << index.ubound () << "], mv): ";
576  maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
577  "of columns (size_t) in the input MultiVector 'A' overflows int.");
579  maxInt < mv.getNumVectors (), std::range_error, os.str () << "Number "
580  "of columns (size_t) in the output MultiVector 'mv' overflows int.");
581  }
582  // We've already validated the static casts above.
583  const int numColsA = static_cast<int> (A.getNumVectors ());
584  const int numColsMv = static_cast<int> (mv.getNumVectors ());
585  // 'index' indexes into mv; it's the index set of the target.
586  const bool validIndex =
587  index.lbound () >= 0 && index.ubound () < numColsMv;
588  // We can't take more columns out of A than A has.
589  const bool validSource = index.size () <= numColsA;
590 
591  if (! validIndex || ! validSource) {
592  std::ostringstream os;
593  os << "Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
594  << ", " << index.ubound () << "], mv): ";
596  index.lbound() < 0, std::invalid_argument,
597  os.str() << "Range lower bound must be nonnegative.");
599  index.ubound() >= numColsMv, std::invalid_argument,
600  os.str() << "Range upper bound must be less than the number of "
601  "columns " << numColsA << " in the 'mv' output argument.");
603  index.size() > numColsA, std::invalid_argument,
604  os.str() << "Range must have no more elements than the number of "
605  "columns " << numColsA << " in the 'A' input argument.");
607  true, std::logic_error, "Should never get here!");
608  }
609 
610  // View of the relevant column(s) of the target multivector mv.
611  // We avoid view creation overhead by only creating a view if
612  // the index range is different than [0, (# columns in mv) - 1].
613  Teuchos::RCP<MV> mv_view;
614  if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
615  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
616  } else {
617  mv_view = CloneViewNonConst (mv, index);
618  }
619 
620  // View of the relevant column(s) of the source multivector A.
621  // If A has fewer columns than mv_view, then create a view of
622  // the first index.size() columns of A.
623  Teuchos::RCP<const MV> A_view;
624  if (index.size () == numColsA) {
625  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
626  } else {
627  A_view = CloneView (A, Teuchos::Range1D (0, index.size () - 1));
628  }
629 
630  Tpetra::deep_copy (*mv_view, *A_view);
631  }
632 
633  static void Assign (const MV& A, MV& mv)
634  {
635  const char errPrefix[] = "Anasazi::MultiVecTraits::Assign(A, mv): ";
636 
637  // Range1D bounds are signed; size_t is unsigned.
638  // Assignment of Tpetra::MultiVector is a deep copy.
639 
640  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
641  // fair to assume that the number of vectors won't overflow int,
642  // since the typical use case of multivectors involves few
643  // columns, but it's friendly to check just in case.
644  const size_t maxInt =
645  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
646  const bool overflow =
647  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
648  if (overflow) {
650  maxInt < A.getNumVectors(), std::range_error,
651  errPrefix << "Number of columns in the input multivector 'A' "
652  "(a size_t) overflows int.");
654  maxInt < mv.getNumVectors(), std::range_error,
655  errPrefix << "Number of columns in the output multivector 'mv' "
656  "(a size_t) overflows int.");
658  true, std::logic_error, "Should never get here!");
659  }
660  // We've already validated the static casts above.
661  const int numColsA = static_cast<int> (A.getNumVectors ());
662  const int numColsMv = static_cast<int> (mv.getNumVectors ());
663  if (numColsA > numColsMv) {
665  numColsA > numColsMv, std::invalid_argument,
666  errPrefix << "Input multivector 'A' has " << numColsA << " columns, "
667  "but output multivector 'mv' has only " << numColsMv << " columns.");
668  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
669  }
670  if (numColsA == numColsMv) {
671  Tpetra::deep_copy (mv, A);
672  } else {
673  Teuchos::RCP<MV> mv_view =
674  CloneViewNonConst (mv, Teuchos::Range1D (0, numColsA - 1));
675  Tpetra::deep_copy (*mv_view, A);
676  }
677  }
678 
679  static void MvRandom (MV& mv) {
680  mv.randomize ();
681  }
682 
683  static void
684  MvInit (MV& mv, const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
685  {
686  mv.putScalar (alpha);
687  }
688 
689  static void MvPrint (const MV& mv, std::ostream& os) {
690  mv.print (os);
691  }
692 
693 #ifdef HAVE_ANASAZI_TSQR
694  typedef Tpetra::TsqrAdaptor<Tpetra::MultiVector<Scalar, LO, GO, Node> > tsqr_adaptor_type;
697 #endif // HAVE_ANASAZI_TSQR
698  };
699 
700 
702  template <class Scalar, class LO, class GO, class Node>
703  class OperatorTraits<Scalar,
704  Tpetra::MultiVector<Scalar,LO,GO,Node>,
705  Tpetra::Operator<Scalar,LO,GO,Node> >
706  {
707  public:
708  static void
709  Apply (const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
710  const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
711  Tpetra::MultiVector<Scalar,LO,GO,Node>& Y)
712  {
713  Op.apply (X, Y, Teuchos::NO_TRANS);
714  }
715  };
716 
717 } // end of Anasazi namespace
718 
719 #endif
720 // end of file ANASAZI_TPETRA_ADAPTER_HPP
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvDot(const MV &A, const MV &B, std::vector< Scalar > &dots)
For all columns j of A, set dots[j] := A[j]^T * B[j].
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
static Teuchos::RCP< MV > Clone(const MV &X, const int numVecs)
Create a new MultiVector with numVecs columns.
ScalarType * values() const
static void MvAddMv(Scalar alpha, const MV &A, Scalar beta, const MV &B, MV &mv)
mv := alpha*A + beta*B
Declaration of basic traits for the multivector type.
static RCP< Time > getNewCounter(const std::string &name)
static RCP< Time > lookupCounter(const std::string &name)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
static void Assign(const MV &A, MV &mv)
mv := A
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const std::vector< int > &index)
Create and return a deep copy of the given columns of mv.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
bool is_null() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
OrdinalType numRows() const
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > &normvec)
For all columns j of mv, set normvec[j] = norm(mv[j]).
static Teuchos::RCP< MV > CloneCopy(const MV &X)
Create and return a deep copy of X.
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &B, Teuchos::SerialDenseMatrix< int, ScalarType > &C)
Compute C := alpha * A^H B.
Virtual base class which defines basic traits for the operator type.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Ordinal ubound() const
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
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).
Ordinal lbound() const
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 MV &mv)
Obtain the number of vectors in mv.
Ordinal size() const
Types and exceptions used within Anasazi solvers and interfaces.
static void MvPrint(const MV &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
OrdinalType stride() const
OrdinalType numCols() const
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const Teuchos::Range1D &index)
Create and return a deep copy of the given columns of mv.