10#ifndef THYRA_MULTI_VECTOR_STD_OPS_HPP
11#define THYRA_MULTI_VECTOR_STD_OPS_HPP
13#include "Thyra_MultiVectorStdOps_decl.hpp"
14#include "Thyra_VectorStdOps.hpp"
15#include "Thyra_VectorSpaceBase.hpp"
16#include "Thyra_VectorStdOps.hpp"
17#include "Thyra_MultiVectorBase.hpp"
18#include "Thyra_VectorBase.hpp"
19#include "RTOpPack_ROpSum.hpp"
20#include "RTOpPack_ROpNorm1.hpp"
21#include "RTOpPack_ROpNormInf.hpp"
22#include "Teuchos_Assert.hpp"
23#include "Teuchos_Assert.hpp"
32 const int m = V.
domain()->dim();
34 V.
range()->scalarProds(V, V, prods());
35 for (
int j = 0; j < m; ++j )
36 norms[j] = ST::magnitude(ST::squareroot(prods[j]));
51 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
52 const int m = V.
domain()->dim();
53 RTOpPack::ROpSum<Scalar> sum_op;
56 for(
int kc = 0; kc < m; ++kc ) {
57 rcp_op_targs[kc] = sum_op.reduct_obj_create();
58 op_targs[kc] = rcp_op_targs[kc].ptr();
62 for(
int kc = 0; kc < m; ++kc ) {
63 sums[kc] = sum_op(*op_targs[kc]);
72 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
74 RTOpPack::ROpNorm1<Scalar> sum_abs_op;
76 RTOpPack::ROpNormInf<Scalar> max_op;
79 max_targ = max_op.reduct_obj_create();
85 return max_op(*max_targ);
101 bool is_compatible = U.
range()->isCompatible(*a.
space());
104 "update(...), Error, U.range()->isCompatible(*a.space())==false" );
105 is_compatible = U.
range()->isCompatible(*V->
range());
108 "update(...), Error, U.range()->isCompatible((V->range())==false" );
112 "update(...), Error, U.domain().isCompatible(V->domain())==false" );
114 const int m = U.
domain()->dim();
115 for(
int j = 0; j < m; ++j ) {
121template<
class Scalar>
128template<
class Scalar>
136template<
class Scalar>
144template<
class Scalar>
149 bool is_compatible = U.
range()->isCompatible(*V->
range());
152 "update(...), Error, U.range()->isCompatible((V->range())==false");
156 "update(...), Error, U.domain().isCompatible(V->domain())==false");
158 const int m = U.
domain()->dim();
159 for(
int j = 0; j < m; ++j )
164template<
class Scalar>
170 bool is_compatible = U.
range()->isCompatible(*V->
range());
173 "update(...), Error, U.range()->isCompatible((V->range())==false");
177 "update(...), Error, U.domain().isCompatible(V->domain())==false");
179 const int m = U.
domain()->dim();
180 for(
int j = 0; j < m; ++j ) {
187template<
class Scalar>
195 Y->linear_combination(alpha, X, beta);
199template<
class Scalar>
203 const int m = V->
domain()->dim();
204 for(
int j = 0; j < m; ++j )
210template<
class Scalar>
212 const Scalar& alpha )
218template<
class Scalar>
220 const Scalar& alpha )
222 const int m = Z->domain()->dim();
223 for(
int j = 0; j < m; ++j )
224 Vp_S( Z->col(j).ptr(), alpha );
229template<
class Scalar>
233 using Teuchos::tuple;
using Teuchos::ptrInArg;
240template<
class Scalar>
244 using Teuchos::tuple;
using Teuchos::ptrInArg;
247 tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
253template<
class Scalar>
257 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::as;
260 tuple(ST::one(), as<Scalar>(-ST::one())), tuple(ptrInArg(X), ptrInArg(Y)),
266template<
class Scalar>
272 using Teuchos::tuple;
using Teuchos::ptrInArg;
275 tuple(alpha, ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
285#define THYRA_MULTI_VECTOR_STD_OPS_INSTANT(SCALAR) \
287 template void norms( const MultiVectorBase<SCALAR >& V, \
288 const ArrayView<ScalarTraits<SCALAR >::magnitudeType> &norms ); \
290 template void dots( const MultiVectorBase<SCALAR >& V1, const MultiVectorBase<SCALAR >& V2, \
291 const ArrayView<SCALAR > &dots ); \
293 template void sums( const MultiVectorBase<SCALAR >& V, const ArrayView<SCALAR > &sums ); \
295 template Teuchos::ScalarTraits<SCALAR >::magnitudeType \
296 norm_1( const MultiVectorBase<SCALAR >& V ); \
298 template void scale( SCALAR alpha, const Ptr<MultiVectorBase<SCALAR > > &V ); \
300 template void scaleUpdate( const VectorBase<SCALAR >& a, \
301 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
303 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, SCALAR alpha ); \
305 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, \
306 const MultiVectorBase<SCALAR >& U ); \
308 template void update( SCALAR alpha, const MultiVectorBase<SCALAR >& U, \
309 const Ptr<MultiVectorBase<SCALAR > > &V ); \
311 template void update( const ArrayView<const SCALAR > &alpha, SCALAR beta, \
312 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
314 template void update( const MultiVectorBase<SCALAR >& U, \
315 const ArrayView<const SCALAR > &alpha, SCALAR beta, \
316 const Ptr<MultiVectorBase<SCALAR > > &V ); \
318 template void linear_combination( \
319 const ArrayView<const SCALAR > &alpha, \
320 const ArrayView<const Ptr<const MultiVectorBase<SCALAR > > > &X, \
321 const SCALAR &beta, \
322 const Ptr<MultiVectorBase<SCALAR > > &Y \
325 template void randomize( SCALAR l, SCALAR u, \
326 const Ptr<MultiVectorBase<SCALAR > > &V ); \
328 template void Vt_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
329 const SCALAR & alpha ); \
331 template void Vp_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
332 const SCALAR & alpha ); \
334 template void Vp_V( const Ptr<MultiVectorBase<SCALAR > > &Z, \
335 const MultiVectorBase<SCALAR >& X ); \
337 template void V_VpV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
338 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
340 template void V_VmV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
341 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
343 template void V_StVpV( \
344 const Ptr<MultiVectorBase<SCALAR > > &Z, const SCALAR &alpha, \
345 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y \
RCP< const LinearOpBase< Scalar > > scale(const Scalar &scalar, const RCP< const LinearOpBase< Scalar > > &Op, const std::string &label="")
Build an implicit const scaled linear operator.
Thrown if vector spaces are incompatible.
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
Return a smart pointer for the range space for this operator.
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
Return a smart pointer for the domain space for this operator.
Interface for a collection of column vectors called a multi-vector.
void assign(const Ptr< MultiVectorBase< Scalar > > &V, Scalar alpha)
V = alpha.
void sums(const MultiVectorBase< Scalar > &V, const ArrayView< Scalar > &sums)
Multi-vector column sum.
RCP< const VectorBase< Scalar > > col(Ordinal j) const
Calls colImpl().
void applyOp(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
Calls mvMultiReductApplyOpImpl().
void randomize(Scalar l, Scalar u, const Ptr< MultiVectorBase< Scalar > > &V)
Generate a random multi-vector with elements uniformly distributed elements.
void V_VpV(const Ptr< MultiVectorBase< Scalar > > &Z, const MultiVectorBase< Scalar > &X, const MultiVectorBase< Scalar > &Y)
Z(i,j) = X(i,j) + Y(i,j), i = 0...Z->range()->dim()-1, j = 0...Z->domain()->dim()-1.
void norms(const MultiVectorBase< Scalar > &V, const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms)
Column-wise multi-vector natural norm.
void dots(const MultiVectorBase< Scalar > &mv, const ArrayView< Scalar > &prods) const
Column-wise Euclidean dot product.
void applyOp(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=0)
Apply a reduction/transformation operator column by column and return an array of the reduction objec...
void Vt_S(const Ptr< MultiVectorBase< Scalar > > &Z, const Scalar &alpha)
Z(i,j) *= alpha, i = 0...Z->range()->dim()-1, j = 0...Z->domain()->dim()-1.
void update(Scalar alpha, const MultiVectorBase< Scalar > &U, const Ptr< MultiVectorBase< Scalar > > &V)
alpha*U + V -> V.
ScalarTraits< Scalar >::magnitudeType norm_1(const MultiVectorBase< Scalar > &V)
Take the induced matrix one norm of a multi-vector.
void V_VmV(const Ptr< MultiVectorBase< Scalar > > &Z, const MultiVectorBase< Scalar > &X, const MultiVectorBase< Scalar > &Y)
Z(i,j) = X(i,j) - Y(i,j), i = 0...Z->range()->dim()-1, j = 0...Z->domain()->dim()-1.
void dots(const MultiVectorBase< Scalar > &V1, const MultiVectorBase< Scalar > &V2, const ArrayView< Scalar > &dots)
Multi-vector dot product.
void linear_combination(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &mv, const Scalar &beta)
Y.col(j)(i) = beta*Y.col(j)(i) + sum( alpha[k]*X[k].col(j)(i),
void linear_combination(const ArrayView< const Scalar > &alpha, const ArrayView< const Ptr< const MultiVectorBase< Scalar > > > &X, const Scalar &beta, const Ptr< MultiVectorBase< Scalar > > &Y)
Y.col(j)(i) = beta*Y.col(j)(i) + sum( alpha[k]*X[k].col(j)(i), k=0...m-1 ), for i = 0....
void scaleUpdate(const VectorBase< Scalar > &a, const MultiVectorBase< Scalar > &U, const Ptr< MultiVectorBase< Scalar > > &V)
A*U + V -> V (where A is a diagonal matrix with diagonal a).
void update(Scalar alpha, const MultiVectorBase< Scalar > &mv)
void V_StVpV(const Ptr< MultiVectorBase< Scalar > > &Z, const Scalar &alpha, const MultiVectorBase< Scalar > &X, const MultiVectorBase< Scalar > &Y)
Z(i,j) = alpha*X(i,j) + Y(i), i = 0...z->space()->dim()-1, , j = 0...Z->domain()->dim()-1.
void Vp_V(const Ptr< MultiVectorBase< Scalar > > &Z, const MultiVectorBase< Scalar > &X)
Z(i,j) += X(i,j), i = 0...Z->range()->dim()-1, j = 0...Z->domain()->dim()-1.
void Vp_S(const Ptr< MultiVectorBase< Scalar > > &Z, const Scalar &alpha)
Z(i,j) += alpha, i = 0...Z->range()->dim()-1, j = 0...Z->domain()->dim()-1.
void assign(Scalar alpha)
V = alpha.
Abstract interface for finite-dimensional dense vectors.
void Vp_StV(const Ptr< VectorBase< Scalar > > &y, const Scalar &alpha, const VectorBase< Scalar > &x)
AXPY: y(i) = alpha * x(i) + y(i), i = 0...y->space()->dim()-1.
void ele_wise_prod(const Scalar &alpha, const VectorBase< Scalar > &x, const VectorBase< Scalar > &v, const Ptr< VectorBase< Scalar > > &y)
Element-wise product update: y(i) += alpha * x(i) * v(i), i = 0...y->space()->dim()-1.
virtual RCP< const VectorSpaceBase< Scalar > > space() const =0
Return a smart pointer to the vector space that this vector belongs to.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TypeTo as(const TypeFrom &t)