Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_MultiVectorStdOps_def.hpp
1// @HEADER
2// *****************************************************************************
3// Thyra: Interfaces and Support for Abstract Numerical Algorithms
4//
5// Copyright 2004 NTESS and the Thyra contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef THYRA_MULTI_VECTOR_STD_OPS_HPP
11#define THYRA_MULTI_VECTOR_STD_OPS_HPP
12
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"
24#include "Teuchos_as.hpp"
25
26
27template<class Scalar>
30{
32 const int m = V.domain()->dim();
33 Array<Scalar> prods(m);
34 V.range()->scalarProds(V, V, prods());
35 for ( int j = 0; j < m; ++j )
36 norms[j] = ST::magnitude(ST::squareroot(prods[j]));
37}
38
39
40template<class Scalar>
42 const ArrayView<Scalar> &dots )
43{
44 V2.dots(V1, dots);
45}
46
47
48template<class Scalar>
50{
51 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
52 const int m = V.domain()->dim();
53 RTOpPack::ROpSum<Scalar> sum_op;
54 Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
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();
59 }
60 applyOp<Scalar>(sum_op, tuple(ptrInArg(V)),
61 ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null), op_targs);
62 for( int kc = 0; kc < m; ++kc ) {
63 sums[kc] = sum_op(*op_targs[kc]);
64 }
65}
66
67
68template<class Scalar>
71{
72 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
73 // Primary column-wise reduction (sum of absolute values)
74 RTOpPack::ROpNorm1<Scalar> sum_abs_op;
75 // Secondary reduction (max over all columns = induced norm_1 matrix norm)
76 RTOpPack::ROpNormInf<Scalar> max_op;
77 // Reduction object (must be same for both sum_abs and max_targ objects)
79 max_targ = max_op.reduct_obj_create();
80 // Perform the reductions
81 Thyra::applyOp<Scalar>(sum_abs_op, max_op, tuple(ptrInArg(V))(),
83 max_targ.ptr());
84 // Return the final value
85 return max_op(*max_targ);
86}
87
88
89template<class Scalar>
90void Thyra::scale( Scalar alpha, const Ptr<MultiVectorBase<Scalar> > &V )
91{
92 V->scale(alpha);
93}
94
95
96template<class Scalar>
99{
100#ifdef TEUCHOS_DEBUG
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" );
109 is_compatible = U.domain()->isCompatible(*V->domain());
112 "update(...), Error, U.domain().isCompatible(V->domain())==false" );
113#endif
114 const int m = U.domain()->dim();
115 for( int j = 0; j < m; ++j ) {
116 ele_wise_prod<Scalar>( 1.0, a, *U.col(j), V->col(j).ptr() );
117 }
118}
119
120
121template<class Scalar>
122void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V, Scalar alpha )
123{
124 V->assign(alpha);
125}
126
127
128template<class Scalar>
130 const MultiVectorBase<Scalar>& U )
131{
132 V->assign(U);
133}
134
135
136template<class Scalar>
137void Thyra::update( Scalar alpha, const MultiVectorBase<Scalar>& U,
138 const Ptr<MultiVectorBase<Scalar> > &V )
139{
140 V->update(alpha, U);
141}
142
143
144template<class Scalar>
145void Thyra::update( const ArrayView<const Scalar> &alpha, Scalar beta,
147{
148#ifdef TEUCHOS_DEBUG
149 bool is_compatible = U.range()->isCompatible(*V->range());
152 "update(...), Error, U.range()->isCompatible((V->range())==false");
153 is_compatible = U.domain()->isCompatible(*V->domain());
156 "update(...), Error, U.domain().isCompatible(V->domain())==false");
157#endif
158 const int m = U.domain()->dim();
159 for( int j = 0; j < m; ++j )
160 Vp_StV<Scalar>( V->col(j).ptr(), alpha[j]*beta, *U.col(j) );
161}
162
163
164template<class Scalar>
166 const ArrayView<const Scalar> &alpha, Scalar beta,
167 const Ptr<MultiVectorBase<Scalar> > &V )
168{
169#ifdef TEUCHOS_DEBUG
170 bool is_compatible = U.range()->isCompatible(*V->range());
173 "update(...), Error, U.range()->isCompatible((V->range())==false");
174 is_compatible = U.domain()->isCompatible(*V->domain());
177 "update(...), Error, U.domain().isCompatible(V->domain())==false");
178#endif
179 const int m = U.domain()->dim();
180 for( int j = 0; j < m; ++j ) {
181 Vt_S<Scalar>( V->col(j).ptr(), alpha[j]*beta );
182 Vp_StV<Scalar>( V->col(j).ptr(), 1.0, *U.col(j) );
183 }
184}
185
186
187template<class Scalar>
189 const ArrayView<const Scalar> &alpha,
190 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &X,
191 const Scalar &beta,
193 )
194{
195 Y->linear_combination(alpha, X, beta);
196}
197
198
199template<class Scalar>
200void Thyra::randomize( Scalar l, Scalar u,
201 const Ptr<MultiVectorBase<Scalar> > &V )
202{
203 const int m = V->domain()->dim();
204 for( int j = 0; j < m; ++j )
205 randomize( l, u, V->col(j).ptr() );
206 // Todo: call applyOp(...) directly!
207}
208
209
210template<class Scalar>
212 const Scalar& alpha )
213{
214 Z->scale(alpha);
215}
216
217
218template<class Scalar>
220 const Scalar& alpha )
221{
222 const int m = Z->domain()->dim();
223 for( int j = 0; j < m; ++j )
224 Vp_S( Z->col(j).ptr(), alpha );
225 // Todo: call applyOp(...) directly!
226}
227
228
229template<class Scalar>
231 const MultiVectorBase<Scalar>& X )
232{
233 using Teuchos::tuple; using Teuchos::ptrInArg;
235 linear_combination<Scalar>( tuple(ST::one()), tuple(ptrInArg(X)),
236 ST::one(), Z );
237}
238
239
240template<class Scalar>
243{
244 using Teuchos::tuple; using Teuchos::ptrInArg;
247 tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
248 ST::zero(), Z
249 );
250}
251
252
253template<class Scalar>
256{
257 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::as;
260 tuple(ST::one(), as<Scalar>(-ST::one())), tuple(ptrInArg(X), ptrInArg(Y)),
261 ST::zero(), Z
262 );
263}
264
265
266template<class Scalar>
268 const Ptr<MultiVectorBase<Scalar> > &Z, const Scalar &alpha,
270 )
271{
272 using Teuchos::tuple; using Teuchos::ptrInArg;
275 tuple(alpha, ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
276 ST::zero(), Z
277 );
278}
279
280
281//
282// Explicit instant macro
283//
284
285#define THYRA_MULTI_VECTOR_STD_OPS_INSTANT(SCALAR) \
286 \
287 template void norms( const MultiVectorBase<SCALAR >& V, \
288 const ArrayView<ScalarTraits<SCALAR >::magnitudeType> &norms ); \
289 \
290 template void dots( const MultiVectorBase<SCALAR >& V1, const MultiVectorBase<SCALAR >& V2, \
291 const ArrayView<SCALAR > &dots ); \
292 \
293 template void sums( const MultiVectorBase<SCALAR >& V, const ArrayView<SCALAR > &sums ); \
294 \
295 template Teuchos::ScalarTraits<SCALAR >::magnitudeType \
296 norm_1( const MultiVectorBase<SCALAR >& V ); \
297 \
298 template void scale( SCALAR alpha, const Ptr<MultiVectorBase<SCALAR > > &V ); \
299 \
300 template void scaleUpdate( const VectorBase<SCALAR >& a, \
301 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
302 \
303 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, SCALAR alpha ); \
304 \
305 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, \
306 const MultiVectorBase<SCALAR >& U ); \
307 \
308 template void update( SCALAR alpha, const MultiVectorBase<SCALAR >& U, \
309 const Ptr<MultiVectorBase<SCALAR > > &V ); \
310 \
311 template void update( const ArrayView<const SCALAR > &alpha, SCALAR beta, \
312 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
313 \
314 template void update( const MultiVectorBase<SCALAR >& U, \
315 const ArrayView<const SCALAR > &alpha, SCALAR beta, \
316 const Ptr<MultiVectorBase<SCALAR > > &V ); \
317 \
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 \
323 ); \
324 \
325 template void randomize( SCALAR l, SCALAR u, \
326 const Ptr<MultiVectorBase<SCALAR > > &V ); \
327 \
328 template void Vt_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
329 const SCALAR & alpha ); \
330 \
331 template void Vp_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
332 const SCALAR & alpha ); \
333 \
334 template void Vp_V( const Ptr<MultiVectorBase<SCALAR > > &Z, \
335 const MultiVectorBase<SCALAR >& X ); \
336 \
337 template void V_VpV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
338 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
339 \
340 template void V_VmV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
341 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
342 \
343 template void V_StVpV( \
344 const Ptr<MultiVectorBase<SCALAR > > &Z, const SCALAR &alpha, \
345 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y \
346 ); \
347
348
349#endif // THYRA_MULTI_VECTOR_STD_OPS_HPP
Ptr< T > ptr() const
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)