42#ifndef THYRA_MULTI_VECTOR_DEFAULT_BASE_DEF_HPP
43#define THYRA_MULTI_VECTOR_DEFAULT_BASE_DEF_HPP
46#include "Thyra_MultiVectorDefaultBase_decl.hpp"
47#include "Thyra_LinearOpDefaultBase.hpp"
48#include "Thyra_MultiVectorStdOps.hpp"
49#include "Thyra_VectorSpaceFactoryBase.hpp"
50#include "Thyra_VectorSpaceBase.hpp"
51#include "Thyra_VectorBase.hpp"
52#include "Thyra_AssertOp.hpp"
53#include "Thyra_DefaultColumnwiseMultiVector.hpp"
54#include "RTOpPack_TOpAssignScalar.hpp"
55#include "RTOpPack_ROpDotProd.hpp"
56#include "RTOpPack_ROpNorm1.hpp"
57#include "RTOpPack_ROpNorm2.hpp"
58#include "RTOpPack_ROpNormInf.hpp"
59#include "RTOpPack_TOpAssignVectors.hpp"
60#include "RTOpPack_TOpAXPY.hpp"
61#include "RTOpPack_TOpLinearCombination.hpp"
62#include "RTOpPack_TOpScaleVector.hpp"
63#include "Teuchos_Workspace.hpp"
64#include "Teuchos_Assert.hpp"
65#include "Teuchos_as.hpp"
75RCP<MultiVectorBase<Scalar> >
79 &l_domain = *this->domain(),
80 &l_range = *this->range();
82 copy = createMembers(l_range,l_domain.
dim());
83 ::Thyra::assign<Scalar>(copy.ptr(), *
this);
96 using Teuchos::tuple;
using Teuchos::null;
98 Thyra::applyOp<Scalar>(assign_scalar_op,
103template<
class Scalar>
106 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
107 RTOpPack::TOpAssignVectors<Scalar> assign_vectors_op;
108 Thyra::applyOp<Scalar>(assign_vectors_op, tuple(ptrInArg(mv)),
113template<
class Scalar>
116 using Teuchos::tuple;
using Teuchos::null;
118 if (alpha==ST::zero()) {
119 this->assign(ST::zero());
123 if (alpha==ST::one()) {
127 Thyra::applyOp<Scalar>(scale_vector_op,
133template<
class Scalar>
139 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
141 Thyra::applyOp<Scalar>( axpy_op, tuple(ptrInArg(mv)),
146template<
class Scalar>
153 using Teuchos::tuple;
using Teuchos::null;
157 const int m = alpha.
size();
159 this->update(alpha[0], *mv[0]);
167 Thyra::applyOp<Scalar>(lin_comb_op, mv,
172template<
class Scalar>
177 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
178 const int m = this->domain()->dim();
179 RTOpPack::ROpNorm1<Scalar> norm_op;
182 for(
int kc = 0; kc < m; ++kc) {
183 rcp_op_targs[kc] = norm_op.reduct_obj_create();
184 op_targs[kc] = rcp_op_targs[kc].ptr();
186 Thyra::applyOp<Scalar>(norm_op,
190 for(
int kc = 0; kc < m; ++kc) {
191 norms[kc] = norm_op(*op_targs[kc]);
196template<
class Scalar>
201 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
202 const int m = this->domain()->dim();
206 for(
int kc = 0; kc < m; ++kc) {
208 op_targs[kc] = rcp_op_targs[kc].ptr();
210 Thyra::applyOp<Scalar>(norm_op,
214 for(
int kc = 0; kc < m; ++kc) {
215 norms[kc] = norm_op(*op_targs[kc]);
220template<
class Scalar>
226 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
227 const int m = this->domain()->dim();
228 RTOpPack::ROpDotProd<Scalar> dot_op;
231 for(
int kc = 0; kc < m; ++kc ) {
232 rcp_dot_targs[kc] = dot_op.reduct_obj_create();
233 dot_targs[kc] = rcp_dot_targs[kc].ptr();
235 Thyra::applyOp<Scalar>( dot_op,
239 for(
int kc = 0; kc < m; ++kc ) {
240 prods[kc] = dot_op(*dot_targs[kc]);
245template<
class Scalar>
250 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
251 const int m = this->domain()->dim();
252 RTOpPack::ROpNormInf<Scalar> norm_op;
255 for(
int kc = 0; kc < m; ++kc) {
256 rcp_op_targs[kc] = norm_op.reduct_obj_create();
257 op_targs[kc] = rcp_op_targs[kc].ptr();
259 Thyra::applyOp<Scalar>(norm_op,
263 for(
int kc = 0; kc < m; ++kc) {
264 norms[kc] = norm_op(*op_targs[kc]);
269template<
class Scalar>
279 const Range1D colRng = Teuchos::full_range(colRng_in,0,dimDomain-1);
280 if( colRng.
lbound() == 0 && as<Ordinal>(colRng.
ubound()) == dimDomain-1 )
282 if( colRng.
size() ) {
284 Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,colRng.
size());
286 col_vecs[j-colRng.
lbound()] = Teuchos::rcp_const_cast<VectorBase<Scalar> >(this->col(j));
293 return Teuchos::null;
297template<
class Scalar>
307 const Range1D colRng = Teuchos::full_range(colRng_in,0,dimDomain-1);
308 if( colRng.
lbound() == 0 && as<Ordinal>(colRng.
ubound()) == dimDomain-1 )
310 if( colRng.
size() ) {
312 Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,colRng.
size());
314 col_vecs[j-colRng.
lbound()] = this->col(j);
321 return Teuchos::null;
325template<
class Scalar>
334 const int numCols = cols.
size();
338 const char msg_err[] =
"MultiVectorDefaultBase<Scalar>::subView(numCols,cols[]): Error!";
342 Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,numCols);
343 for(
int k = 0; k < numCols; ++k ) {
344 const int col_k = cols[k];
347 !( 0 <= col_k && col_k < dimDomain ), std::invalid_argument
348 ,msg_err <<
" col["<<k<<
"] = " << col_k <<
" is not in the range [0,"<<(dimDomain-1)<<
"]!"
351 col_vecs[k] = Teuchos::rcp_const_cast<VectorBase<Scalar> >(this->col(col_k));
355 this->range(), l_range.
smallVecSpcFcty()->createVecSpc(numCols), col_vecs
361template<
class Scalar>
370 const int numCols = cols.
size();
374 const char msg_err[] =
"MultiVectorDefaultBase<Scalar>::subView(numCols,cols[]): Error!";
378 Workspace< RCP< VectorBase<Scalar> > > col_vecs(wss,numCols);
379 for(
int k = 0; k < numCols; ++k ) {
380 const int col_k = cols[k];
383 !( 0 <= col_k && col_k < dimDomain ), std::invalid_argument
384 ,msg_err <<
" col["<<k<<
"] = " << col_k <<
" is not in the range [0,"<<(dimDomain-1)<<
"]!"
387 col_vecs[k] = this->col(col_k);
391 this->range(), l_range.
smallVecSpcFcty()->createVecSpc(numCols), col_vecs
397template<
class Scalar>
403 const Ordinal prim_global_offset_in
411 const int num_multi_vecs = multi_vecs.size();
412 const int num_targ_multi_vecs = targ_multi_vecs.size();
427 Workspace<RCP<const VectorBase<Scalar> > > vecs_s(wss, num_multi_vecs);
428 Workspace<Ptr<const VectorBase<Scalar> > > vecs(wss, num_multi_vecs);
429 Workspace<RCP<VectorBase<Scalar> > > targ_vecs_s(wss, num_targ_multi_vecs);
430 Workspace<Ptr<VectorBase<Scalar> > > targ_vecs(wss, num_targ_multi_vecs);
432 for(
Ordinal j = 0; j < sec_dim; ++j) {
434 {
for(
Ordinal k = 0; k < as<Ordinal>(num_multi_vecs); ++k) {
435 vecs_s[k] = multi_vecs[k]->col(j);
436 vecs[k] = vecs_s[k].ptr();
438 {
for(
Ordinal k = 0; k < as<Ordinal>(num_targ_multi_vecs); ++k) {
439 targ_vecs_s[k] = targ_multi_vecs[k]->col(j);
440 targ_vecs[k] = targ_vecs_s[k].ptr();
446 targ_vecs().getConst(),
448 prim_global_offset_in);
456template<
class Scalar>
463 const Ordinal prim_global_offset_in
479 const int reduct_objs_size = (!is_null(reduct_obj) ? sec_dim : 0);
480 Workspace<RCP<RTOpPack::ReductTarget> > rcp_reduct_objs(wss, reduct_objs_size);
481 Workspace<Ptr<RTOpPack::ReductTarget> > reduct_objs(wss, reduct_objs_size);
482 if (!is_null(reduct_obj)) {
483 for(
Ordinal k = 0; k < sec_dim; ++k) {
485 reduct_objs[k] = rcp_reduct_objs[k].ptr();
491 prim_op, multi_vecs, targ_multi_vecs, reduct_objs,
492 prim_global_offset_in);
496 if (!is_null(reduct_obj)) {
497 for (
Ordinal k = 0; k < sec_dim; ++k) {
504template<
class Scalar>
512 rangeDim = this->range()->dim(),
513 domainDim = this->domain()->dim();
519 !(rowRng.
ubound() < rangeDim), std::out_of_range
520 ,
"MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(...): Error, rowRng = ["
521 <<rowRng.
lbound()<<
","<<rowRng.
ubound()<<
"] is not in the range = [0,"
525 !(colRng.ubound() < domainDim), std::out_of_range
526 ,
"MultiVectorDefaultBase<Scalar>::acquireDetachedMultiVectorViewImpl(...): Error, colRng = ["
527 <<colRng.lbound()<<
","<<colRng.ubound()<<
"] is not in the range = [0,"
528 <<(domainDim-1)<<
"]!"
535 for(
int k = colRng.lbound(); k <= colRng.ubound(); ++k ) {
537 col_k->acquireDetachedView( rowRng, &sv );
538 for(
int i = 0; i < rowRng.
size(); ++i )
539 values[ i + k*rowRng.
size() ] = sv[i];
540 col_k->releaseDetachedView( &sv );
554template<
class Scalar>
564template<
class Scalar>
585template<
class Scalar>
592 sub_mv==NULL, std::logic_error,
593 "MultiVectorDefaultBase<Scalar>::commitNonconstDetachedMultiVectorViewImpl(...): Error!"
599 for(
int k = sub_mv->
colOffset(); k < sub_mv->numSubCols(); ++k ) {
601 col_k->acquireDetachedView( rowRng, &msv );
602 for(
int i = 0; i < rowRng.
size(); ++i )
603 msv[i] = sub_mv->
values()[ i + k*rowRng.
size() ];
604 col_k->commitDetachedView( &msv );
Ordinal colOffset() const
void initialize(Ordinal globalOffset_in, Ordinal subDim_in, Ordinal colOffset_in, Ordinal numSubCols_in, const ArrayRCP< const Scalar > &values_in, Ordinal leadingDim_in)
Ordinal globalOffset() const
Teuchos::RCP< ReductTarget > reduct_obj_create() const
void reduce_reduct_objs(const ReductTarget &in_reduct_obj, const Ptr< ReductTarget > &inout_reduct_obj) const
const ArrayRCP< Scalar > values() const
Default subclass for MultiVectorBase implemented using columns of separate abstract vectors.
Interface for a collection of column vectors called a multi-vector.
virtual void dotsImpl(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
Default implementation of dots using RTOps.
virtual void norms2Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
Default implementation of norms_2 using RTOps.
virtual void linearCombinationImpl(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
Default implementation of linear_combination using RTOps.
virtual void releaseDetachedMultiVectorViewImpl(RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
virtual RCP< MultiVectorBase< Scalar > > clone_mv() const
RCP< const MultiVectorBase< Scalar > > contigSubViewImpl(const Range1D &colRng) const
virtual void mvSingleReductApplyOpImpl(const RTOpPack::RTOpT< Scalar > &primary_op, const RTOpPack::RTOpT< Scalar > &secondary_op, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &multi_vecs, const ArrayView< const Ptr< MultiVectorBase< Scalar > > > &targ_multi_vecs, const Ptr< RTOpPack::ReductTarget > &reduct_obj, const Ordinal primary_global_offset) const
virtual void acquireNonconstDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
virtual void normsInfImpl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
Default implementation of norms_inf using RTOps.
virtual void assignMultiVecImpl(const MultiVectorBase< Scalar > &mv)
Default implementation of assign(MV) using RTOps.
virtual void scaleImpl(Scalar alpha)
Default implementation of scale using RTOps.
virtual void norms1Impl(const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms) const
Default implementation of norms_1 using RTOps.
virtual void updateImpl(Scalar alpha, const MultiVectorBase< Scalar > &mv)
Default implementation of update using RTOps.
virtual void mvMultiReductApplyOpImpl(const RTOpPack::RTOpT< Scalar > &primary_op, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &multi_vecs, const ArrayView< const Ptr< MultiVectorBase< Scalar > > > &targ_multi_vecs, const ArrayView< const Ptr< RTOpPack::ReductTarget > > &reduct_objs, const Ordinal primary_global_offset) const
virtual void commitNonconstDetachedMultiVectorViewImpl(RTOpPack::SubMultiVectorView< Scalar > *sub_mv)
virtual void assignImpl(Scalar alpha)
Default implementation of assign(scalar) using RTOps.
virtual void acquireDetachedMultiVectorViewImpl(const Range1D &rowRng, const Range1D &colRng, RTOpPack::ConstSubMultiVectorView< Scalar > *sub_mv) const
RCP< const MultiVectorBase< Scalar > > nonContigSubViewImpl(const ArrayView< const int > &cols) const
RCP< MultiVectorBase< Scalar > > nonconstContigSubViewImpl(const Range1D &colRng)
RCP< MultiVectorBase< Scalar > > nonconstNonContigSubViewImpl(const ArrayView< const int > &cols)
Abstract interface for objects that represent a space for vectors.
virtual Ordinal dim() const =0
Return the dimension of the vector space.
virtual RCP< const VectorSpaceFactoryBase< Scalar > > smallVecSpcFcty() const =0
Return a VectorSpaceFactoryBase object for the creation of (usually serial) vector spaces with a smal...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
#define TEUCHOS_ASSERT_EQUALITY(val1, val2)
Teuchos::Ordinal Ordinal
Type for the dimension of a vector space. `*.
TypeTo as(const TypeFrom &t)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()