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 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, A.getMap ()->getNode ());
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  A.getMap ()->getNode ()));
490  // create local multivector to hold the result
491  const bool INIT_TO_ZERO = true;
492  MV C_mv (LocalMap, numColsC, INIT_TO_ZERO);
493 
494  // multiply result into local multivector
495  C_mv.multiply (Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, alpha, A, B,
497 
498  // create arrayview encapsulating the Teuchos::SerialDenseMatrix
499  Teuchos::ArrayView<Scalar> C_view (C.values (), strideC * numColsC);
500 
501  // No accumulation to do; simply extract the multivector data
502  // into C. Extract a copy of the result into the array view
503  // (and therefore, the SerialDenseMatrix).
504  C_mv.get1dCopy (C_view, strideC);
505  }
506 
508  static void
509  MvDot (const MV& A, const MV& B, std::vector<Scalar> &dots)
510  {
511  const size_t numVecs = A.getNumVectors ();
513  numVecs != B.getNumVectors (), std::invalid_argument,
514  "Anasazi::MultiVecTraits::MvDot(A,B,dots): "
515  "A and B must have the same number of columns. "
516  "A has " << numVecs << " column(s), "
517  "but B has " << B.getNumVectors () << " column(s).");
518 #ifdef HAVE_TPETRA_DEBUG
520  dots.size() < numVecs, std::invalid_argument,
521  "Anasazi::MultiVecTraits::MvDot(A,B,dots): "
522  "The output array 'dots' must have room for all dot products. "
523  "A and B each have " << numVecs << " column(s), "
524  "but 'dots' only has " << dots.size() << " entry(/ies).");
525 #endif // HAVE_TPETRA_DEBUG
526 
527  Teuchos::ArrayView<Scalar> av (dots);
528  A.dot (B, av (0, numVecs));
529  }
530 
532  static void
533  MvNorm (const MV& mv,
534  std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
535  {
536  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
537 #ifdef HAVE_TPETRA_DEBUG
539  normvec.size () < static_cast<std::vector<int>::size_type> (mv.getNumVectors ()),
540  std::invalid_argument,
541  "Anasazi::MultiVecTraits::MvNorm(mv,normvec): The normvec output "
542  "argument must have at least as many entries as the number of vectors "
543  "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
544  << " < mv.getNumVectors() = " << mv.getNumVectors () << ".");
545 #endif // HAVE_TPETRA_DEBUG
547  mv.norm2 (av (0, mv.getNumVectors ()));
548  }
549 
550  static void
551  SetBlock (const MV& A, const std::vector<int>& index, MV& mv)
552  {
553  using Teuchos::Range1D;
554  using Teuchos::RCP;
555  const size_t inNumVecs = A.getNumVectors ();
556 #ifdef HAVE_TPETRA_DEBUG
558  inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
559  "Anasazi::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must "
560  "have no more entries as the number of columns in the input MultiVector"
561  " A. A.getNumVectors() = " << inNumVecs << " < index.size () = "
562  << index.size () << ".");
563 #endif // HAVE_TPETRA_DEBUG
564  RCP<MV> mvsub = CloneViewNonConst (mv, index);
565  if (inNumVecs > static_cast<size_t> (index.size ())) {
566  RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
567  Tpetra::deep_copy (*mvsub, *Asub);
568  } else {
569  Tpetra::deep_copy (*mvsub, A);
570  }
571  }
572 
573  static void
574  SetBlock (const MV& A, const Teuchos::Range1D& index, MV& mv)
575  {
576  // Range1D bounds are signed; size_t is unsigned.
577  // Assignment of Tpetra::MultiVector is a deep copy.
578 
579  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
580  // fair to assume that the number of vectors won't overflow int,
581  // since the typical use case of multivectors involves few
582  // columns, but it's friendly to check just in case.
583  const size_t maxInt =
584  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
585  const bool overflow =
586  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
587  if (overflow) {
588  std::ostringstream os;
589  os << "Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
590  << ", " << index.ubound () << "], mv): ";
592  maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
593  "of columns (size_t) in the input MultiVector 'A' overflows int.");
595  maxInt < mv.getNumVectors (), std::range_error, os.str () << "Number "
596  "of columns (size_t) in the output MultiVector 'mv' overflows int.");
597  }
598  // We've already validated the static casts above.
599  const int numColsA = static_cast<int> (A.getNumVectors ());
600  const int numColsMv = static_cast<int> (mv.getNumVectors ());
601  // 'index' indexes into mv; it's the index set of the target.
602  const bool validIndex =
603  index.lbound () >= 0 && index.ubound () < numColsMv;
604  // We can't take more columns out of A than A has.
605  const bool validSource = index.size () <= numColsA;
606 
607  if (! validIndex || ! validSource) {
608  std::ostringstream os;
609  os << "Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
610  << ", " << index.ubound () << "], mv): ";
612  index.lbound() < 0, std::invalid_argument,
613  os.str() << "Range lower bound must be nonnegative.");
615  index.ubound() >= numColsMv, std::invalid_argument,
616  os.str() << "Range upper bound must be less than the number of "
617  "columns " << numColsA << " in the 'mv' output argument.");
619  index.size() > numColsA, std::invalid_argument,
620  os.str() << "Range must have no more elements than the number of "
621  "columns " << numColsA << " in the 'A' input argument.");
623  true, std::logic_error, "Should never get here!");
624  }
625 
626  // View of the relevant column(s) of the target multivector mv.
627  // We avoid view creation overhead by only creating a view if
628  // the index range is different than [0, (# columns in mv) - 1].
629  Teuchos::RCP<MV> mv_view;
630  if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
631  mv_view = Teuchos::rcpFromRef (mv); // Non-const, non-owning RCP
632  } else {
633  mv_view = CloneViewNonConst (mv, index);
634  }
635 
636  // View of the relevant column(s) of the source multivector A.
637  // If A has fewer columns than mv_view, then create a view of
638  // the first index.size() columns of A.
639  Teuchos::RCP<const MV> A_view;
640  if (index.size () == numColsA) {
641  A_view = Teuchos::rcpFromRef (A); // Const, non-owning RCP
642  } else {
643  A_view = CloneView (A, Teuchos::Range1D (0, index.size () - 1));
644  }
645 
646  Tpetra::deep_copy (*mv_view, *A_view);
647  }
648 
649  static void Assign (const MV& A, MV& mv)
650  {
651  const char errPrefix[] = "Anasazi::MultiVecTraits::Assign(A, mv): ";
652 
653  // Range1D bounds are signed; size_t is unsigned.
654  // Assignment of Tpetra::MultiVector is a deep copy.
655 
656  // Tpetra::MultiVector::getNumVectors() returns size_t. It's
657  // fair to assume that the number of vectors won't overflow int,
658  // since the typical use case of multivectors involves few
659  // columns, but it's friendly to check just in case.
660  const size_t maxInt =
661  static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
662  const bool overflow =
663  maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
664  if (overflow) {
666  maxInt < A.getNumVectors(), std::range_error,
667  errPrefix << "Number of columns in the input multivector 'A' "
668  "(a size_t) overflows int.");
670  maxInt < mv.getNumVectors(), std::range_error,
671  errPrefix << "Number of columns in the output multivector 'mv' "
672  "(a size_t) overflows int.");
674  true, std::logic_error, "Should never get here!");
675  }
676  // We've already validated the static casts above.
677  const int numColsA = static_cast<int> (A.getNumVectors ());
678  const int numColsMv = static_cast<int> (mv.getNumVectors ());
679  if (numColsA > numColsMv) {
681  numColsA > numColsMv, std::invalid_argument,
682  errPrefix << "Input multivector 'A' has " << numColsA << " columns, "
683  "but output multivector 'mv' has only " << numColsMv << " columns.");
684  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Should never get here!");
685  }
686  if (numColsA == numColsMv) {
687  Tpetra::deep_copy (mv, A);
688  } else {
689  Teuchos::RCP<MV> mv_view =
690  CloneViewNonConst (mv, Teuchos::Range1D (0, numColsA - 1));
691  Tpetra::deep_copy (*mv_view, A);
692  }
693  }
694 
695  static void MvRandom (MV& mv) {
696  mv.randomize ();
697  }
698 
699  static void
700  MvInit (MV& mv, const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
701  {
702  mv.putScalar (alpha);
703  }
704 
705  static void MvPrint (const MV& mv, std::ostream& os) {
706  mv.print (os);
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
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.
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
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.
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
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
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).
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.
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.
OrdinalType stride() const
OrdinalType numCols() const
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.