42 #ifndef ANASAZI_TPETRA_ADAPTER_HPP 43 #define ANASAZI_TPETRA_ADAPTER_HPP 82 #include <Tpetra_MultiVector.hpp> 83 #include <Tpetra_Operator.hpp> 86 #include <Teuchos_Assert.hpp> 87 #include <Teuchos_DefaultSerialComm.hpp> 88 #include <Teuchos_CommHelpers.hpp> 90 #include <Teuchos_FancyOStream.hpp> 98 #ifdef HAVE_ANASAZI_TSQR 99 # include <Tpetra_TsqrAdaptor.hpp> 100 #endif // HAVE_ANASAZI_TSQR 116 template<
class Scalar,
class LO,
class GO,
class Node>
118 typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
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.");
170 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
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 178 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
179 columns[j] = index[j];
198 const bool validRange = index.
size() > 0 &&
202 std::ostringstream os;
203 os <<
"Anasazi::MultiVecTraits::CloneCopy(mv,index=[" 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 " 214 os.str() <<
"Index range exceeds number of vectors " 215 << mv.getNumVectors() <<
" in the input multivector.");
217 os.str() <<
"Should never get here!");
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.");
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 244 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
245 columns[j] = index[j];
261 const int numCols =
static_cast<int> (mv.getNumVectors());
262 const bool validRange = index.
size() > 0 &&
265 std::ostringstream os;
266 os <<
"Belos::MultiVecTraits::CloneViewNonConst(mv,index=[" 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 " 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!");
289 CloneView (
const MV& mv,
const std::vector<int>& index)
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 308 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
309 columns[j] = index[j];
315 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (
Teuchos::View);
325 const int numCols =
static_cast<int> (mv.getNumVectors());
326 const bool validRange = index.
size() > 0 &&
329 std::ostringstream os;
330 os <<
"Anasazi::MultiVecTraits::CloneView(mv, index=[" 333 os.str() <<
"Empty index range is not allowed.");
335 os.str() <<
"Index range includes negative index/ices, which is not " 338 index.
ubound() >= numCols, std::invalid_argument,
339 os.str() <<
"Index range exceeds number of vectors " << numCols
340 <<
" in the input multivector.");
342 os.str() <<
"Should never get here!");
345 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (
Teuchos::View);
350 return static_cast<ptrdiff_t
> (mv.getGlobalLength ());
354 return static_cast<int> (mv.getNumVectors ());
366 using Teuchos::rcpFromRef;
367 typedef Tpetra::Map<LO, GO, Node> map_type;
369 #ifdef HAVE_ANASAZI_TPETRA_TIMERS 370 const std::string timerName (
"Anasazi::MVT::MvTimesMatAddMv");
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.");
388 mv.update (alpha*B(0,0), A, beta);
394 map_type LocalMap (B.
numRows (), A.getMap ()->getIndexBase (),
395 rcpFromRef<const Comm<int> > (serialComm),
396 Tpetra::LocallyReplicated, A.getMap ()->getNode ());
400 MV B_mv (rcpFromRef (LocalMap), Bvalues, B.
stride (), B.
numCols ());
421 static void MvScale (MV& mv, Scalar alpha) {
425 static void MvScale (MV& mv,
const std::vector<Scalar>& alphas) {
435 using Tpetra::LocallyReplicated;
440 typedef Tpetra::Map<LO,GO,Node> map_type;
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);
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.");
454 TimeMonitor timeMon (*timer);
455 #endif // HAVE_ANASAZI_TPETRA_TIMERS 465 const int numRowsC = C.
numRows ();
466 const int numColsC = C.
numCols ();
467 const int strideC = C.
stride ();
470 if (numRowsC == 1 && numColsC == 1) {
484 RCP<const Comm<int> > pcomm = A.getMap ()->getComm ();
487 RCP<const map_type> LocalMap =
488 rcp (
new map_type (numRowsC, 0, pcomm, LocallyReplicated,
489 A.getMap ()->getNode ()));
491 const bool INIT_TO_ZERO =
true;
492 MV C_mv (LocalMap, numColsC, INIT_TO_ZERO);
504 C_mv.get1dCopy (C_view, strideC);
509 MvDot (
const MV& A,
const MV& B, std::vector<Scalar> &dots)
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 528 A.dot (B, av (0, numVecs));
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 ()));
551 SetBlock (
const MV& A,
const std::vector<int>& index, MV& mv)
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 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);
569 Tpetra::deep_copy (*mvsub, A);
583 const size_t maxInt =
585 const bool overflow =
586 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
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.");
599 const int numColsA =
static_cast<int> (A.getNumVectors ());
600 const int numColsMv =
static_cast<int> (mv.getNumVectors ());
602 const bool validIndex =
605 const bool validSource = index.
size () <= numColsA;
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!");
630 if (index.
lbound () == 0 && index.
ubound () + 1 == numColsMv) {
631 mv_view = Teuchos::rcpFromRef (mv);
640 if (index.
size () == numColsA) {
641 A_view = Teuchos::rcpFromRef (A);
646 Tpetra::deep_copy (*mv_view, *A_view);
649 static void Assign (
const MV& A, MV& mv)
651 const char errPrefix[] =
"Anasazi::MultiVecTraits::Assign(A, mv): ";
660 const size_t maxInt =
662 const bool overflow =
663 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
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!");
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.");
686 if (numColsA == numColsMv) {
687 Tpetra::deep_copy (mv, A);
691 Tpetra::deep_copy (*mv_view, A);
702 mv.putScalar (alpha);
705 static void MvPrint (
const MV& mv, std::ostream& os) {
709 #ifdef HAVE_ANASAZI_TSQR 710 typedef Tpetra::TsqrAdaptor<Tpetra::MultiVector<Scalar, LO, GO, Node> > tsqr_adaptor_type;
713 #endif // HAVE_ANASAZI_TSQR 718 template <
class Scalar,
class LO,
class GO,
class Node>
720 Tpetra::MultiVector<Scalar,LO,GO,Node>,
721 Tpetra::Operator<Scalar,LO,GO,Node> >
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)
734 template<
class ST,
class LO,
class GO,
class NT>
736 typedef Tpetra::Operator<ST, LO, GO, NT> operator_type;
739 getOutputStream (
const operator_type& op,
int rootRank = 0)
745 const int myRank = comm.
is_null () ? 0 : comm->getRank ();
746 const int numProcs = comm.
is_null () ? 1 : comm->getSize ();
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.
#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...
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.
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).
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.
OrdinalType stride() const
OrdinalType numCols() const
Anasazi'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.