29 #ifndef ANASAZI_TPETRA_ADAPTER_HPP 30 #define ANASAZI_TPETRA_ADAPTER_HPP 69 #include <Tpetra_MultiVector.hpp> 70 #include <Tpetra_Operator.hpp> 72 #include <Teuchos_Array.hpp> 73 #include <Teuchos_Assert.hpp> 74 #include <Teuchos_DefaultSerialComm.hpp> 75 #include <Teuchos_ScalarTraits.hpp> 82 #ifdef HAVE_ANASAZI_TSQR 83 # include <Tpetra_TsqrAdaptor.hpp> 84 #endif // HAVE_ANASAZI_TSQR 100 template<
class Scalar,
class LO,
class GO,
class Node>
102 typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
109 static Teuchos::RCP<MV>
Clone (
const MV& X,
const int numVecs) {
110 Teuchos::RCP<MV> Y (
new MV (X.getMap (), numVecs,
false));
111 Y->setCopyOrView (Teuchos::View);
121 Teuchos::RCP<MV> X_copy (
new MV (X, Teuchos::Copy));
128 X_copy->setCopyOrView (Teuchos::View);
143 static Teuchos::RCP<MV>
146 #ifdef HAVE_TPETRA_DEBUG 147 const char fnName[] =
"Anasazi::MultiVecTraits::CloneCopy(mv,index)";
148 const size_t inNumVecs = mv.getNumVectors ();
149 TEUCHOS_TEST_FOR_EXCEPTION(
150 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
151 std::runtime_error, fnName <<
": All indices must be nonnegative.");
152 TEUCHOS_TEST_FOR_EXCEPTION(
154 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= inNumVecs,
156 fnName <<
": All indices must be strictly less than the number of " 157 "columns " << inNumVecs <<
" of the input multivector mv.");
158 #endif // HAVE_TPETRA_DEBUG 161 Teuchos::Array<size_t> columns (index.size ());
162 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
163 columns[j] = index[j];
168 Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
169 X_copy->setCopyOrView (Teuchos::View);
179 static Teuchos::RCP<MV>
182 const bool validRange = index.size() > 0 &&
183 index.lbound() >= 0 &&
186 std::ostringstream os;
187 os <<
"Anasazi::MultiVecTraits::CloneCopy(mv,index=[" 188 << index.lbound() <<
"," << index.ubound() <<
"]): ";
189 TEUCHOS_TEST_FOR_EXCEPTION(
190 index.size() == 0, std::invalid_argument,
191 os.str() <<
"Empty index range is not allowed.");
192 TEUCHOS_TEST_FOR_EXCEPTION(
193 index.lbound() < 0, std::invalid_argument,
194 os.str() <<
"Index range includes negative index/ices, which is not " 196 TEUCHOS_TEST_FOR_EXCEPTION(
198 os.str() <<
"Index range exceeds number of vectors " 199 << mv.getNumVectors() <<
" in the input multivector.");
200 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
201 os.str() <<
"Should never get here!");
203 Teuchos::RCP<MV> X_copy = mv.subCopy (index);
204 X_copy->setCopyOrView (Teuchos::View);
208 static Teuchos::RCP<MV>
211 #ifdef HAVE_TPETRA_DEBUG 212 const char fnName[] =
"Anasazi::MultiVecTraits::CloneViewNonConst(mv,index)";
213 const size_t numVecs = mv.getNumVectors ();
214 TEUCHOS_TEST_FOR_EXCEPTION(
215 index.size () > 0 && *std::min_element (index.begin (), index.end ()) < 0,
216 std::invalid_argument,
217 fnName <<
": All indices must be nonnegative.");
218 TEUCHOS_TEST_FOR_EXCEPTION(
220 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
221 std::invalid_argument,
222 fnName <<
": All indices must be strictly less than the number of " 223 "columns " << numVecs <<
" in the input MultiVector mv.");
224 #endif // HAVE_TPETRA_DEBUG 227 Teuchos::Array<size_t> columns (index.size ());
228 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
229 columns[j] = index[j];
234 Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
235 X_view->setCopyOrView (Teuchos::View);
239 static Teuchos::RCP<MV>
245 const int numCols =
static_cast<int> (mv.getNumVectors());
246 const bool validRange = index.size() > 0 &&
247 index.lbound() >= 0 && index.ubound() < numCols;
249 std::ostringstream os;
250 os <<
"Belos::MultiVecTraits::CloneViewNonConst(mv,index=[" 251 << index.lbound() <<
", " << index.ubound() <<
"]): ";
252 TEUCHOS_TEST_FOR_EXCEPTION(
253 index.size() == 0, std::invalid_argument,
254 os.str() <<
"Empty index range is not allowed.");
255 TEUCHOS_TEST_FOR_EXCEPTION(
256 index.lbound() < 0, std::invalid_argument,
257 os.str() <<
"Index range includes negative inde{x,ices}, which is " 259 TEUCHOS_TEST_FOR_EXCEPTION(
260 index.ubound() >= numCols, std::invalid_argument,
261 os.str() <<
"Index range exceeds number of vectors " << numCols
262 <<
" in the input multivector.");
263 TEUCHOS_TEST_FOR_EXCEPTION(
264 true, std::logic_error,
265 os.str() <<
"Should never get here!");
267 Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
268 X_view->setCopyOrView (Teuchos::View);
272 static Teuchos::RCP<const MV>
273 CloneView (
const MV& mv,
const std::vector<int>& index)
275 #ifdef HAVE_TPETRA_DEBUG 276 const char fnName[] =
"Belos::MultiVecTraits<Scalar, " 277 "Tpetra::MultiVector<...> >::CloneView(mv,index)";
278 const size_t numVecs = mv.getNumVectors ();
279 TEUCHOS_TEST_FOR_EXCEPTION(
280 *std::min_element (index.begin (), index.end ()) < 0,
281 std::invalid_argument,
282 fnName <<
": All indices must be nonnegative.");
283 TEUCHOS_TEST_FOR_EXCEPTION(
284 static_cast<size_t> (*std::max_element (index.begin (), index.end ())) >= numVecs,
285 std::invalid_argument,
286 fnName <<
": All indices must be strictly less than the number of " 287 "columns " << numVecs <<
" in the input MultiVector mv.");
288 #endif // HAVE_TPETRA_DEBUG 291 Teuchos::Array<size_t> columns (index.size ());
292 for (std::vector<int>::size_type j = 0; j < index.size (); ++j) {
293 columns[j] = index[j];
298 Teuchos::RCP<const MV> X_view = mv.subView (columns);
299 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
303 static Teuchos::RCP<const MV>
304 CloneView (
const MV& mv,
const Teuchos::Range1D& index)
309 const int numCols =
static_cast<int> (mv.getNumVectors());
310 const bool validRange = index.size() > 0 &&
311 index.lbound() >= 0 && index.ubound() < numCols;
313 std::ostringstream os;
314 os <<
"Anasazi::MultiVecTraits::CloneView(mv, index=[" 315 << index.lbound () <<
", " << index.ubound() <<
"]): ";
316 TEUCHOS_TEST_FOR_EXCEPTION(index.size() == 0, std::invalid_argument,
317 os.str() <<
"Empty index range is not allowed.");
318 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
319 os.str() <<
"Index range includes negative index/ices, which is not " 321 TEUCHOS_TEST_FOR_EXCEPTION(
322 index.ubound() >= numCols, std::invalid_argument,
323 os.str() <<
"Index range exceeds number of vectors " << numCols
324 <<
" in the input multivector.");
325 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
326 os.str() <<
"Should never get here!");
328 Teuchos::RCP<const MV> X_view = mv.subView (index);
329 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
334 return static_cast<ptrdiff_t
> (mv.getGlobalLength ());
338 return static_cast<int> (mv.getNumVectors ());
344 const Teuchos::SerialDenseMatrix<int, Scalar>& B,
348 using Teuchos::ArrayView;
350 using Teuchos::rcpFromRef;
351 typedef Tpetra::Map<LO, GO, Node> map_type;
353 #ifdef HAVE_ANASAZI_TPETRA_TIMERS 354 const std::string timerName (
"Anasazi::MVT::MvTimesMatAddMv");
355 Teuchos::RCP<Teuchos::Time> timer =
356 Teuchos::TimeMonitor::lookupCounter (timerName);
357 if (timer.is_null ()) {
358 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
360 TEUCHOS_TEST_FOR_EXCEPTION(
361 timer.is_null (), std::logic_error,
362 "Anasazi::MultiVecTraits::MvTimesMatAddMv: " 363 "Failed to look up timer \"" << timerName <<
"\". " 364 "Please report this bug to the Belos developers.");
367 Teuchos::TimeMonitor timeMon (*timer);
371 if (B.numRows () == 1 && B.numCols () == 1) {
372 mv.update (alpha*B(0,0), A, beta);
377 Teuchos::SerialComm<int> serialComm;
378 map_type LocalMap (B.numRows (), A.getMap ()->getIndexBase (),
379 rcpFromRef<const Comm<int> > (serialComm),
380 Tpetra::LocallyReplicated, A.getMap ()->getNode ());
382 ArrayView<const Scalar> Bvalues (B.values (), B.stride () * B.numCols ());
384 MV B_mv (rcpFromRef (LocalMap), Bvalues, B.stride (), B.numCols ());
385 mv.multiply (Teuchos::NO_TRANS, Teuchos::NO_TRANS, alpha, A, B_mv, beta);
402 mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
405 static void MvScale (MV& mv, Scalar alpha) {
409 static void MvScale (MV& mv,
const std::vector<Scalar>& alphas) {
417 Teuchos::SerialDenseMatrix<int,Scalar>& C)
419 using Tpetra::LocallyReplicated;
423 using Teuchos::TimeMonitor;
424 typedef Tpetra::Map<LO,GO,Node> map_type;
426 #ifdef HAVE_ANASAZI_TPETRA_TIMERS 427 const std::string timerName (
"Anasazi::MVT::MvTransMv");
428 RCP<Teuchos::Time> timer = TimeMonitor::lookupCounter (timerName);
429 if (timer.is_null ()) {
430 timer = TimeMonitor::getNewCounter (timerName);
432 TEUCHOS_TEST_FOR_EXCEPTION(
433 timer.is_null (), std::logic_error,
"Anasazi::MvTransMv: " 434 "Failed to look up timer \"" << timerName <<
"\". " 435 "Please report this bug to the Belos developers.");
438 TimeMonitor timeMon (*timer);
439 #endif // HAVE_ANASAZI_TPETRA_TIMERS 449 const int numRowsC = C.numRows ();
450 const int numColsC = C.numCols ();
451 const int strideC = C.stride ();
454 if (numRowsC == 1 && numColsC == 1) {
455 if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
460 A.dot (B, Teuchos::ArrayView<Scalar> (C.values (), 1));
461 if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
468 RCP<const Comm<int> > pcomm = A.getMap ()->getComm ();
471 RCP<const map_type> LocalMap =
472 rcp (
new map_type (numRowsC, 0, pcomm, LocallyReplicated,
473 A.getMap ()->getNode ()));
475 const bool INIT_TO_ZERO =
true;
476 MV C_mv (LocalMap, numColsC, INIT_TO_ZERO);
479 C_mv.multiply (Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, alpha, A, B,
480 Teuchos::ScalarTraits<Scalar>::zero ());
483 Teuchos::ArrayView<Scalar> C_view (C.values (), strideC * numColsC);
488 C_mv.get1dCopy (C_view, strideC);
493 MvDot (
const MV& A,
const MV& B, std::vector<Scalar> &dots)
495 const size_t numVecs = A.getNumVectors ();
496 TEUCHOS_TEST_FOR_EXCEPTION(
497 numVecs != B.getNumVectors (), std::invalid_argument,
498 "Anasazi::MultiVecTraits::MvDot(A,B,dots): " 499 "A and B must have the same number of columns. " 500 "A has " << numVecs <<
" column(s), " 501 "but B has " << B.getNumVectors () <<
" column(s).");
502 #ifdef HAVE_TPETRA_DEBUG 503 TEUCHOS_TEST_FOR_EXCEPTION(
504 dots.size() < numVecs, std::invalid_argument,
505 "Anasazi::MultiVecTraits::MvDot(A,B,dots): " 506 "The output array 'dots' must have room for all dot products. " 507 "A and B each have " << numVecs <<
" column(s), " 508 "but 'dots' only has " << dots.size() <<
" entry(/ies).");
509 #endif // HAVE_TPETRA_DEBUG 511 Teuchos::ArrayView<Scalar> av (dots);
512 A.dot (B, av (0, numVecs));
518 std::vector<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
520 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
521 #ifdef HAVE_TPETRA_DEBUG 522 TEUCHOS_TEST_FOR_EXCEPTION(
523 normvec.size () <
static_cast<std::vector<int>::size_type
> (mv.getNumVectors ()),
524 std::invalid_argument,
525 "Anasazi::MultiVecTraits::MvNorm(mv,normvec): The normvec output " 526 "argument must have at least as many entries as the number of vectors " 527 "(columns) in the MultiVector mv. normvec.size() = " << normvec.size ()
528 <<
" < mv.getNumVectors() = " << mv.getNumVectors () <<
".");
529 #endif // HAVE_TPETRA_DEBUG 530 Teuchos::ArrayView<magnitude_type> av (normvec);
531 mv.norm2 (av (0, mv.getNumVectors ()));
535 SetBlock (
const MV& A,
const std::vector<int>& index, MV& mv)
537 using Teuchos::Range1D;
539 const size_t inNumVecs = A.getNumVectors ();
540 #ifdef HAVE_TPETRA_DEBUG 541 TEUCHOS_TEST_FOR_EXCEPTION(
542 inNumVecs < static_cast<size_t> (index.size ()), std::invalid_argument,
543 "Anasazi::MultiVecTraits::SetBlock(A,index,mv): 'index' argument must " 544 "have no more entries as the number of columns in the input MultiVector" 545 " A. A.getNumVectors() = " << inNumVecs <<
" < index.size () = " 546 << index.size () <<
".");
547 #endif // HAVE_TPETRA_DEBUG 549 if (inNumVecs > static_cast<size_t> (index.size ())) {
550 RCP<const MV> Asub = A.subView (Range1D (0, index.size () - 1));
551 Tpetra::deep_copy (*mvsub, *Asub);
553 Tpetra::deep_copy (*mvsub, A);
558 SetBlock (
const MV& A,
const Teuchos::Range1D& index, MV& mv)
567 const size_t maxInt =
568 static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
569 const bool overflow =
570 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
572 std::ostringstream os;
573 os <<
"Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
574 <<
", " << index.ubound () <<
"], mv): ";
575 TEUCHOS_TEST_FOR_EXCEPTION(
576 maxInt < A.getNumVectors (), std::range_error, os.str () <<
"Number " 577 "of columns (size_t) in the input MultiVector 'A' overflows int.");
578 TEUCHOS_TEST_FOR_EXCEPTION(
579 maxInt < mv.getNumVectors (), std::range_error, os.str () <<
"Number " 580 "of columns (size_t) in the output MultiVector 'mv' overflows int.");
583 const int numColsA =
static_cast<int> (A.getNumVectors ());
584 const int numColsMv =
static_cast<int> (mv.getNumVectors ());
586 const bool validIndex =
587 index.lbound () >= 0 && index.ubound () < numColsMv;
589 const bool validSource = index.size () <= numColsA;
591 if (! validIndex || ! validSource) {
592 std::ostringstream os;
593 os <<
"Anasazi::MultiVecTraits::SetBlock(A, index=[" << index.lbound ()
594 <<
", " << index.ubound () <<
"], mv): ";
595 TEUCHOS_TEST_FOR_EXCEPTION(
596 index.lbound() < 0, std::invalid_argument,
597 os.str() <<
"Range lower bound must be nonnegative.");
598 TEUCHOS_TEST_FOR_EXCEPTION(
599 index.ubound() >= numColsMv, std::invalid_argument,
600 os.str() <<
"Range upper bound must be less than the number of " 601 "columns " << numColsA <<
" in the 'mv' output argument.");
602 TEUCHOS_TEST_FOR_EXCEPTION(
603 index.size() > numColsA, std::invalid_argument,
604 os.str() <<
"Range must have no more elements than the number of " 605 "columns " << numColsA <<
" in the 'A' input argument.");
606 TEUCHOS_TEST_FOR_EXCEPTION(
607 true, std::logic_error,
"Should never get here!");
613 Teuchos::RCP<MV> mv_view;
614 if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
615 mv_view = Teuchos::rcpFromRef (mv);
623 Teuchos::RCP<const MV> A_view;
624 if (index.size () == numColsA) {
625 A_view = Teuchos::rcpFromRef (A);
627 A_view =
CloneView (A, Teuchos::Range1D (0, index.size () - 1));
630 Tpetra::deep_copy (*mv_view, *A_view);
633 static void Assign (
const MV& A, MV& mv)
635 const char errPrefix[] =
"Anasazi::MultiVecTraits::Assign(A, mv): ";
644 const size_t maxInt =
645 static_cast<size_t> (Teuchos::OrdinalTraits<int>::max ());
646 const bool overflow =
647 maxInt < A.getNumVectors () && maxInt < mv.getNumVectors ();
649 TEUCHOS_TEST_FOR_EXCEPTION(
650 maxInt < A.getNumVectors(), std::range_error,
651 errPrefix <<
"Number of columns in the input multivector 'A' " 652 "(a size_t) overflows int.");
653 TEUCHOS_TEST_FOR_EXCEPTION(
654 maxInt < mv.getNumVectors(), std::range_error,
655 errPrefix <<
"Number of columns in the output multivector 'mv' " 656 "(a size_t) overflows int.");
657 TEUCHOS_TEST_FOR_EXCEPTION(
658 true, std::logic_error,
"Should never get here!");
661 const int numColsA =
static_cast<int> (A.getNumVectors ());
662 const int numColsMv =
static_cast<int> (mv.getNumVectors ());
663 if (numColsA > numColsMv) {
664 TEUCHOS_TEST_FOR_EXCEPTION(
665 numColsA > numColsMv, std::invalid_argument,
666 errPrefix <<
"Input multivector 'A' has " << numColsA <<
" columns, " 667 "but output multivector 'mv' has only " << numColsMv <<
" columns.");
668 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Should never get here!");
670 if (numColsA == numColsMv) {
671 Tpetra::deep_copy (mv, A);
673 Teuchos::RCP<MV> mv_view =
675 Tpetra::deep_copy (*mv_view, A);
684 MvInit (MV& mv,
const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
686 mv.putScalar (alpha);
689 static void MvPrint (
const MV& mv, std::ostream& os) {
693 #ifdef HAVE_ANASAZI_TSQR 694 typedef Tpetra::TsqrAdaptor<Tpetra::MultiVector<Scalar, LO, GO, Node> > tsqr_adaptor_type;
697 #endif // HAVE_ANASAZI_TSQR 702 template <
class Scalar,
class LO,
class GO,
class Node>
704 Tpetra::MultiVector<Scalar,LO,GO,Node>,
705 Tpetra::Operator<Scalar,LO,GO,Node> >
709 Apply (
const Tpetra::Operator<Scalar,LO,GO,Node>& Op,
710 const Tpetra::MultiVector<Scalar,LO,GO,Node>& X,
711 Tpetra::MultiVector<Scalar,LO,GO,Node>& Y)
713 Op.apply (X, Y, Teuchos::NO_TRANS);
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
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.
static void MvPrint(const MV &mv, std::ostream &os)
Print the mv multi-vector to the os output stream.
static Teuchos::RCP< MV > CloneCopy(const MV &mv, const Teuchos::Range1D &index)
Create and return a deep copy of the given columns of mv.