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>
28void Thyra::norms( const MultiVectorBase<Scalar>& V,
29 const ArrayView<typename ScalarTraits<Scalar>::magnitudeType> &norms )
30{
31 typedef Teuchos::ScalarTraits<Scalar> ST;
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>
41void Thyra::dots( const MultiVectorBase<Scalar>& V1, const MultiVectorBase<Scalar>& V2,
42 const ArrayView<Scalar> &dots )
43{
44 V2.dots(V1, dots);
45}
46
47
48template<class Scalar>
49void Thyra::sums( const MultiVectorBase<Scalar>& V, const ArrayView<Scalar> &sums )
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);
55 Array<Ptr<RTOpPack::ReductTarget> > 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>
69typename Teuchos::ScalarTraits<Scalar>::magnitudeType
70Thyra::norm_1( const MultiVectorBase<Scalar>& V )
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)
78 RCP<RTOpPack::ReductTarget>
79 max_targ = max_op.reduct_obj_create();
80 // Perform the reductions
81 Thyra::applyOp<Scalar>(sum_abs_op, max_op, tuple(ptrInArg(V))(),
82 ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null),
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>
97void Thyra::scaleUpdate( const VectorBase<Scalar>& a,
98 const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
99{
100#ifdef TEUCHOS_DEBUG
101 bool is_compatible = U.range()->isCompatible(*a.space());
102 TEUCHOS_TEST_FOR_EXCEPTION(
103 !is_compatible, Exceptions::IncompatibleVectorSpaces,
104 "update(...), Error, U.range()->isCompatible(*a.space())==false" );
105 is_compatible = U.range()->isCompatible(*V->range());
106 TEUCHOS_TEST_FOR_EXCEPTION(
107 !is_compatible, Exceptions::IncompatibleVectorSpaces,
108 "update(...), Error, U.range()->isCompatible((V->range())==false" );
109 is_compatible = U.domain()->isCompatible(*V->domain());
110 TEUCHOS_TEST_FOR_EXCEPTION(
111 !is_compatible, Exceptions::IncompatibleVectorSpaces,
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>
129void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V,
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,
146 const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
147{
148#ifdef TEUCHOS_DEBUG
149 bool is_compatible = U.range()->isCompatible(*V->range());
150 TEUCHOS_TEST_FOR_EXCEPTION(
151 !is_compatible, Exceptions::IncompatibleVectorSpaces,
152 "update(...), Error, U.range()->isCompatible((V->range())==false");
153 is_compatible = U.domain()->isCompatible(*V->domain());
154 TEUCHOS_TEST_FOR_EXCEPTION(
155 !is_compatible, Exceptions::IncompatibleVectorSpaces,
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>
165void Thyra::update( const MultiVectorBase<Scalar>& U,
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());
171 TEUCHOS_TEST_FOR_EXCEPTION(
172 !is_compatible, Exceptions::IncompatibleVectorSpaces,
173 "update(...), Error, U.range()->isCompatible((V->range())==false");
174 is_compatible = U.domain()->isCompatible(*V->domain());
175 TEUCHOS_TEST_FOR_EXCEPTION(
176 !is_compatible, Exceptions::IncompatibleVectorSpaces,
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>
188void Thyra::linear_combination(
189 const ArrayView<const Scalar> &alpha,
190 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &X,
191 const Scalar &beta,
192 const Ptr<MultiVectorBase<Scalar> > &Y
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>
211void Thyra::Vt_S( const Ptr<MultiVectorBase<Scalar> > &Z,
212 const Scalar& alpha )
213{
214 Z->scale(alpha);
215}
216
217
218template<class Scalar>
219void Thyra::Vp_S( const Ptr<MultiVectorBase<Scalar> > &Z,
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>
230void Thyra::Vp_V( const Ptr<MultiVectorBase<Scalar> > &Z,
231 const MultiVectorBase<Scalar>& X )
232{
233 using Teuchos::tuple; using Teuchos::ptrInArg;
234 typedef Teuchos::ScalarTraits<Scalar> ST;
235 linear_combination<Scalar>( tuple(ST::one()), tuple(ptrInArg(X)),
236 ST::one(), Z );
237}
238
239
240template<class Scalar>
241void Thyra::V_VpV( const Ptr<MultiVectorBase<Scalar> > &Z,
242 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
243{
244 using Teuchos::tuple; using Teuchos::ptrInArg;
245 typedef Teuchos::ScalarTraits<Scalar> ST;
246 linear_combination<Scalar>(
247 tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
248 ST::zero(), Z
249 );
250}
251
252
253template<class Scalar>
254void Thyra::V_VmV( const Ptr<MultiVectorBase<Scalar> > &Z,
255 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
256{
257 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::as;
258 typedef Teuchos::ScalarTraits<Scalar> ST;
259 linear_combination<Scalar>(
260 tuple(ST::one(), as<Scalar>(-ST::one())), tuple(ptrInArg(X), ptrInArg(Y)),
261 ST::zero(), Z
262 );
263}
264
265
266template<class Scalar>
267void Thyra::V_StVpV(
268 const Ptr<MultiVectorBase<Scalar> > &Z, const Scalar &alpha,
269 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y
270 )
271{
272 using Teuchos::tuple; using Teuchos::ptrInArg;
273 typedef Teuchos::ScalarTraits<Scalar> ST;
274 linear_combination<Scalar>(
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
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::as
TypeTo as(const TypeFrom &t)
Teuchos::dyn_cast
T_To & dyn_cast(T_From &from)
Generated by doxygen 1.9.8