118 typedef Tpetra::MultiVector<Scalar, LO, GO, Node> MV;
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);
137 Teuchos::RCP<MV> X_copy (
new MV (X, Teuchos::Copy));
144 X_copy->setCopyOrView (Teuchos::View);
159 static Teuchos::RCP<MV>
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(
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.");
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];
184 Teuchos::RCP<MV> X_copy = mv.subCopy (columns ());
185 X_copy->setCopyOrView (Teuchos::View);
195 static Teuchos::RCP<MV>
198 const bool validRange = index.size() > 0 &&
199 index.lbound() >= 0 &&
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 "
212 TEUCHOS_TEST_FOR_EXCEPTION(
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!");
219 Teuchos::RCP<MV> X_copy = mv.subCopy (index);
220 X_copy->setCopyOrView (Teuchos::View);
224 static Teuchos::RCP<MV>
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(
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.");
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];
250 Teuchos::RCP<MV> X_view = mv.subViewNonConst (columns ());
251 X_view->setCopyOrView (Teuchos::View);
255 static Teuchos::RCP<MV>
261 const int numCols =
static_cast<int> (mv.getNumVectors());
262 const bool validRange = index.size() > 0 &&
263 index.lbound() >= 0 && index.ubound() < numCols;
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 "
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!");
283 Teuchos::RCP<MV> X_view = mv.subViewNonConst (index);
284 X_view->setCopyOrView (Teuchos::View);
288 static Teuchos::RCP<const MV>
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 ();
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.");
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];
314 Teuchos::RCP<const MV> X_view = mv.subView (columns);
315 Teuchos::rcp_const_cast<MV> (X_view)->setCopyOrView (Teuchos::View);
319 static Teuchos::RCP<const MV>
320 CloneView (
const MV& mv,
const Teuchos::Range1D& index)
325 const int numCols =
static_cast<int> (mv.getNumVectors());
326 const bool validRange = index.size() > 0 &&
327 index.lbound() >= 0 && index.ubound() < numCols;
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 "
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!");
344 Teuchos::RCP<const MV> X_view = mv.subView (index);
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 ());
360 const Teuchos::SerialDenseMatrix<int, Scalar>& B,
364 using Teuchos::ArrayView;
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");
371 Teuchos::RCP<Teuchos::Time> timer =
372 Teuchos::TimeMonitor::lookupCounter (timerName);
373 if (timer.is_null ()) {
374 timer = Teuchos::TimeMonitor::getNewCounter (timerName);
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.");
383 Teuchos::TimeMonitor timeMon (*timer);
387 if (B.numRows () == 1 && B.numCols () == 1) {
388 mv.update (alpha*B(0,0), A, beta);
393 Teuchos::SerialComm<int> serialComm;
394 map_type LocalMap (B.numRows (), A.getMap ()->getIndexBase (),
395 rcpFromRef<
const Comm<int> > (serialComm),
396 Tpetra::LocallyReplicated);
398 ArrayView<const Scalar> Bvalues (B.values (), B.stride () * B.numCols ());
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);
418 mv.update (alpha, A, beta, B, Teuchos::ScalarTraits<Scalar>::zero ());
421 static void MvScale (MV& mv, Scalar alpha) {
425 static void MvScale (MV& mv,
const std::vector<Scalar>& alphas) {
433 Teuchos::SerialDenseMatrix<int,Scalar>& C)
435 using Tpetra::LocallyReplicated;
439 using Teuchos::TimeMonitor;
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);
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.");
454 TimeMonitor timeMon (*timer);
465 const int numRowsC = C.numRows ();
466 const int numColsC = C.numCols ();
467 const int strideC = C.stride ();
470 if (numRowsC == 1 && numColsC == 1) {
471 if (alpha == Teuchos::ScalarTraits<Scalar>::zero ()) {
476 A.dot (B, Teuchos::ArrayView<Scalar> (C.values (), 1));
477 if (alpha != Teuchos::ScalarTraits<Scalar>::one ()) {
484 RCP<const Comm<int> > pcomm = A.getMap ()->getComm ();
487 RCP<const map_type> LocalMap =
488 rcp (
new map_type (numRowsC, 0, pcomm, LocallyReplicated));
490 const bool INIT_TO_ZERO =
true;
491 MV C_mv (LocalMap, numColsC, INIT_TO_ZERO);
494 C_mv.multiply (Teuchos::CONJ_TRANS, Teuchos::NO_TRANS, alpha, A, B,
495 Teuchos::ScalarTraits<Scalar>::zero ());
498 Teuchos::ArrayView<Scalar> C_view (C.values (), strideC * numColsC);
503 C_mv.get1dCopy (C_view, strideC);
508 MvDot (
const MV& A,
const MV& B, std::vector<Scalar> &dots)
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).");
526 Teuchos::ArrayView<Scalar> av (dots);
527 A.dot (B, av (0, numVecs));
533 std::vector<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType> &normvec)
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 () <<
".");
545 Teuchos::ArrayView<magnitude_type> av (normvec);
546 mv.norm2 (av (0, mv.getNumVectors ()));
550 SetBlock (
const MV& A,
const std::vector<int>& index, MV& mv)
552 using Teuchos::Range1D;
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 () <<
".");
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);
568 Tpetra::deep_copy (*mvsub, A);
573 SetBlock (
const MV& A,
const Teuchos::Range1D& index, MV& mv)
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 ();
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.");
598 const int numColsA =
static_cast<int> (A.getNumVectors ());
599 const int numColsMv =
static_cast<int> (mv.getNumVectors ());
601 const bool validIndex =
602 index.lbound () >= 0 && index.ubound () < numColsMv;
604 const bool validSource = index.size () <= numColsA;
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!");
628 Teuchos::RCP<MV> mv_view;
629 if (index.lbound () == 0 && index.ubound () + 1 == numColsMv) {
630 mv_view = Teuchos::rcpFromRef (mv);
638 Teuchos::RCP<const MV> A_view;
639 if (index.size () == numColsA) {
640 A_view = Teuchos::rcpFromRef (A);
642 A_view =
CloneView (A, Teuchos::Range1D (0, index.size () - 1));
645 Tpetra::deep_copy (*mv_view, *A_view);
648 static void Assign (
const MV& A, MV& mv)
650 const char errPrefix[] =
"Anasazi::MultiVecTraits::Assign(A, mv): ";
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 ();
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!");
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!");
685 if (numColsA == numColsMv) {
686 Tpetra::deep_copy (mv, A);
688 Teuchos::RCP<MV> mv_view =
690 Tpetra::deep_copy (*mv_view, A);
699 MvInit (MV& mv,
const Scalar alpha = Teuchos::ScalarTraits<Scalar>::zero ())
701 mv.putScalar (alpha);
704 static void MvPrint (
const MV& mv, std::ostream& os) {
708#ifdef HAVE_ANASAZI_TSQR
711 typedef Tpetra::TsqrAdaptor<Tpetra::MultiVector<Scalar, LO, GO, Node> > tsqr_adaptor_type;