42 #ifndef THYRA_MULTI_VECTOR_STD_OPS_HPP 43 #define THYRA_MULTI_VECTOR_STD_OPS_HPP 45 #include "Thyra_MultiVectorStdOps_decl.hpp" 46 #include "Thyra_VectorStdOps.hpp" 47 #include "Thyra_VectorSpaceBase.hpp" 48 #include "Thyra_VectorStdOps.hpp" 49 #include "Thyra_MultiVectorBase.hpp" 50 #include "Thyra_VectorBase.hpp" 51 #include "RTOpPack_ROpSum.hpp" 52 #include "RTOpPack_ROpDotProd.hpp" 53 #include "RTOpPack_ROpNorm1.hpp" 54 #include "RTOpPack_ROpNormInf.hpp" 55 #include "RTOpPack_TOpAssignVectors.hpp" 56 #include "RTOpPack_TOpAXPY.hpp" 57 #include "RTOpPack_TOpLinearCombination.hpp" 58 #include "RTOpPack_TOpScaleVector.hpp" 59 #include "Teuchos_Assert.hpp" 60 #include "Teuchos_Assert.hpp" 61 #include "Teuchos_as.hpp" 64 template<
class Scalar>
65 void Thyra::norms(
const MultiVectorBase<Scalar>& V,
66 const ArrayView<
typename ScalarTraits<Scalar>::magnitudeType> &
norms )
68 typedef Teuchos::ScalarTraits<Scalar> ST;
69 const int m = V.domain()->dim();
70 Array<Scalar> prods(m);
71 V.range()->scalarProds(V, V, prods());
72 for (
int j = 0; j < m; ++j )
73 norms[j] = ST::magnitude(ST::squareroot(prods[j]));
77 template<
class Scalar>
78 void Thyra::dots(
const MultiVectorBase<Scalar>& V1,
const MultiVectorBase<Scalar>& V2,
79 const ArrayView<Scalar> &
dots )
81 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
82 const int m = V1.domain()->dim();
83 RTOpPack::ROpDotProd<Scalar> dot_op;
84 Array<RCP<RTOpPack::ReductTarget> > rcp_dot_targs(m);
85 Array<Ptr<RTOpPack::ReductTarget> > dot_targs(m);
86 for(
int kc = 0; kc < m; ++kc ) {
87 rcp_dot_targs[kc] = dot_op.reduct_obj_create();
88 dot_targs[kc] = rcp_dot_targs[kc].ptr();
90 applyOp<Scalar>( dot_op, tuple(ptrInArg(V1), ptrInArg(V2)),
91 ArrayView<Ptr<MultiVectorBase<Scalar> > >(null),
93 for(
int kc = 0; kc < m; ++kc ) {
94 dots[kc] = dot_op(*dot_targs[kc]);
99 template<
class Scalar>
100 void Thyra::sums(
const MultiVectorBase<Scalar>& V,
const ArrayView<Scalar> &
sums )
102 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
103 const int m = V.domain()->dim();
104 RTOpPack::ROpSum<Scalar> sum_op;
105 Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
106 Array<Ptr<RTOpPack::ReductTarget> > op_targs(m);
107 for(
int kc = 0; kc < m; ++kc ) {
108 rcp_op_targs[kc] = sum_op.reduct_obj_create();
109 op_targs[kc] = rcp_op_targs[kc].ptr();
111 applyOp<Scalar>(sum_op, tuple(ptrInArg(V)),
112 ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null), op_targs);
113 for(
int kc = 0; kc < m; ++kc ) {
114 sums[kc] = sum_op(*op_targs[kc]);
119 template<
class Scalar>
120 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
121 Thyra::norm_1(
const MultiVectorBase<Scalar>& V )
123 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
125 RTOpPack::ROpNorm1<Scalar> sum_abs_op;
127 RTOpPack::ROpNormInf<Scalar> max_op;
129 RCP<RTOpPack::ReductTarget>
130 max_targ = max_op.reduct_obj_create();
132 Thyra::applyOp<Scalar>(sum_abs_op, max_op, tuple(ptrInArg(V))(),
133 ArrayView<
const Ptr<MultiVectorBase<Scalar> > >(null),
136 return max_op(*max_targ);
140 template<
class Scalar>
141 void Thyra::scale( Scalar alpha,
const Ptr<MultiVectorBase<Scalar> > &V )
143 using Teuchos::tuple;
using Teuchos::null;
144 typedef ScalarTraits<Scalar> ST;
145 if (alpha==ST::zero()) {
149 if (alpha==ST::one()) {
153 applyOp<Scalar>(scale_vector_op,
154 ArrayView<Ptr<const MultiVectorBase<Scalar> > >(null),
159 template<
class Scalar>
160 void Thyra::scaleUpdate(
const VectorBase<Scalar>& a,
161 const MultiVectorBase<Scalar>& U,
const Ptr<MultiVectorBase<Scalar> > &V )
164 bool is_compatible = U.range()->isCompatible(*a.space());
165 TEUCHOS_TEST_FOR_EXCEPTION(
166 !is_compatible, Exceptions::IncompatibleVectorSpaces,
167 "update(...), Error, U.range()->isCompatible(*a.space())==false" );
168 is_compatible = U.range()->isCompatible(*V->range());
169 TEUCHOS_TEST_FOR_EXCEPTION(
170 !is_compatible, Exceptions::IncompatibleVectorSpaces,
171 "update(...), Error, U.range()->isCompatible((V->range())==false" );
172 is_compatible = U.domain()->isCompatible(*V->domain());
173 TEUCHOS_TEST_FOR_EXCEPTION(
174 !is_compatible, Exceptions::IncompatibleVectorSpaces,
175 "update(...), Error, U.domain().isCompatible(V->domain())==false" );
177 const int m = U.domain()->dim();
178 for(
int j = 0; j < m; ++j ) {
179 ele_wise_prod<Scalar>( 1.0, a, *U.col(j), V->col(j).ptr() );
184 template<
class Scalar>
185 void Thyra::assign(
const Ptr<MultiVectorBase<Scalar> > &V, Scalar alpha )
191 template<
class Scalar>
192 void Thyra::assign(
const Ptr<MultiVectorBase<Scalar> > &V,
193 const MultiVectorBase<Scalar>& U )
195 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
196 RTOpPack::TOpAssignVectors<Scalar> assign_vectors_op;
197 applyOp<Scalar>( assign_vectors_op, tuple(ptrInArg(U)), tuple(V), null );
201 template<
class Scalar>
202 void Thyra::update( Scalar alpha,
const MultiVectorBase<Scalar>& U,
203 const Ptr<MultiVectorBase<Scalar> > &V )
205 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::null;
207 applyOp<Scalar>( axpy_op, tuple(ptrInArg(U)), tuple(V), null );
211 template<
class Scalar>
212 void Thyra::update(
const ArrayView<const Scalar> &alpha, Scalar beta,
213 const MultiVectorBase<Scalar>& U,
const Ptr<MultiVectorBase<Scalar> > &V )
216 bool is_compatible = U.range()->isCompatible(*V->range());
217 TEUCHOS_TEST_FOR_EXCEPTION(
218 !is_compatible, Exceptions::IncompatibleVectorSpaces,
219 "update(...), Error, U.range()->isCompatible((V->range())==false");
220 is_compatible = U.domain()->isCompatible(*V->domain());
221 TEUCHOS_TEST_FOR_EXCEPTION(
222 !is_compatible, Exceptions::IncompatibleVectorSpaces,
223 "update(...), Error, U.domain().isCompatible(V->domain())==false");
225 const int m = U.domain()->dim();
226 for(
int j = 0; j < m; ++j )
227 Vp_StV<Scalar>( V->col(j).ptr(), alpha[j]*beta, *U.col(j) );
231 template<
class Scalar>
232 void Thyra::update(
const MultiVectorBase<Scalar>& U,
233 const ArrayView<const Scalar> &alpha, Scalar beta,
234 const Ptr<MultiVectorBase<Scalar> > &V )
237 bool is_compatible = U.range()->isCompatible(*V->range());
238 TEUCHOS_TEST_FOR_EXCEPTION(
239 !is_compatible, Exceptions::IncompatibleVectorSpaces,
240 "update(...), Error, U.range()->isCompatible((V->range())==false");
241 is_compatible = U.domain()->isCompatible(*V->domain());
242 TEUCHOS_TEST_FOR_EXCEPTION(
243 !is_compatible, Exceptions::IncompatibleVectorSpaces,
244 "update(...), Error, U.domain().isCompatible(V->domain())==false");
246 const int m = U.domain()->dim();
247 for(
int j = 0; j < m; ++j ) {
248 Vt_S<Scalar>( V->col(j).ptr(), alpha[j]*beta );
249 Vp_StV<Scalar>( V->col(j).ptr(), 1.0, *U.col(j) );
254 template<
class Scalar>
255 void Thyra::linear_combination(
256 const ArrayView<const Scalar> &alpha,
257 const ArrayView<
const Ptr<
const MultiVectorBase<Scalar> > > &X,
259 const Ptr<MultiVectorBase<Scalar> > &Y
262 using Teuchos::tuple;
using Teuchos::null;
264 TEUCHOS_ASSERT_EQUALITY(alpha.size(),X.size());
266 const int m = alpha.size();
267 if ( beta == ScalarTraits<Scalar>::one() && m == 1 ) {
268 update( alpha[0], *X[0], Y );
276 Thyra::applyOp<Scalar>( lin_comb_op, X, tuple(Y), null );
280 template<
class Scalar>
281 void Thyra::randomize( Scalar l, Scalar u,
282 const Ptr<MultiVectorBase<Scalar> > &V )
284 const int m = V->domain()->dim();
285 for(
int j = 0; j < m; ++j )
291 template<
class Scalar>
292 void Thyra::Vt_S(
const Ptr<MultiVectorBase<Scalar> > &Z,
293 const Scalar& alpha )
295 const int m = Z->domain()->dim();
296 for(
int j = 0; j < m; ++j )
297 Vt_S( Z->col(j).ptr(), alpha );
302 template<
class Scalar>
303 void Thyra::Vp_S(
const Ptr<MultiVectorBase<Scalar> > &Z,
304 const Scalar& alpha )
306 const int m = Z->domain()->dim();
307 for(
int j = 0; j < m; ++j )
308 Vp_S( Z->col(j).ptr(), alpha );
313 template<
class Scalar>
314 void Thyra::Vp_V(
const Ptr<MultiVectorBase<Scalar> > &Z,
315 const MultiVectorBase<Scalar>& X )
317 using Teuchos::tuple;
using Teuchos::ptrInArg;
318 typedef Teuchos::ScalarTraits<Scalar> ST;
319 linear_combination<Scalar>( tuple(ST::one()), tuple(ptrInArg(X)),
324 template<
class Scalar>
325 void Thyra::V_VpV(
const Ptr<MultiVectorBase<Scalar> > &Z,
326 const MultiVectorBase<Scalar>& X,
const MultiVectorBase<Scalar>& Y )
328 using Teuchos::tuple;
using Teuchos::ptrInArg;
329 typedef Teuchos::ScalarTraits<Scalar> ST;
330 linear_combination<Scalar>(
331 tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
337 template<
class Scalar>
338 void Thyra::V_VmV(
const Ptr<MultiVectorBase<Scalar> > &Z,
339 const MultiVectorBase<Scalar>& X,
const MultiVectorBase<Scalar>& Y )
341 using Teuchos::tuple;
using Teuchos::ptrInArg;
using Teuchos::as;
342 typedef Teuchos::ScalarTraits<Scalar> ST;
343 linear_combination<Scalar>(
344 tuple(ST::one(), as<Scalar>(-ST::one())), tuple(ptrInArg(X), ptrInArg(Y)),
350 template<
class Scalar>
352 const Ptr<MultiVectorBase<Scalar> > &Z,
const Scalar &alpha,
353 const MultiVectorBase<Scalar>& X,
const MultiVectorBase<Scalar>& Y
356 using Teuchos::tuple;
using Teuchos::ptrInArg;
357 typedef Teuchos::ScalarTraits<Scalar> ST;
358 linear_combination<Scalar>(
359 tuple(alpha, ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
369 #define THYRA_MULTI_VECTOR_STD_OPS_INSTANT(SCALAR) \ 371 template void norms( const MultiVectorBase<SCALAR >& V, \ 372 const ArrayView<ScalarTraits<SCALAR >::magnitudeType> &norms ); \ 374 template void dots( const MultiVectorBase<SCALAR >& V1, const MultiVectorBase<SCALAR >& V2, \ 375 const ArrayView<SCALAR > &dots ); \ 377 template void sums( const MultiVectorBase<SCALAR >& V, const ArrayView<SCALAR > &sums ); \ 379 template Teuchos::ScalarTraits<SCALAR >::magnitudeType \ 380 norm_1( const MultiVectorBase<SCALAR >& V ); \ 382 template void scale( SCALAR alpha, const Ptr<MultiVectorBase<SCALAR > > &V ); \ 384 template void scaleUpdate( const VectorBase<SCALAR >& a, \ 385 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \ 387 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, SCALAR alpha ); \ 389 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, \ 390 const MultiVectorBase<SCALAR >& U ); \ 392 template void update( SCALAR alpha, const MultiVectorBase<SCALAR >& U, \ 393 const Ptr<MultiVectorBase<SCALAR > > &V ); \ 395 template void update( const ArrayView<const SCALAR > &alpha, SCALAR beta, \ 396 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \ 398 template void update( const MultiVectorBase<SCALAR >& U, \ 399 const ArrayView<const SCALAR > &alpha, SCALAR beta, \ 400 const Ptr<MultiVectorBase<SCALAR > > &V ); \ 402 template void linear_combination( \ 403 const ArrayView<const SCALAR > &alpha, \ 404 const ArrayView<const Ptr<const MultiVectorBase<SCALAR > > > &X, \ 405 const SCALAR &beta, \ 406 const Ptr<MultiVectorBase<SCALAR > > &Y \ 409 template void randomize( SCALAR l, SCALAR u, \ 410 const Ptr<MultiVectorBase<SCALAR > > &V ); \ 412 template void Vt_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \ 413 const SCALAR & alpha ); \ 415 template void Vp_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \ 416 const SCALAR & alpha ); \ 418 template void Vp_V( const Ptr<MultiVectorBase<SCALAR > > &Z, \ 419 const MultiVectorBase<SCALAR >& X ); \ 421 template void V_VpV( const Ptr<MultiVectorBase<SCALAR > > &Z, \ 422 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \ 424 template void V_VmV( const Ptr<MultiVectorBase<SCALAR > > &Z, \ 425 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \ 427 template void V_StVpV( \ 428 const Ptr<MultiVectorBase<SCALAR > > &Z, const SCALAR &alpha, \ 429 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y \ 433 #endif // THYRA_MULTI_VECTOR_STD_OPS_HPP
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(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 LinearOpBase< Scalar > > scale(const Scalar &scalar, const RCP< const LinearOpBase< Scalar > > &Op, const std::string &label="")
Build an implicit const scaled linear operator.
void norms(const MultiVectorBase< Scalar > &V, const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms)
Column-wise multi-vector natural norm.
void randomize(Scalar l, Scalar u, const Ptr< MultiVectorBase< Scalar > > &V)
Generate a random multi-vector with elements uniformly distributed elements.
void dots(const MultiVectorBase< Scalar > &V1, const MultiVectorBase< Scalar > &V2, const ArrayView< Scalar > &dots)
Multi-vector dot product.
void update(Scalar alpha, const MultiVectorBase< Scalar > &U, const Ptr< MultiVectorBase< Scalar > > &V)
alpha*U + V -> V.
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.