Anasazi  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
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 Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef ANASAZI_TPETRA_ADAPTER_HPP
43 #define ANASAZI_TPETRA_ADAPTER_HPP
44 
81 
82 #include <Tpetra_MultiVector.hpp>
83 #include <Tpetra_Operator.hpp>
84 
85 #include <Teuchos_Array.hpp>
86 #include <Teuchos_Assert.hpp>
87 #include <Teuchos_DefaultSerialComm.hpp>
88 #include <Teuchos_CommHelpers.hpp>
89 #include <Teuchos_ScalarTraits.hpp>
90 #include <Teuchos_FancyOStream.hpp>
91 
92 #include <AnasaziConfigDefs.hpp>
93 #include <AnasaziTypes.hpp>
97 
98 #ifdef HAVE_ANASAZI_TSQR
99 # include <Tpetra_TsqrAdaptor.hpp>
100 #endif // HAVE_ANASAZI_TSQR
101 
102 
103 namespace Anasazi {
104 
116  template<class Scalar, class LO, class GO, class Node>
117  class MultiVecTraits<Scalar, Tpetra::MultiVector<Scalar,LO,GO,Node> > {
118  typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
119  public:
125  static Teuchos::RCP<MV> Clone (const MV& X, const int numVecs) {
126  Teuchos::RCP<MV> Y (new MV (X.getMap (), numVecs, false));
127  Y->setCopyOrView (Teuchos::View);
128  return Y;
129  }
130 
132  static Teuchos::RCP<MV> CloneCopy (const MV& X)
133  {
134  // Make a deep copy of X. The one-argument copy constructor
135  // does a shallow copy by default; the second argument tells it
136  // to do a deep copy.
137  Teuchos::RCP<MV> X_copy (new MV (X, Teuchos::Copy));
138  // Make Tpetra::MultiVector use the new view semantics. This is
139  // a no-op for the Kokkos refactor version of Tpetra; it only
140  // does something for the "classic" version of Tpetra. This
141  // shouldn't matter because Belos only handles MV through RCP
142  // and through this interface anyway, but it doesn't hurt to set
143  // it and make sure that it works.
144  X_copy->setCopyOrView (Teuchos::View);
145  return X_copy;
146  }
147 
159  static Teuchos::RCP<MV>
160  CloneCopy (const MV& mv, const std::vector<int>& index)
161  {
162 #ifdef HAVE_TPETRA_DEBUG
163  const char fnName[] = "Anasazi::MultiVecTraits::CloneCopy(mv,index)";
164  const size_t inNumVecs = mv.getNumVectors ();
166  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
167  std::runtime_error, fnName << ": All indices must be nonnegative.");
169  index.size () > 0 &&
170  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
171  std::runtime_error,
172  fnName << ": All indices must be strictly less than the number of "
173  "columns " << inNumVecs << " of the input multivector mv.");
174 #endif // HAVE_TPETRA_DEBUG
175 
176  // Tpetra wants an array of size_t, not of int.
177  Teuchos::Array<size_t> columns (index.size ());
178  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
179  columns[j] = index[j];
180  }
181  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
182  // continuous column index range in MultiVector::subCopy, so we
183  // don't have to check here.
184  Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
185  X_copy->setCopyOrView (Teuchos::View);
186  return X_copy;
187  }
188 
195  static Teuchos::RCP<MV>
196  CloneCopy (const MV& mv, const Teuchos::Range1D& index)
197  {
198  const bool validRange = index.size() > 0 &&
199  index.lbound() >= 0 &&
200  index.ubound() < GetNumberVecs(mv);
201  if (! validRange) { // invalid range; generate error message
202  std::ostringstream os;
203  os << "Anasazi::MultiVecTraits::CloneCopy(mv,index=["
204  << index.lbound() << "," << index.ubound() << "]): ";
206  index.size() == 0, std::invalid_argument,
207  os.str() << "Empty index range is not allowed.");
209  index.lbound() < 0, std::invalid_argument,
210  os.str() << "Index range includes negative index/ices, which is not "
211  "allowed.");
213  index.ubound() >= GetNumberVecs(mv), std::invalid_argument,
214  os.str() << "Index range exceeds number of vectors "
215  << mv.getNumVectors() << " in the input multivector.");
216  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
217  os.str() << "Should never get here!");
218  }
219  Teuchos::RCP<MV> X_copy = mv.subCopy (index);
220  X_copy->setCopyOrView (Teuchos::View);
221  return X_copy;
222  }
223 
224  static Teuchos::RCP<MV>
225  CloneViewNonConst (MV& mv, const std::vector<int>& index)
226  {
227 #ifdef HAVE_TPETRA_DEBUG
228  const char fnName[] = "Anasazi::MultiVecTraits::CloneViewNonConst(mv,index)";
229  const size_t numVecs = mv.getNumVectors ();
231  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
232  std::invalid_argument,
233  fnName << ": All indices must be nonnegative.");
235  index.size () > 0 &&
236  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
237  std::invalid_argument,
238  fnName << ": All indices must be strictly less than the number of "
239  "columns " << numVecs << " in the input MultiVector mv.");
240 #endif // HAVE_TPETRA_DEBUG
241 
242  // Tpetra wants an array of size_t, not of int.
243  Teuchos::Array<size_t> columns (index.size ());
244  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
245  columns[j] = index[j];
246  }
247  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
248  // continuous column index range in
249  // MultiVector::subViewNonConst, so we don't have to check here.
250  Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
251  X_view->setCopyOrView (Teuchos::View);
252  return X_view;
253  }
254 
255  static Teuchos::RCP<MV>
256  CloneViewNonConst (MV& mv, const Teuchos::Range1D& index)
257  {
258  // NOTE (mfh 11 Jan 2011) We really should check for possible
259  // overflow of int here. However, the number of columns in a
260  // multivector typically fits in an int.
261  const int numCols = static_cast<int> (mv.getNumVectors());
262  const bool validRange = index.size() > 0 &&
263  index.lbound() >= 0 && index.ubound() < numCols;
264  if (! validRange) {
265  std::ostringstream os;
266  os << "Belos::MultiVecTraits::CloneViewNonConst(mv,index=["
267  << index.lbound() << ", " << index.ubound() << "]): ";
269  index.size() == 0, std::invalid_argument,
270  os.str() << "Empty index range is not allowed.");
272  index.lbound() < 0, std::invalid_argument,
273  os.str() << "Index range includes negative inde{x,ices}, which is "
274  "not allowed.");
276  index.ubound() >= numCols, std::invalid_argument,
277  os.str() << "Index range exceeds number of vectors " << numCols
278  << " in the input multivector.");
280  true, std::logic_error,
281  os.str() << "Should never get here!");
282  }
283  Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
284  X_view->setCopyOrView (Teuchos::View);
285  return X_view;
286  }
287 
289  CloneView (const MV& mv, const std::vector<int>& index)
290  {
291 #ifdef HAVE_TPETRA_DEBUG
292  const char fnName[] = "Belos::MultiVecTraits<Scalar, "
293  "Tpetra::MultiVector<...> >::CloneView(mv,index)";
294  const size_t numVecs = mv.getNumVectors ();
296  *std::min_element (index.begin (), index.end ()) < 0,
297  std::invalid_argument,
298  fnName << ": All indices must be nonnegative.");
300  static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
301  std::invalid_argument,
302  fnName << ": All indices must be strictly less than the number of "
303  "columns " << numVecs << " in the input MultiVector mv.");
304 #endif // HAVE_TPETRA_DEBUG
305 
306  // Tpetra wants an array of size_t, not of int.
307  Teuchos::Array<size_t> columns (index.size ());
308  for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
309  columns[j] = index[j];
310  }
311  // mfh 14 Aug 2014: Tpetra already detects and optimizes for a
312  // continuous column index range in MultiVector::subView, so we
313  // don't have to check here.
314  Teuchos::RCP<const MV> X_view = mv.subView (columns);
315  Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
316  return X_view;
317  }
318 
320  CloneView (const MV& mv, const Teuchos::Range1D& index)
321  {
322  // NOTE (mfh 11 Jan 2011) We really should check for possible
323  // overflow of int here. However, the number of columns in a
324  // multivector typically fits in an int.
325  const int numCols = static_cast<int> (mv.getNumVectors());
326  const bool validRange = index.size() > 0 &&
327  index.lbound() >= 0 && index.ubound() < numCols;
328  if (! validRange) {
329  std::ostringstream os;
330  os << "Anasazi::MultiVecTraits::CloneView(mv, index=["
331  << index.lbound () << ", " << index.ubound() << "]): ";
332  TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
333  os.str() << "Empty index range is not allowed.");
334  TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
335  os.str() << "Index range includes negative index/ices, which is not "
336  "allowed.");
338  index.ubound() >= numCols, std::invalid_argument,
339  os.str() << "Index range exceeds number of vectors " << numCols
340  << " in the input multivector.");
341  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
342  os.str() << "Should never get here!");
343  }
344  Teuchos::RCP<const MV> X_view = mv.subView (index);
345  Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
346  return X_view;
347  }
348 
349  static ptrdiff_t GetGlobalLength( const MV& mv ) {
350  return static_cast<ptrdiff_t> (mv.getGlobalLength ());
351  }
352 
353  static int GetNumberVecs (const MV& mv) {
354  return static_cast<int> (mv.getNumVectors ());
355  }
356 
357  static void
358  MvTimesMatAddMv (Scalar alpha,
359  const MV& A,
361  Scalar beta,
362  MV& mv)
363  {
364  using Teuchos::ArrayView;
365  using Teuchos::Comm;
366  using Teuchos::rcpFromRef;
367  typedef Tpetra::Map<LO, GO, Node> map_type;
368 
369 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
370  const std::string timerName ("Anasazi::MVT::MvTimesMatAddMv");
373  if (timer.is_null ()) {
374  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
375  }
377  timer.is_null (), std::logic_error,
378  "Anasazi::MultiVecTraits::MvTimesMatAddMv: "
379  "Failed to look up timer \"" << timerName << "\". "
380  "Please report this bug to the Belos developers.");
381 
382  // This starts the timer. It will be stopped on scope exit.
383  Teuchos::TimeMonitor timeMon (*timer);
384 #endif
385 
386  // Check if B is 1-by-1, in which case we can just call update()
387  if (B.numRows () == 1 && B.numCols () == 1) {
388  mv.update (alpha*B(0,0), A, beta);
389  return;
390  }
391 
392  // Create local map
393  Teuchos::SerialComm<int> serialComm;
394  map_type LocalMap (B.numRows (), A.getMap ()->getIndexBase (),
395  rcpFromRef<const Comm<int> > (serialComm),
396  Tpetra::LocallyReplicated);
397  // encapsulate Teuchos::SerialDenseMatrix data in ArrayView
398  ArrayView<const Scalar> Bvalues (B.values (), B.stride () * B.numCols ());
399  // create locally replicated MultiVector with a copy of this data
400  MV B_mv (rcpFromRef (LocalMap), Bvalues, B.stride (), B.numCols ());
401  mv.multiply (Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
402  }
403 
411  static void
412  MvAddMv (Scalar alpha,
413  const MV& A,
414  Scalar beta,
415  const MV& B,
416  MV& mv)
417  {
418  mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
419  }
420 
421  static void MvScale (MV& mv, Scalar alpha) {
422  mv.scale (alpha);
423  }
424 
425  static void MvScale (MV& mv, const std::vector<Scalar>& alphas) {
426  mv.scale (alphas);
427  }
428 
429  static void
430  MvTransMv (const Scalar alpha,
431  const MV& A,
432  const MV& B,
434  {
435  using Tpetra::LocallyReplicated;
436  using Teuchos::Comm;
437  using Teuchos::RCP;
438  using Teuchos::rcp;
439  using Teuchos::TimeMonitor;
440  typedef Tpetra::Map<LO,GO,Node> map_type;
441 
442 #ifdef HAVE_ANASAZI_TPETRA_TIMERS
443  const std::string timerName ("Anasazi::MVT::MvTransMv");
444  RCP<Teuchos::Time> timer = TimeMonitor::lookupCounter (timerName);
445  if (timer.is_null ()) {
446  timer = TimeMonitor::getNewCounter (timerName);
447  }
449  timer.is_null (), std::logic_error, "Anasazi::MvTransMv: "
450  "Failed to look up timer \"" << timerName << "\". "
451  "Please report this bug to the Belos developers.");
452 
453  // This starts the timer. It will be stopped on scope exit.
454  TimeMonitor timeMon (*timer);
455 #endif // HAVE_ANASAZI_TPETRA_TIMERS
456 
457  // Form alpha * A^H * B, then copy into the SerialDenseMatrix.
458  // We will create a multivector C_mv from a a local map. This
459  // map has a serial comm, the purpose being to short-circuit the
460  // MultiVector::reduce() call at the end of
461  // MultiVector::multiply(). Otherwise, the reduced multivector
462  // data would be copied back to the GPU, only to turn around and
463  // have to get it back here. This saves us a round trip for
464  // this data.
465  const int numRowsC = C.numRows ();
466  const int numColsC = C.numCols ();
467  const int strideC = C.stride ();
468 
469  // Check if numRowsC == numColsC == 1, in which case we can call dot()
470  if (numRowsC == 1 && numColsC == 1) {
471  if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
472  // Short-circuit, as required by BLAS semantics.
473  C(0,0) = alpha;
474  return;
475  }
476  A.dot (B, Teuchos::ArrayView<Scalar> (C.values (), 1));
477  if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
478  C(0,0) *= alpha;
479  }
480  return;
481  }
482 
483  // get comm
484  RCP<const Comm<int> > pcomm = A.getMap ()->getComm ();
485 
486  // create local map with comm
487  RCP<const map_type> LocalMap =
488  rcp (new map_type (numRowsC, 0, pcomm, LocallyReplicated));
489  // create local multivector to hold the result
490  const bool INIT_TO_ZERO = true;
491  MV C_mv (LocalMap, numColsC, INIT_TO_ZERO);
492 
493  // multiply result into local multivector
494  C_mv.multiply (Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, alpha, A, B,
496 
497  // create arrayview encapsulating the Teuchos::SerialDenseMatrix
498  Teuchos::ArrayView<Scalar> C_view (C.values (), strideC * numColsC);
499 
500  // No accumulation to do; simply extract the multivector data
501  // into C. Extract a copy of the result into the array view
502  // (and therefore, the SerialDenseMatrix).
503  C_mv.get1dCopy (C_view, strideC);
504  }
505 
507  static void
508  MvDot (const MV& A, const MV& B, std::vector<Scalar> &dots)
509  {
510  const size_t numVecs = A.getNumVectors ();
512  numVecs != B.getNumVectors (), std::invalid_argument,
513  "Anasazi::MultiVecTraits::MvDot(A,B,dots): "
514  "A and B must have the same number of columns. "
515  "A has " << numVecs << " column(s), "
516  "but B has " << B.getNumVectors () << " column(s).");
517 #ifdef HAVE_TPETRA_DEBUG
519  dots.size() < numVecs, std::invalid_argument,
520  "Anasazi::MultiVecTraits::MvDot(A,B,dots): "
521  "The output array 'dots' must have room for all dot products. "
522  "A and B each have " << numVecs << " column(s), "
523  "but 'dots' only has " << dots.size() << " entry(/ies).");
524 #endif // HAVE_TPETRA_DEBUG
525 
526  Teuchos::ArrayView<Scalar> av (dots);
527  A.dot (B, av (0, numVecs));
528  }
529 
531  static void
532  MvNorm (const MV& mv,
533  std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
534  {
535  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
536 #ifdef HAVE_TPETRA_DEBUG
538  normvec.size () < static_cast<std::vector<int>::size_type> (mv.getNumVectors ()),
539  std::invalid_argument,
540  "Anasazi::MultiVecTraits::MvNorm(mv,normvec): The normvec output "
541  "argument must have at least as many entries as the number of vectors "
542  "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
543  << " < mv.getNumVectors() = " << mv.getNumVectors () << ".");
544 #endif // HAVE_TPETRA_DEBUG
546  mv.norm2 (av (0, mv.getNumVectors ()));
547  }
548 
549  static void
550  SetBlock (const MV& A, const std::vector<int>& index, MV& mv)
551  {
552  using Teuchos::Range1D;
553  using Teuchos::RCP;
554  const size_t inNumVecs = A.getNumVectors ();
555 #ifdef HAVE_TPETRA_DEBUG
557  inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
558  "Anasazi::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must "
559  "have no more entries as the number of columns in the input MultiVector"
560  " A. A.getNumVectors() = " << inNumVecs << " < index.size () = "
561  << index.size () << ".");
562 #endif // HAVE_TPETRA_DEBUG
563  RCP<MV> mvsub = CloneViewNonConst (mv, index);
564  if (inNumVecs > static_cast<size_t> (index.size ())) {
565  RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
566  Tpetra::deep_copy (*mvsub, *Asub);
567  } else {
568  Tpetra::deep_copy (*mvsub, A);
569  }
570  }
571 
572  static void
573  SetBlock (const MV& A, const Teuchos::Range1D& index, MV& mv)
574  {
575  // Range1D bounds are signed; size_t is unsigned.
576  // Assignment of Tpetra::MultiVector is a deep copy.
577 
578  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
579  // fair to assume that the number of vectors won't overflow int,
580  // since the typical use case of multivectors involves few
581  // columns, but it's friendly to check just in case.
582  const size_t maxInt =
583  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
584  const bool overflow =
585  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
586  if (overflow) {
587  std::ostringstream os;
588  os << "Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
589  << ", " << index.ubound () << "], mv): ";
591  maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
592  "of columns (size_t) in the input MultiVector 'A' overflows int.");
594  maxInt < mv.getNumVectors (), std::range_error, os.str () << "Number "
595  "of columns (size_t) in the output MultiVector 'mv' overflows int.");
596  }
597  // We've already validated the static casts above.
598  const int numColsA = static_cast<int> (A.getNumVectors ());
599  const int numColsMv = static_cast<int> (mv.getNumVectors ());
600  // 'index' indexes into mv; it's the index set of the target.
601  const bool validIndex =
602  index.lbound () >= 0 && index.ubound () < numColsMv;
603  // We can't take more columns out of A than A has.
604  const bool validSource = index.size () <= numColsA;
605 
606  if (! validIndex || ! validSource) {
607  std::ostringstream os;
608  os << "Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
609  << ", " << index.ubound () << "], mv): ";
611  index.lbound() < 0, std::invalid_argument,
612  os.str() << "Range lower bound must be nonnegative.");
614  index.ubound() >= numColsMv, std::invalid_argument,
615  os.str() << "Range upper bound must be less than the number of "
616  "columns " << numColsA << " in the 'mv' output argument.");
618  index.size() > numColsA, std::invalid_argument,
619  os.str() << "Range must have no more elements than the number of "
620  "columns " << numColsA << " in the 'A' input argument.");
622  true, std::logic_error, "Should never get here!");
623  }
624 
625  // View of the relevant column(s) of the target multivector mv.
626  // We avoid view creation overhead by only creating a view if
627  // the index range is different than [0, (# columns in mv) - 1].
628  Teuchos::RCP<MV> mv_view;
629  if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
630  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
631  } else {
632  mv_view = CloneViewNonConst (mv, index);
633  }
634 
635  // View of the relevant column(s) of the source multivector A.
636  // If A has fewer columns than mv_view, then create a view of
637  // the first index.size() columns of A.
638  Teuchos::RCP<const MV> A_view;
639  if (index.size () == numColsA) {
640  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
641  } else {
642  A_view = CloneView (A, Teuchos::Range1D (0, index.size () - 1));
643  }
644 
645  Tpetra::deep_copy (*mv_view, *A_view);
646  }
647 
648  static void Assign (const MV& A, MV& mv)
649  {
650  const char errPrefix[] = "Anasazi::MultiVecTraits::Assign(A, mv): ";
651 
652  // Range1D bounds are signed; size_t is unsigned.
653  // Assignment of Tpetra::MultiVector is a deep copy.
654 
655  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
656  // fair to assume that the number of vectors won't overflow int,
657  // since the typical use case of multivectors involves few
658  // columns, but it's friendly to check just in case.
659  const size_t maxInt =
660  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
661  const bool overflow =
662  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
663  if (overflow) {
665  maxInt < A.getNumVectors(), std::range_error,
666  errPrefix << "Number of columns in the input multivector 'A' "
667  "(a size_t) overflows int.");
669  maxInt < mv.getNumVectors(), std::range_error,
670  errPrefix << "Number of columns in the output multivector 'mv' "
671  "(a size_t) overflows int.");
673  true, std::logic_error, "Should never get here!");
674  }
675  // We've already validated the static casts above.
676  const int numColsA = static_cast<int> (A.getNumVectors ());
677  const int numColsMv = static_cast<int> (mv.getNumVectors ());
678  if (numColsA > numColsMv) {
680  numColsA > numColsMv, std::invalid_argument,
681  errPrefix << "Input multivector 'A' has " << numColsA << " columns, "
682  "but output multivector 'mv' has only " << numColsMv << " columns.");
683  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
684  }
685  if (numColsA == numColsMv) {
686  Tpetra::deep_copy (mv, A);
687  } else {
688  Teuchos::RCP<MV> mv_view =
689  CloneViewNonConst (mv, Teuchos::Range1D (0, numColsA - 1));
690  Tpetra::deep_copy (*mv_view, A);
691  }
692  }
693 
694  static void MvRandom (MV& mv) {
695  mv.randomize ();
696  }
697 
698  static void
699  MvInit (MV& mv, const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
700  {
701  mv.putScalar (alpha);
702  }
703 
704  static void MvPrint (const MV& mv, std::ostream& os) {
705  Teuchos::FancyOStream fos (Teuchos::rcpFromRef (os));
706  mv.describe (fos, Teuchos::VERB_EXTREME);
707  }
708 
709 #ifdef HAVE_ANASAZI_TSQR
710  typedef Tpetra::TsqrAdaptor<Tpetra::MultiVector<Scalar, LO, GO, Node> > tsqr_adaptor_type;
713 #endif // HAVE_ANASAZI_TSQR
714  };
715 
716 
718  template <class Scalar, class LO, class GO, class Node>
719  class OperatorTraits<Scalar,
720  Tpetra::MultiVector<Scalar,LO,GO,Node>,
721  Tpetra::Operator<Scalar,LO,GO,Node> >
722  {
723  public:
724  static void
725  Apply (const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
726  const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
727  Tpetra::MultiVector<Scalar,LO,GO,Node>& Y)
728  {
729  Op.apply (X, Y, Teuchos::NO_TRANS);
730  }
731  };
732 
733 
734 template<class ST, class LO, class GO, class NT>
735 struct OutputStreamTraits<Tpetra::Operator<ST, LO, GO, NT> > {
736  typedef Tpetra::Operator<ST, LO, GO, NT> operator_type;
737 
739  getOutputStream (const operator_type& op, int rootRank = 0)
740  {
741  Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
742  Teuchos::RCP<const Teuchos::Comm<int> > comm = (op.getDomainMap())->getComm ();
743 
744  // Select minimum MPI rank as the root rank for printing.
745  const int myRank = comm.is_null () ? 0 : comm->getRank ();
746  const int numProcs = comm.is_null () ? 1 : comm->getSize ();
747  if (rootRank < 0)
748  {
749  Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,1,&myRank,&rootRank);
750  }
751 
752  // This is irreversible, but that's only a problem if the input std::ostream
753  // is actually a Teuchos::FancyOStream on which this method has been
754  // called before, with a different root rank.
755  fos->setProcRankAndSize (myRank, numProcs);
756  fos->setOutputToRootOnly (rootRank);
757  return fos;
758  }
759 };
760 
761 
762 } // end of Anasazi namespace
763 
764 #endif
765 // end of file ANASAZI_TPETRA_ADAPTER_HPP
ScalarType * values() const
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
basic_FancyOStream & setProcRankAndSize(const int procRank, const int numProcs)
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.
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
Output managers remove the need for the eigensolver to know any information about the required output...
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.
Ordinal ubound() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
basic_FancyOStream & setOutputToRootOnly(const int rootRank)
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 MV &mv)
Obtain the number of vectors in mv.
Ordinal lbound() const
OrdinalType numCols() const
Ordinal size() const
Types and exceptions used within Anasazi solvers and interfaces.
Abstract class definition for Anasazi output stream.
static void MvPrint(const MV &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
Anasazi&#39;s templated virtual class for constructing an operator that can interface with the OperatorTr...
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const Teuchos::Range1D &index)
Create and return a deep copy of the given columns of mv.
OrdinalType stride() const
OrdinalType numRows() const
bool is_null() const