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 ();
165  TEUCHOS_TEST_FOR_EXCEPTION(
166  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
167  std::runtime_error, fnName << ": All indices must be nonnegative.");
168  TEUCHOS_TEST_FOR_EXCEPTION(
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() << "]): ";
205  TEUCHOS_TEST_FOR_EXCEPTION(
206  index.size() == 0, std::invalid_argument,
207  os.str() << "Empty index range is not allowed.");
208  TEUCHOS_TEST_FOR_EXCEPTION(
209  index.lbound() < 0, std::invalid_argument,
210  os.str() << "Index range includes negative index/ices, which is not "
211  "allowed.");
212  TEUCHOS_TEST_FOR_EXCEPTION(
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 ();
230  TEUCHOS_TEST_FOR_EXCEPTION(
231  index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
232  std::invalid_argument,
233  fnName << ": All indices must be nonnegative.");
234  TEUCHOS_TEST_FOR_EXCEPTION(
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() << "]): ";
268  TEUCHOS_TEST_FOR_EXCEPTION(
269  index.size() == 0, std::invalid_argument,
270  os.str() << "Empty index range is not allowed.");
271  TEUCHOS_TEST_FOR_EXCEPTION(
272  index.lbound() < 0, std::invalid_argument,
273  os.str() << "Index range includes negative inde{x,ices}, which is "
274  "not allowed.");
275  TEUCHOS_TEST_FOR_EXCEPTION(
276  index.ubound() >= numCols, std::invalid_argument,
277  os.str() << "Index range exceeds number of vectors " << numCols
278  << " in the input multivector.");
279  TEUCHOS_TEST_FOR_EXCEPTION(
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 
288  static Teuchos::RCP<const MV>
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 ();
295  TEUCHOS_TEST_FOR_EXCEPTION(
296  *std::min_element (index.begin (), index.end ()) < 0,
297  std::invalid_argument,
298  fnName << ": All indices must be nonnegative.");
299  TEUCHOS_TEST_FOR_EXCEPTION(
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 
319  static Teuchos::RCP<const MV>
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.");
337  TEUCHOS_TEST_FOR_EXCEPTION(
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,
360  const Teuchos::SerialDenseMatrix<int, Scalar>& B,
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");
371  Teuchos::RCP<Teuchos::Time> timer =
372  Teuchos::TimeMonitor::lookupCounter (timerName);
373  if (timer.is_null ()) {
374  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
375  }
376  TEUCHOS_TEST_FOR_EXCEPTION(
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,
433  Teuchos::SerialDenseMatrix<int,Scalar>& C)
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  }
448  TEUCHOS_TEST_FOR_EXCEPTION(
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,
495  Teuchos::ScalarTraits<Scalar>::zero ());
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 ();
511  TEUCHOS_TEST_FOR_EXCEPTION(
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
518  TEUCHOS_TEST_FOR_EXCEPTION(
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
537  TEUCHOS_TEST_FOR_EXCEPTION(
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
545  Teuchos::ArrayView<magnitude_type> av (normvec);
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
556  TEUCHOS_TEST_FOR_EXCEPTION(
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): ";
590  TEUCHOS_TEST_FOR_EXCEPTION(
591  maxInt < A.getNumVectors (), std::range_error, os.str () << "Number "
592  "of columns (size_t) in the input MultiVector 'A' overflows int.");
593  TEUCHOS_TEST_FOR_EXCEPTION(
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): ";
610  TEUCHOS_TEST_FOR_EXCEPTION(
611  index.lbound() < 0, std::invalid_argument,
612  os.str() << "Range lower bound must be nonnegative.");
613  TEUCHOS_TEST_FOR_EXCEPTION(
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.");
617  TEUCHOS_TEST_FOR_EXCEPTION(
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.");
621  TEUCHOS_TEST_FOR_EXCEPTION(
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) {
664  TEUCHOS_TEST_FOR_EXCEPTION(
665  maxInt < A.getNumVectors(), std::range_error,
666  errPrefix << "Number of columns in the input multivector 'A' "
667  "(a size_t) overflows int.");
668  TEUCHOS_TEST_FOR_EXCEPTION(
669  maxInt < mv.getNumVectors(), std::range_error,
670  errPrefix << "Number of columns in the output multivector 'mv' "
671  "(a size_t) overflows int.");
672  TEUCHOS_TEST_FOR_EXCEPTION(
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) {
679  TEUCHOS_TEST_FOR_EXCEPTION(
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  mv.print (os);
706  }
707 
708 #ifdef HAVE_ANASAZI_TSQR
709  typedef Tpetra::TsqrAdaptor<Tpetra::MultiVector<Scalar, LO, GO, Node> > tsqr_adaptor_type;
712 #endif // HAVE_ANASAZI_TSQR
713  };
714 
715 
717  template <class Scalar, class LO, class GO, class Node>
718  class OperatorTraits<Scalar,
719  Tpetra::MultiVector<Scalar,LO,GO,Node>,
720  Tpetra::Operator<Scalar,LO,GO,Node> >
721  {
722  public:
723  static void
724  Apply (const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
725  const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
726  Tpetra::MultiVector<Scalar,LO,GO,Node>& Y)
727  {
728  Op.apply (X, Y, Teuchos::NO_TRANS);
729  }
730  };
731 
732 
733 template<class ST, class LO, class GO, class NT>
734 struct OutputStreamTraits<Tpetra::Operator<ST, LO, GO, NT> > {
735  typedef Tpetra::Operator<ST, LO, GO, NT> operator_type;
736 
737  static Teuchos::RCP<Teuchos::FancyOStream>
738  getOutputStream (const operator_type& op, int rootRank = 0)
739  {
740  Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
741  Teuchos::RCP<const Teuchos::Comm<int> > comm = (op.getDomainMap())->getComm ();
742 
743  // Select minimum MPI rank as the root rank for printing.
744  const int myRank = comm.is_null () ? 0 : comm->getRank ();
745  const int numProcs = comm.is_null () ? 1 : comm->getSize ();
746  if (rootRank < 0)
747  {
748  Teuchos::reduceAll(*comm,Teuchos::REDUCE_SUM,1,&myRank,&rootRank);
749  }
750 
751  // This is irreversible, but that's only a problem if the input std::ostream
752  // is actually a Teuchos::FancyOStream on which this method has been
753  // called before, with a different root rank.
754  fos->setProcRankAndSize (myRank, numProcs);
755  fos->setOutputToRootOnly (rootRank);
756  return fos;
757  }
758 };
759 
760 
761 } // end of Anasazi namespace
762 
763 #endif
764 // 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.
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.
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...
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...
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.
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.