52#ifndef BELOS_THYRA_ADAPTER_HPP
53#define BELOS_THYRA_ADAPTER_HPP
59#include <Thyra_DetachedMultiVectorView.hpp>
60#include <Thyra_MultiVectorBase.hpp>
61#include <Thyra_MultiVectorStdOps.hpp>
63# include <Thyra_TsqrAdaptor.hpp>
66#include <Teuchos_TimeMonitor.hpp>
82 template<
class ScalarType>
86 typedef Thyra::MultiVectorBase<ScalarType> TMVB;
87 typedef Teuchos::ScalarTraits<ScalarType> ST;
88 typedef typename ST::magnitudeType magType;
99 static Teuchos::RCP<TMVB>
Clone(
const TMVB& mv,
const int numvecs )
101 Teuchos::RCP<TMVB> c = Thyra::createMembers( mv.range(), numvecs );
111 int numvecs = mv.domain()->dim();
113 Teuchos::RCP< TMVB > cc = Thyra::createMembers( mv.range(), numvecs );
115 Thyra::assign(cc.ptr(), mv);
124 static Teuchos::RCP<TMVB>
CloneCopy(
const TMVB& mv,
const std::vector<int>& index )
126 int numvecs = index.size();
128 Teuchos::RCP<TMVB> cc = Thyra::createMembers( mv.range(), numvecs );
130 Teuchos::RCP<const TMVB> view = mv.subView(index);
132 Thyra::assign(cc.ptr(), *view);
136 static Teuchos::RCP<TMVB>
137 CloneCopy (
const TMVB& mv,
const Teuchos::Range1D& index)
139 const int numVecs = index.size();
141 Teuchos::RCP<TMVB> cc = Thyra::createMembers (mv.range(), numVecs);
143 Teuchos::RCP<const TMVB> view = mv.subView (index);
145 Thyra::assign (cc.ptr(), *view);
156 int numvecs = index.size();
172 for (
int i=0; i<numvecs; i++) {
173 if (lb+i != index[i]) contig =
false;
176 Teuchos::RCP< TMVB > cc;
178 const Thyra::Range1D rng(lb,lb+numvecs-1);
180 cc = mv.subView(rng);
184 cc = mv.subView(index);
189 static Teuchos::RCP<TMVB>
197 return mv.subView (index);
206 static Teuchos::RCP<const TMVB>
CloneView(
const TMVB& mv,
const std::vector<int>& index )
208 int numvecs = index.size();
224 for (
int i=0; i<numvecs; i++) {
225 if (lb+i != index[i]) contig =
false;
228 Teuchos::RCP< const TMVB > cc;
230 const Thyra::Range1D rng(lb,lb+numvecs-1);
232 cc = mv.subView(rng);
236 cc = mv.subView(index);
241 static Teuchos::RCP<const TMVB>
242 CloneView (
const TMVB& mv,
const Teuchos::Range1D& index)
249 return mv.subView (index);
259 return Teuchos::as<ptrdiff_t>(mv.range()->dim());
264 {
return mv.domain()->dim(); }
274 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
275 const ScalarType beta, TMVB& mv )
277 using Teuchos::arrayView;
using Teuchos::arcpFromArrayView;
279 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"Belos::MVT::MvTimesMatAddMv")));
281 const int m = B.numRows();
282 const int n = B.numCols();
283 auto vs = A.domain();
285 Teuchos::RCP< const TMVB >
286 B_thyra = vs->createCachedMembersView(
289 arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
293 Thyra::apply<ScalarType>(A, Thyra::NOTRANS, *B_thyra, Teuchos::outArg(mv), alpha, beta);
298 static void MvAddMv(
const ScalarType alpha,
const TMVB& A,
299 const ScalarType beta,
const TMVB& B, TMVB& mv )
301 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::inoutArg;
303 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"Belos::MVT::MvAddMv")));
305 Thyra::linear_combination<ScalarType>(
306 tuple(alpha, beta)(), tuple(ptrInArg(A), ptrInArg(B))(), Teuchos::ScalarTraits<ScalarType>::zero(), inoutArg(mv));
311 static void MvScale ( TMVB& mv,
const ScalarType alpha )
313 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"Belos::MVT::MvScale")));
315 Thyra::scale(alpha, Teuchos::inoutArg(mv));
320 static void MvScale (TMVB& mv,
const std::vector<ScalarType>& alpha)
322 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"Belos::MVT::MvScale")));
324 for (
unsigned int i=0; i<alpha.size(); i++) {
325 Thyra::scale<ScalarType> (alpha[i], mv.col(i).ptr());
331 static void MvTransMv(
const ScalarType alpha,
const TMVB& A,
const TMVB& mv,
332 Teuchos::SerialDenseMatrix<int,ScalarType>& B )
334 using Teuchos::arrayView;
using Teuchos::arcpFromArrayView;
336 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"Belos::MVT::MvTransMv")));
339 int m = A.domain()->dim();
340 int n = mv.domain()->dim();
341 auto vs = A.domain();
344 B_thyra = vs->createCachedMembersView(
347 arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
350 Thyra::apply<ScalarType>(A, Thyra::CONJTRANS, mv, B_thyra.ptr(), alpha);
356 static void MvDot(
const TMVB& mv,
const TMVB& A, std::vector<ScalarType>& b )
358 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"Belos::MVT::MvDot")));
360 Thyra::dots(mv, A, Teuchos::arrayViewFromVector(b));
371 static void MvNorm(
const TMVB& mv, std::vector<magType>& normvec,
373 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"Belos::MVT::MvNorm")));
376 Thyra::norms_2(mv, Teuchos::arrayViewFromVector(normvec));
378 Thyra::norms_1(mv, Teuchos::arrayViewFromVector(normvec));
380 Thyra::norms_inf(mv, Teuchos::arrayViewFromVector(normvec));
382 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
383 "Belos::MultiVecTraits::MvNorm (Thyra specialization): "
384 "invalid norm type. Must be either TwoNorm, OneNorm or InfNorm");
394 static void SetBlock(
const TMVB& A,
const std::vector<int>& index, TMVB& mv )
397 int numvecs = index.size();
398 std::vector<int> indexA(numvecs);
399 int numAcols = A.domain()->dim();
400 for (
int i=0; i<numvecs; i++) {
405 if ( numAcols < numvecs ) {
410 else if ( numAcols > numvecs ) {
412 indexA.resize( numAcols );
415 Teuchos::RCP< const TMVB > relsource = A.subView(indexA);
417 Teuchos::RCP< TMVB > reldest = mv.subView(index);
419 Thyra::assign(reldest.ptr(), *relsource);
423 SetBlock (
const TMVB& A,
const Teuchos::Range1D& index, TMVB& mv)
425 const int numColsA = A.domain()->dim();
426 const int numColsMv = mv.domain()->dim();
428 const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
430 const bool validSource = index.size() <= numColsA;
432 if (! validIndex || ! validSource)
434 std::ostringstream os;
435 os <<
"Belos::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar> "
436 ">::SetBlock(A, [" << index.lbound() <<
", " << index.ubound()
438 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
439 os.str() <<
"Range lower bound must be nonnegative.");
440 TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
441 os.str() <<
"Range upper bound must be less than "
442 "the number of columns " << numColsA <<
" in the "
443 "'mv' output argument.");
444 TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
445 os.str() <<
"Range must have no more elements than"
446 " the number of columns " << numColsA <<
" in the "
447 "'A' input argument.");
448 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Should never get here!");
454 Teuchos::RCP<TMVB> mv_view;
455 if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
456 mv_view = Teuchos::rcpFromRef (mv);
458 mv_view = mv.subView (index);
463 Teuchos::RCP<const TMVB> A_view;
464 if (index.size() == numColsA)
465 A_view = Teuchos::rcpFromRef (A);
467 A_view = A.subView (Teuchos::Range1D(0, index.size()-1));
470 Thyra::assign(mv_view.ptr(), *A_view);
474 Assign (
const TMVB& A, TMVB& mv)
476 Teuchos::TimeMonitor tM(*Teuchos::TimeMonitor::getNewTimer(std::string(
"Belos::MVT::Assign")));
478 const int numColsA = A.domain()->dim();
479 const int numColsMv = mv.domain()->dim();
480 if (numColsA > numColsMv)
482 std::ostringstream os;
483 os <<
"Belos::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar>"
484 " >::Assign(A, mv): ";
485 TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
486 os.str() <<
"Input multivector 'A' has "
487 << numColsA <<
" columns, but output multivector "
488 "'mv' has only " << numColsMv <<
" columns.");
489 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Should never get here!");
492 if (numColsA == numColsMv) {
493 Thyra::assign (Teuchos::outArg (mv), A);
495 Teuchos::RCP<TMVB> mv_view =
497 Thyra::assign (mv_view.ptr(), A);
507 Thyra::randomize<ScalarType>(
508 -Teuchos::ScalarTraits<ScalarType>::one(),
509 Teuchos::ScalarTraits<ScalarType>::one(),
510 Teuchos::outArg(mv));
515 MvInit (TMVB& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero())
517 Thyra::assign (Teuchos::outArg (mv), alpha);
527 static void MvPrint(
const TMVB& mv, std::ostream& os )
528 { os << describe(mv,Teuchos::VERB_EXTREME); }
532#ifdef HAVE_BELOS_TSQR
538 typedef Thyra::TsqrAdaptor< ScalarType > tsqr_adaptor_type;
555 template<
class ScalarType>
557 Thyra::MultiVectorBase<ScalarType>,
558 Thyra::LinearOpBase<ScalarType> >
561 typedef Thyra::MultiVectorBase<ScalarType> TMVB;
562 typedef Thyra::LinearOpBase<ScalarType> TLOB;
586 Thyra::EOpTransp whichOp;
594 whichOp = Thyra::NOTRANS;
595 else if (trans ==
TRANS)
596 whichOp = Thyra::TRANS;
598 whichOp = Thyra::CONJTRANS;
600 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
601 "Belos::OperatorTraits::Apply (Thyra specialization): "
602 "'trans' argument must be neither NOTRANS=" <<
NOTRANS
604 <<
", but instead has an invalid value of " << trans <<
".");
605 Thyra::apply<ScalarType>(Op, whichOp, x, Teuchos::outArg(y));
611 typedef Teuchos::ScalarTraits<ScalarType> STS;
620 return Op.opSupported (Thyra::TRANS) &&
621 (! STS::isComplex || Op.opSupported (Thyra::CONJTRANS));
static void MvPrint(const TMVB &mv, std::ostream &os)
Print the mv multi-std::vector to the os output stream.
static void SetBlock(const TMVB &A, const std::vector< int > &index, TMVB &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index.
static void MvAddMv(const ScalarType alpha, const TMVB &A, const ScalarType beta, const TMVB &B, TMVB &mv)
Replace mv with .
static int GetNumberVecs(const TMVB &mv)
Obtain the number of vectors in mv.
static Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv)
Creates a new MultiVectorBase and copies contents of mv into the new std::vector (deep copy).
static ptrdiff_t GetGlobalLength(const TMVB &mv)
Obtain the std::vector length of mv.
static void MvDot(const TMVB &mv, const TMVB &A, std::vector< ScalarType > &b)
Compute a std::vector b where the components are the individual dot-products of the i-th columns of A...
static Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv, const std::vector< int > &index)
Creates a new MultiVectorBase and copies the selected contents of mv into the new std::vector (deep c...
static Teuchos::RCP< TMVB > CloneViewNonConst(TMVB &mv, const std::vector< int > &index)
Creates a new MultiVectorBase that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< TMVB > Clone(const TMVB &mv, const int numvecs)
Creates a new empty MultiVectorBase containing numvecs columns.
static Teuchos::RCP< const TMVB > CloneView(const TMVB &mv, const std::vector< int > &index)
Creates a new const MultiVectorBase that shares the selected contents of mv (shallow copy).
static void MvScale(TMVB &mv, const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
static void MvTransMv(const ScalarType alpha, const TMVB &A, const TMVB &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
static void MvScale(TMVB &mv, const std::vector< ScalarType > &alpha)
Scale each element of the i-th vector in *this with alpha[i].
static void MvTimesMatAddMv(const ScalarType alpha, const TMVB &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, TMVB &mv)
Update mv with .
static void MvRandom(TMVB &mv)
Replace the vectors in mv with random vectors.
static void MvNorm(const TMVB &mv, std::vector< magType > &normvec, NormType type=TwoNorm)
Compute the 2-norm of each individual std::vector of mv. Upon return, normvec[i] holds the value of ,...
static void MvInit(TMVB &mv, ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
static void Assign(const MV &A, MV &mv)
static bool HasApplyTranspose(const TLOB &Op)
Whether the operator implements applying the transpose.
static void Apply(const TLOB &Op, const TMVB &x, TMVB &y, ETrans trans=NOTRANS)
Apply Op to x, storing the result in y.