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_MultiVectorBase.hpp"
17#include "Thyra_VectorBase.hpp"
18#include "RTOpPack_ROpSum.hpp"
19#include "RTOpPack_ROpNorm1.hpp"
20#include "RTOpPack_ROpNormInf.hpp"
21#include "Teuchos_Assert.hpp"
22#include "Teuchos_as.hpp"
23
24
25template<class Scalar>
26void Thyra::norms( const MultiVectorBase<Scalar>& V,
27 const ArrayView<typename ScalarTraits<Scalar>::magnitudeType> &norms )
28{
29 typedef Teuchos::ScalarTraits<Scalar> ST;
30 const int m = V.domain()->dim();
31 Array<Scalar> prods(m);
32 V.range()->scalarProds(V, V, prods());
33 for ( int j = 0; j < m; ++j )
34 norms[j] = ST::magnitude(ST::squareroot(prods[j]));
35}
36
37
38template<class Scalar>
39void Thyra::dots( const MultiVectorBase<Scalar>& V1, const MultiVectorBase<Scalar>& V2,
40 const ArrayView<Scalar> &dots )
41{
42 V2.dots(V1, dots);
43}
44
45
46template<class Scalar>
47void Thyra::sums( const MultiVectorBase<Scalar>& V, const ArrayView<Scalar> &sums )
48{
49 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
50 const int m = V.domain()->dim();
51 RTOpPack::ROpSum<Scalar> sum_op;
52 Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
53 Array<Ptr<RTOpPack::ReductTarget> > op_targs(m);
54 for( int kc = 0; kc < m; ++kc ) {
55 rcp_op_targs[kc] = sum_op.reduct_obj_create();
56 op_targs[kc] = rcp_op_targs[kc].ptr();
57 }
58 applyOp<Scalar>(sum_op, tuple(ptrInArg(V)),
59 ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null), op_targs);
60 for( int kc = 0; kc < m; ++kc ) {
61 sums[kc] = sum_op(*op_targs[kc]);
62 }
63}
64
65
66template<class Scalar>
67typename Teuchos::ScalarTraits<Scalar>::magnitudeType
68Thyra::norm_1( const MultiVectorBase<Scalar>& V )
69{
70 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
71 // Primary column-wise reduction (sum of absolute values)
72 RTOpPack::ROpNorm1<Scalar> sum_abs_op;
73 // Secondary reduction (max over all columns = induced norm_1 matrix norm)
74 RTOpPack::ROpNormInf<Scalar> max_op;
75 // Reduction object (must be same for both sum_abs and max_targ objects)
76 RCP<RTOpPack::ReductTarget>
77 max_targ = max_op.reduct_obj_create();
78 // Perform the reductions
79 Thyra::applyOp<Scalar>(sum_abs_op, max_op, tuple(ptrInArg(V))(),
80 ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null),
81 max_targ.ptr());
82 // Return the final value
83 return max_op(*max_targ);
84}
85
86
87template<class Scalar>
88void Thyra::scale( Scalar alpha, const Ptr<MultiVectorBase<Scalar> > &V )
89{
90 V->scale(alpha);
91}
92
93
94template<class Scalar>
95void Thyra::scaleUpdate( const VectorBase<Scalar>& a,
96 const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
97{
98#ifdef TEUCHOS_DEBUG
99 bool is_compatible = U.range()->isCompatible(*a.space());
100 TEUCHOS_TEST_FOR_EXCEPTION(
101 !is_compatible, Exceptions::IncompatibleVectorSpaces,
102 "update(...), Error, U.range()->isCompatible(*a.space())==false" );
103 is_compatible = U.range()->isCompatible(*V->range());
104 TEUCHOS_TEST_FOR_EXCEPTION(
105 !is_compatible, Exceptions::IncompatibleVectorSpaces,
106 "update(...), Error, U.range()->isCompatible((V->range())==false" );
107 is_compatible = U.domain()->isCompatible(*V->domain());
108 TEUCHOS_TEST_FOR_EXCEPTION(
109 !is_compatible, Exceptions::IncompatibleVectorSpaces,
110 "update(...), Error, U.domain().isCompatible(V->domain())==false" );
111#endif
112 const int m = U.domain()->dim();
113 for( int j = 0; j < m; ++j ) {
114 ele_wise_prod<Scalar>( 1.0, a, *U.col(j), V->col(j).ptr() );
115 }
116}
117
118
119template<class Scalar>
120void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V, Scalar alpha )
121{
122 V->assign(alpha);
123}
124
125
126template<class Scalar>
127void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V,
128 const MultiVectorBase<Scalar>& U )
129{
130 V->assign(U);
131}
132
133
134template<class Scalar>
135void Thyra::update( Scalar alpha, const MultiVectorBase<Scalar>& U,
136 const Ptr<MultiVectorBase<Scalar> > &V )
137{
138 V->update(alpha, U);
139}
140
141
142template<class Scalar>
143void Thyra::update( const ArrayView<const Scalar> &alpha, Scalar beta,
144 const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
145{
146#ifdef TEUCHOS_DEBUG
147 bool is_compatible = U.range()->isCompatible(*V->range());
148 TEUCHOS_TEST_FOR_EXCEPTION(
149 !is_compatible, Exceptions::IncompatibleVectorSpaces,
150 "update(...), Error, U.range()->isCompatible((V->range())==false");
151 is_compatible = U.domain()->isCompatible(*V->domain());
152 TEUCHOS_TEST_FOR_EXCEPTION(
153 !is_compatible, Exceptions::IncompatibleVectorSpaces,
154 "update(...), Error, U.domain().isCompatible(V->domain())==false");
155#endif
156 const int m = U.domain()->dim();
157 for( int j = 0; j < m; ++j )
158 Vp_StV<Scalar>( V->col(j).ptr(), alpha[j]*beta, *U.col(j) );
159}
160
161
162template<class Scalar>
163void Thyra::update( const MultiVectorBase<Scalar>& U,
164 const ArrayView<const Scalar> &alpha, Scalar beta,
165 const Ptr<MultiVectorBase<Scalar> > &V )
166{
167#ifdef TEUCHOS_DEBUG
168 bool 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");
176#endif
177 const int m = U.domain()->dim();
178 for( int j = 0; j < m; ++j ) {
179 Vt_S<Scalar>( V->col(j).ptr(), alpha[j]*beta );
180 Vp_StV<Scalar>( V->col(j).ptr(), 1.0, *U.col(j) );
181 }
182}
183
184
185template<class Scalar>
186void Thyra::linear_combination(
187 const ArrayView<const Scalar> &alpha,
188 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &X,
189 const Scalar &beta,
190 const Ptr<MultiVectorBase<Scalar> > &Y
191 )
192{
193 Y->linear_combination(alpha, X, beta);
194}
195
196
197template<class Scalar>
198void Thyra::randomize( Scalar l, Scalar u,
199 const Ptr<MultiVectorBase<Scalar> > &V )
200{
201 const int m = V->domain()->dim();
202 for( int j = 0; j < m; ++j )
203 randomize( l, u, V->col(j).ptr() );
204 // Todo: call applyOp(...) directly!
205}
206
207
208template<class Scalar>
209void Thyra::Vt_S( const Ptr<MultiVectorBase<Scalar> > &Z,
210 const Scalar& alpha )
211{
212 Z->scale(alpha);
213}
214
215
216template<class Scalar>
217void Thyra::Vp_S( const Ptr<MultiVectorBase<Scalar> > &Z,
218 const Scalar& alpha )
219{
220 const int m = Z->domain()->dim();
221 for( int j = 0; j < m; ++j )
222 Vp_S( Z->col(j).ptr(), alpha );
223 // Todo: call applyOp(...) directly!
224}
225
226
227template<class Scalar>
228void Thyra::Vp_V( const Ptr<MultiVectorBase<Scalar> > &Z,
229 const MultiVectorBase<Scalar>& X )
230{
231 using Teuchos::tuple; using Teuchos::ptrInArg;
232 typedef Teuchos::ScalarTraits<Scalar> ST;
233 linear_combination<Scalar>( tuple(ST::one()), tuple(ptrInArg(X)),
234 ST::one(), Z );
235}
236
237
238template<class Scalar>
239void Thyra::V_VpV( const Ptr<MultiVectorBase<Scalar> > &Z,
240 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
241{
242 using Teuchos::tuple; using Teuchos::ptrInArg;
243 typedef Teuchos::ScalarTraits<Scalar> ST;
244 linear_combination<Scalar>(
245 tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
246 ST::zero(), Z
247 );
248}
249
250
251template<class Scalar>
252void Thyra::V_VmV( const Ptr<MultiVectorBase<Scalar> > &Z,
253 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
254{
255 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::as;
256 typedef Teuchos::ScalarTraits<Scalar> ST;
257 linear_combination<Scalar>(
258 tuple(ST::one(), as<Scalar>(-ST::one())), tuple(ptrInArg(X), ptrInArg(Y)),
259 ST::zero(), Z
260 );
261}
262
263
264template<class Scalar>
265void Thyra::V_StVpV(
266 const Ptr<MultiVectorBase<Scalar> > &Z, const Scalar &alpha,
267 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y
268 )
269{
270 using Teuchos::tuple; using Teuchos::ptrInArg;
271 typedef Teuchos::ScalarTraits<Scalar> ST;
272 linear_combination<Scalar>(
273 tuple(alpha, ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
274 ST::zero(), Z
275 );
276}
277
278
279//
280// Explicit instant macro
281//
282
283#define THYRA_MULTI_VECTOR_STD_OPS_INSTANT(SCALAR) \
284 \
285 template void norms( const MultiVectorBase<SCALAR >& V, \
286 const ArrayView<ScalarTraits<SCALAR >::magnitudeType> &norms ); \
287 \
288 template void dots( const MultiVectorBase<SCALAR >& V1, const MultiVectorBase<SCALAR >& V2, \
289 const ArrayView<SCALAR > &dots ); \
290 \
291 template void sums( const MultiVectorBase<SCALAR >& V, const ArrayView<SCALAR > &sums ); \
292 \
293 template Teuchos::ScalarTraits<SCALAR >::magnitudeType \
294 norm_1( const MultiVectorBase<SCALAR >& V ); \
295 \
296 template void scale( SCALAR alpha, const Ptr<MultiVectorBase<SCALAR > > &V ); \
297 \
298 template void scaleUpdate( const VectorBase<SCALAR >& a, \
299 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
300 \
301 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, SCALAR alpha ); \
302 \
303 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, \
304 const MultiVectorBase<SCALAR >& U ); \
305 \
306 template void update( SCALAR alpha, const MultiVectorBase<SCALAR >& U, \
307 const Ptr<MultiVectorBase<SCALAR > > &V ); \
308 \
309 template void update( const ArrayView<const SCALAR > &alpha, SCALAR beta, \
310 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
311 \
312 template void update( const MultiVectorBase<SCALAR >& U, \
313 const ArrayView<const SCALAR > &alpha, SCALAR beta, \
314 const Ptr<MultiVectorBase<SCALAR > > &V ); \
315 \
316 template void linear_combination( \
317 const ArrayView<const SCALAR > &alpha, \
318 const ArrayView<const Ptr<const MultiVectorBase<SCALAR > > > &X, \
319 const SCALAR &beta, \
320 const Ptr<MultiVectorBase<SCALAR > > &Y \
321 ); \
322 \
323 template void randomize( SCALAR l, SCALAR u, \
324 const Ptr<MultiVectorBase<SCALAR > > &V ); \
325 \
326 template void Vt_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
327 const SCALAR & alpha ); \
328 \
329 template void Vp_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
330 const SCALAR & alpha ); \
331 \
332 template void Vp_V( const Ptr<MultiVectorBase<SCALAR > > &Z, \
333 const MultiVectorBase<SCALAR >& X ); \
334 \
335 template void V_VpV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
336 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
337 \
338 template void V_VmV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
339 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
340 \
341 template void V_StVpV( \
342 const Ptr<MultiVectorBase<SCALAR > > &Z, const SCALAR &alpha, \
343 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y \
344 ); \
345
346
347#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)
Generated by doxygen 1.9.8