52 #ifndef BELOS_THYRA_ADAPTER_HPP 53 #define BELOS_THYRA_ADAPTER_HPP 55 #include "BelosConfigDefs.hpp" 56 #include "BelosMultiVecTraits.hpp" 57 #include "BelosOperatorTraits.hpp" 59 #include <Thyra_DetachedMultiVectorView.hpp> 60 #include <Thyra_MultiVectorBase.hpp> 61 #include <Thyra_MultiVectorStdOps.hpp> 62 #ifdef HAVE_BELOS_TSQR 64 #endif // HAVE_BELOS_TSQR 80 template<
class ScalarType>
81 class MultiVecTraits< ScalarType,
Thyra::MultiVectorBase<ScalarType> >
84 typedef Thyra::MultiVectorBase<ScalarType>
TMVB;
85 typedef Teuchos::ScalarTraits<ScalarType>
ST;
86 typedef typename ST::magnitudeType
magType;
97 static Teuchos::RCP<TMVB>
Clone(
const TMVB& mv,
const int numvecs )
99 Teuchos::RCP<TMVB> c = Thyra::createMembers( mv.range(), numvecs );
109 int numvecs = mv.domain()->dim();
111 Teuchos::RCP< TMVB > cc = Thyra::createMembers( mv.range(), numvecs );
113 Thyra::assign(cc.ptr(), mv);
122 static Teuchos::RCP<TMVB>
CloneCopy(
const TMVB& mv,
const std::vector<int>& index )
124 int numvecs = index.size();
126 Teuchos::RCP<TMVB> cc = Thyra::createMembers( mv.range(), numvecs );
128 Teuchos::RCP<const TMVB> view = mv.subView(index);
130 Thyra::assign(cc.ptr(), *view);
134 static Teuchos::RCP<TMVB>
137 const int numVecs = index.size();
139 Teuchos::RCP<TMVB> cc = Thyra::createMembers (mv.range(), numVecs);
141 Teuchos::RCP<const TMVB> view = mv.subView (index);
143 Thyra::assign (cc.ptr(), *view);
154 int numvecs = index.size();
170 for (
int i=0; i<numvecs; i++) {
171 if (lb+i != index[i]) contig =
false;
174 Teuchos::RCP< TMVB > cc;
176 const Thyra::Range1D rng(lb,lb+numvecs-1);
178 cc = mv.subView(rng);
182 cc = mv.subView(index);
187 static Teuchos::RCP<TMVB>
195 return mv.subView (index);
204 static Teuchos::RCP<const TMVB>
CloneView(
const TMVB& mv,
const std::vector<int>& index )
206 int numvecs = index.size();
222 for (
int i=0; i<numvecs; i++) {
223 if (lb+i != index[i]) contig =
false;
226 Teuchos::RCP< const TMVB > cc;
228 const Thyra::Range1D rng(lb,lb+numvecs-1);
230 cc = mv.subView(rng);
234 cc = mv.subView(index);
239 static Teuchos::RCP<const TMVB>
247 return mv.subView (index);
257 return Teuchos::as<ptrdiff_t>(mv.range()->dim());
262 {
return mv.domain()->dim(); }
272 const Teuchos::SerialDenseMatrix<int,ScalarType>& B,
273 const ScalarType beta,
TMVB& mv )
275 using Teuchos::arrayView;
using Teuchos::arcpFromArrayView;
276 const int m = B.numRows();
277 const int n = B.numCols();
279 Teuchos::RCP< const TMVB >
280 B_thyra = Thyra::createMembersView(
284 arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
288 Thyra::apply<ScalarType>(A, Thyra::NOTRANS, *B_thyra, Teuchos::outArg(mv), alpha, beta);
294 const ScalarType beta,
const TMVB& B,
TMVB& mv )
296 typedef Teuchos::ScalarTraits<ScalarType>
ST;
297 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::inoutArg;
299 Thyra::linear_combination<ScalarType>(
300 tuple(alpha, beta)(), tuple(ptrInArg(A), ptrInArg(B))(), ST::zero(), inoutArg(mv));
306 { Thyra::scale(alpha, Teuchos::inoutArg(mv)); }
310 static void MvScale (
TMVB& mv,
const std::vector<ScalarType>& alpha)
312 for (
unsigned int i=0; i<alpha.size(); i++) {
313 Thyra::scale<ScalarType> (alpha[i], mv.col(i).ptr());
320 Teuchos::SerialDenseMatrix<int,ScalarType>& B )
322 using Teuchos::arrayView;
using Teuchos::arcpFromArrayView;
324 int m = A.domain()->dim();
325 int n = mv.domain()->dim();
328 B_thyra = Thyra::createMembersView(
332 arcpFromArrayView(arrayView(&B(0,0), B.stride()*B.numCols())), B.stride()
335 Thyra::apply<ScalarType>(A, Thyra::CONJTRANS, mv, B_thyra.ptr(), alpha);
341 static void MvDot(
const TMVB& mv,
const TMVB& A, std::vector<ScalarType>& b )
342 { Thyra::dots(mv, A, Teuchos::arrayViewFromVector(b)); }
352 static void MvNorm(
const TMVB& mv, std::vector<magType>& normvec,
353 NormType type = TwoNorm ) {
355 Thyra::norms_2(mv, Teuchos::arrayViewFromVector(normvec));
356 else if(type == OneNorm)
357 Thyra::norms_1(mv, Teuchos::arrayViewFromVector(normvec));
358 else if(type == InfNorm)
359 Thyra::norms_inf(mv, Teuchos::arrayViewFromVector(normvec));
361 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
362 "Belos::MultiVecTraits::MvNorm (Thyra specialization): " 363 "invalid norm type. Must be either TwoNorm, OneNorm or InfNorm");
376 int numvecs = index.size();
377 std::vector<int> indexA(numvecs);
378 int numAcols = A.domain()->dim();
379 for (
int i=0; i<numvecs; i++) {
384 if ( numAcols < numvecs ) {
389 else if ( numAcols > numvecs ) {
391 indexA.resize( numAcols );
394 Teuchos::RCP< const TMVB > relsource = A.subView(indexA);
396 Teuchos::RCP< TMVB > reldest = mv.subView(index);
398 Thyra::assign(reldest.ptr(), *relsource);
404 const int numColsA = A.domain()->dim();
405 const int numColsMv = mv.domain()->dim();
407 const bool validIndex = index.lbound() >= 0 && index.ubound() < numColsMv;
409 const bool validSource = index.size() <= numColsA;
411 if (! validIndex || ! validSource)
413 std::ostringstream os;
414 os <<
"Belos::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar> " 415 ">::SetBlock(A, [" << index.lbound() <<
", " << index.ubound()
417 TEUCHOS_TEST_FOR_EXCEPTION(index.lbound() < 0, std::invalid_argument,
418 os.str() <<
"Range lower bound must be nonnegative.");
419 TEUCHOS_TEST_FOR_EXCEPTION(index.ubound() >= numColsMv, std::invalid_argument,
420 os.str() <<
"Range upper bound must be less than " 421 "the number of columns " << numColsA <<
" in the " 422 "'mv' output argument.");
423 TEUCHOS_TEST_FOR_EXCEPTION(index.size() > numColsA, std::invalid_argument,
424 os.str() <<
"Range must have no more elements than" 425 " the number of columns " << numColsA <<
" in the " 426 "'A' input argument.");
427 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Should never get here!");
433 Teuchos::RCP<TMVB> mv_view;
434 if (index.lbound() == 0 && index.ubound()+1 == numColsMv)
435 mv_view = Teuchos::rcpFromRef (mv);
437 mv_view = mv.subView (index);
442 Teuchos::RCP<const TMVB> A_view;
443 if (index.size() == numColsA)
444 A_view = Teuchos::rcpFromRef (A);
446 A_view = A.subView (Teuchos::Range1D(0, index.size()-1));
449 Thyra::assign(mv_view.ptr(), *A_view);
455 const int numColsA = A.domain()->dim();
456 const int numColsMv = mv.domain()->dim();
457 if (numColsA > numColsMv)
459 std::ostringstream os;
460 os <<
"Belos::MultiVecTraits<Scalar, Thyra::MultiVectorBase<Scalar>" 461 " >::Assign(A, mv): ";
462 TEUCHOS_TEST_FOR_EXCEPTION(numColsA > numColsMv, std::invalid_argument,
463 os.str() <<
"Input multivector 'A' has " 464 << numColsA <<
" columns, but output multivector " 465 "'mv' has only " << numColsMv <<
" columns.");
466 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Should never get here!");
469 if (numColsA == numColsMv) {
470 Thyra::assign (Teuchos::outArg (mv), A);
472 Teuchos::RCP<TMVB> mv_view =
473 CloneViewNonConst (mv, Teuchos::Range1D(0, numColsA-1));
474 Thyra::assign (mv_view.ptr(), A);
484 Thyra::randomize<ScalarType>(
485 -Teuchos::ScalarTraits<ScalarType>::one(),
486 Teuchos::ScalarTraits<ScalarType>::one(),
487 Teuchos::outArg(mv));
492 MvInit (
TMVB& mv, ScalarType alpha = Teuchos::ScalarTraits<ScalarType>::zero())
494 Thyra::assign (Teuchos::outArg (mv), alpha);
505 { os << describe(mv,Teuchos::VERB_EXTREME); }
509 #ifdef HAVE_BELOS_TSQR 516 #endif // HAVE_BELOS_TSQR 532 template<
class ScalarType>
533 class OperatorTraits <ScalarType,
534 Thyra::MultiVectorBase<ScalarType>,
535 Thyra::LinearOpBase<ScalarType> >
538 typedef Thyra::MultiVectorBase<ScalarType>
TMVB;
539 typedef Thyra::LinearOpBase<ScalarType>
TLOB;
561 ETrans trans = NOTRANS)
563 Thyra::EOpTransp whichOp;
571 whichOp = Thyra::NOTRANS;
572 else if (trans ==
TRANS)
573 whichOp = Thyra::TRANS;
575 whichOp = Thyra::CONJTRANS;
577 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
578 "Belos::OperatorTraits::Apply (Thyra specialization): " 579 "'trans' argument must be neither NOTRANS=" <<
NOTRANS 581 <<
", but instead has an invalid value of " << trans <<
".");
582 Thyra::apply<ScalarType>(Op, whichOp, x, Teuchos::outArg(y));
588 typedef Teuchos::ScalarTraits<ScalarType> STS;
597 return Op.opSupported (Thyra::TRANS) &&
598 (! STS::isComplex || Op.opSupported (Thyra::CONJTRANS));
Thyra::LinearOpBase< ScalarType > TLOB
static int GetNumberVecs(const TMVB &mv)
Obtain the number of vectors in mv.
Thyra::MultiVectorBase< ScalarType > TMVB
static ptrdiff_t GetGlobalLength(const TMVB &mv)
Obtain the std::vector length of mv.
static void MvPrint(const TMVB &mv, std::ostream &os)
Print the mv multi-std::vector to the os output stream.
Thyra::MultiVectorBase< ScalarType > TMVB
Stub adaptor from Thyra::MultiVectorBase to TSQR.
Teuchos::ScalarTraits< ScalarType > ST
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 void Assign(const TMVB &A, TMVB &mv)
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 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...
ST::magnitudeType magType
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 void MvScale(TMVB &mv, const ScalarType alpha)
Scale each element of the vectors in *this with alpha.
static void MvTimesMatAddMv(const ScalarType alpha, const TMVB &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, TMVB &mv)
Update mv with .
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 Teuchos::RCP< TMVB > CloneCopy(const TMVB &mv, const Teuchos::Range1D &index)
static void Apply(const TLOB &Op, const TMVB &x, TMVB &y, ETrans trans=NOTRANS)
Apply Op to x, storing the result in y.
static void MvRandom(TMVB &mv)
Replace the vectors in mv with random vectors.
static Teuchos::RCP< TMVB > CloneViewNonConst(TMVB &mv, const Teuchos::Range1D &index)
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 SetBlock(const TMVB &A, const Teuchos::Range1D &index, TMVB &mv)
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 bool HasApplyTranspose(const TLOB &Op)
Whether the operator implements applying the transpose.
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< const TMVB > CloneView(const TMVB &mv, const Teuchos::Range1D &index)
static Teuchos::RCP< TMVB > Clone(const TMVB &mv, const int numvecs)
Creates a new empty MultiVectorBase containing numvecs columns.
static void MvAddMv(const ScalarType alpha, const TMVB &A, const ScalarType beta, const TMVB &B, TMVB &mv)
Replace mv with .