Thyra Version of the Day
Loading...
Searching...
No Matches
Thyra_MultiVectorStdOps_def.hpp
1// @HEADER
2// ***********************************************************************
3//
4// Thyra: Interfaces and Support for Abstract Numerical Algorithms
5// Copyright (2004) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef THYRA_MULTI_VECTOR_STD_OPS_HPP
43#define THYRA_MULTI_VECTOR_STD_OPS_HPP
44
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_ROpNorm1.hpp"
53#include "RTOpPack_ROpNormInf.hpp"
54#include "Teuchos_Assert.hpp"
55#include "Teuchos_Assert.hpp"
56#include "Teuchos_as.hpp"
57
58
59template<class Scalar>
60void Thyra::norms( const MultiVectorBase<Scalar>& V,
61 const ArrayView<typename ScalarTraits<Scalar>::magnitudeType> &norms )
62{
64 const int m = V.domain()->dim();
65 Array<Scalar> prods(m);
66 V.range()->scalarProds(V, V, prods());
67 for ( int j = 0; j < m; ++j )
68 norms[j] = ST::magnitude(ST::squareroot(prods[j]));
69}
70
71
72template<class Scalar>
73void Thyra::dots( const MultiVectorBase<Scalar>& V1, const MultiVectorBase<Scalar>& V2,
74 const ArrayView<Scalar> &dots )
75{
76 V2.dots(V1, dots);
77}
78
79
80template<class Scalar>
81void Thyra::sums( const MultiVectorBase<Scalar>& V, const ArrayView<Scalar> &sums )
82{
83 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
84 const int m = V.domain()->dim();
85 RTOpPack::ROpSum<Scalar> sum_op;
86 Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
87 Array<Ptr<RTOpPack::ReductTarget> > op_targs(m);
88 for( int kc = 0; kc < m; ++kc ) {
89 rcp_op_targs[kc] = sum_op.reduct_obj_create();
90 op_targs[kc] = rcp_op_targs[kc].ptr();
91 }
92 applyOp<Scalar>(sum_op, tuple(ptrInArg(V)),
93 ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null), op_targs);
94 for( int kc = 0; kc < m; ++kc ) {
95 sums[kc] = sum_op(*op_targs[kc]);
96 }
97}
98
99
100template<class Scalar>
102Thyra::norm_1( const MultiVectorBase<Scalar>& V )
103{
104 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
105 // Primary column-wise reduction (sum of absolute values)
106 RTOpPack::ROpNorm1<Scalar> sum_abs_op;
107 // Secondary reduction (max over all columns = induced norm_1 matrix norm)
108 RTOpPack::ROpNormInf<Scalar> max_op;
109 // Reduction object (must be same for both sum_abs and max_targ objects)
110 RCP<RTOpPack::ReductTarget>
111 max_targ = max_op.reduct_obj_create();
112 // Perform the reductions
113 Thyra::applyOp<Scalar>(sum_abs_op, max_op, tuple(ptrInArg(V))(),
114 ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null),
115 max_targ.ptr());
116 // Return the final value
117 return max_op(*max_targ);
118}
119
120
121template<class Scalar>
122void Thyra::scale( Scalar alpha, const Ptr<MultiVectorBase<Scalar> > &V )
123{
124 V->scale(alpha);
125}
126
127
128template<class Scalar>
129void Thyra::scaleUpdate( const VectorBase<Scalar>& a,
130 const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
131{
132#ifdef TEUCHOS_DEBUG
133 bool is_compatible = U.range()->isCompatible(*a.space());
135 !is_compatible, Exceptions::IncompatibleVectorSpaces,
136 "update(...), Error, U.range()->isCompatible(*a.space())==false" );
137 is_compatible = U.range()->isCompatible(*V->range());
139 !is_compatible, Exceptions::IncompatibleVectorSpaces,
140 "update(...), Error, U.range()->isCompatible((V->range())==false" );
141 is_compatible = U.domain()->isCompatible(*V->domain());
143 !is_compatible, Exceptions::IncompatibleVectorSpaces,
144 "update(...), Error, U.domain().isCompatible(V->domain())==false" );
145#endif
146 const int m = U.domain()->dim();
147 for( int j = 0; j < m; ++j ) {
148 ele_wise_prod<Scalar>( 1.0, a, *U.col(j), V->col(j).ptr() );
149 }
150}
151
152
153template<class Scalar>
154void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V, Scalar alpha )
155{
156 V->assign(alpha);
157}
158
159
160template<class Scalar>
161void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V,
162 const MultiVectorBase<Scalar>& U )
163{
164 V->assign(U);
165}
166
167
168template<class Scalar>
169void Thyra::update( Scalar alpha, const MultiVectorBase<Scalar>& U,
170 const Ptr<MultiVectorBase<Scalar> > &V )
171{
172 V->update(alpha, U);
173}
174
175
176template<class Scalar>
177void Thyra::update( const ArrayView<const Scalar> &alpha, Scalar beta,
178 const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
179{
180#ifdef TEUCHOS_DEBUG
181 bool is_compatible = U.range()->isCompatible(*V->range());
183 !is_compatible, Exceptions::IncompatibleVectorSpaces,
184 "update(...), Error, U.range()->isCompatible((V->range())==false");
185 is_compatible = U.domain()->isCompatible(*V->domain());
187 !is_compatible, Exceptions::IncompatibleVectorSpaces,
188 "update(...), Error, U.domain().isCompatible(V->domain())==false");
189#endif
190 const int m = U.domain()->dim();
191 for( int j = 0; j < m; ++j )
192 Vp_StV<Scalar>( V->col(j).ptr(), alpha[j]*beta, *U.col(j) );
193}
194
195
196template<class Scalar>
197void Thyra::update( const MultiVectorBase<Scalar>& U,
198 const ArrayView<const Scalar> &alpha, Scalar beta,
199 const Ptr<MultiVectorBase<Scalar> > &V )
200{
201#ifdef TEUCHOS_DEBUG
202 bool is_compatible = U.range()->isCompatible(*V->range());
204 !is_compatible, Exceptions::IncompatibleVectorSpaces,
205 "update(...), Error, U.range()->isCompatible((V->range())==false");
206 is_compatible = U.domain()->isCompatible(*V->domain());
208 !is_compatible, Exceptions::IncompatibleVectorSpaces,
209 "update(...), Error, U.domain().isCompatible(V->domain())==false");
210#endif
211 const int m = U.domain()->dim();
212 for( int j = 0; j < m; ++j ) {
213 Vt_S<Scalar>( V->col(j).ptr(), alpha[j]*beta );
214 Vp_StV<Scalar>( V->col(j).ptr(), 1.0, *U.col(j) );
215 }
216}
217
218
219template<class Scalar>
220void Thyra::linear_combination(
221 const ArrayView<const Scalar> &alpha,
222 const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &X,
223 const Scalar &beta,
224 const Ptr<MultiVectorBase<Scalar> > &Y
225 )
226{
227 Y->linear_combination(alpha, X, beta);
228}
229
230
231template<class Scalar>
232void Thyra::randomize( Scalar l, Scalar u,
233 const Ptr<MultiVectorBase<Scalar> > &V )
234{
235 const int m = V->domain()->dim();
236 for( int j = 0; j < m; ++j )
237 randomize( l, u, V->col(j).ptr() );
238 // Todo: call applyOp(...) directly!
239}
240
241
242template<class Scalar>
243void Thyra::Vt_S( const Ptr<MultiVectorBase<Scalar> > &Z,
244 const Scalar& alpha )
245{
246 Z->scale(alpha);
247}
248
249
250template<class Scalar>
251void Thyra::Vp_S( const Ptr<MultiVectorBase<Scalar> > &Z,
252 const Scalar& alpha )
253{
254 const int m = Z->domain()->dim();
255 for( int j = 0; j < m; ++j )
256 Vp_S( Z->col(j).ptr(), alpha );
257 // Todo: call applyOp(...) directly!
258}
259
260
261template<class Scalar>
262void Thyra::Vp_V( const Ptr<MultiVectorBase<Scalar> > &Z,
263 const MultiVectorBase<Scalar>& X )
264{
265 using Teuchos::tuple; using Teuchos::ptrInArg;
267 linear_combination<Scalar>( tuple(ST::one()), tuple(ptrInArg(X)),
268 ST::one(), Z );
269}
270
271
272template<class Scalar>
273void Thyra::V_VpV( const Ptr<MultiVectorBase<Scalar> > &Z,
274 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
275{
276 using Teuchos::tuple; using Teuchos::ptrInArg;
278 linear_combination<Scalar>(
279 tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
280 ST::zero(), Z
281 );
282}
283
284
285template<class Scalar>
286void Thyra::V_VmV( const Ptr<MultiVectorBase<Scalar> > &Z,
287 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
288{
289 using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::as;
291 linear_combination<Scalar>(
292 tuple(ST::one(), as<Scalar>(-ST::one())), tuple(ptrInArg(X), ptrInArg(Y)),
293 ST::zero(), Z
294 );
295}
296
297
298template<class Scalar>
299void Thyra::V_StVpV(
300 const Ptr<MultiVectorBase<Scalar> > &Z, const Scalar &alpha,
301 const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y
302 )
303{
304 using Teuchos::tuple; using Teuchos::ptrInArg;
306 linear_combination<Scalar>(
307 tuple(alpha, ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
308 ST::zero(), Z
309 );
310}
311
312
313//
314// Explicit instant macro
315//
316
317#define THYRA_MULTI_VECTOR_STD_OPS_INSTANT(SCALAR) \
318 \
319 template void norms( const MultiVectorBase<SCALAR >& V, \
320 const ArrayView<ScalarTraits<SCALAR >::magnitudeType> &norms ); \
321 \
322 template void dots( const MultiVectorBase<SCALAR >& V1, const MultiVectorBase<SCALAR >& V2, \
323 const ArrayView<SCALAR > &dots ); \
324 \
325 template void sums( const MultiVectorBase<SCALAR >& V, const ArrayView<SCALAR > &sums ); \
326 \
327 template Teuchos::ScalarTraits<SCALAR >::magnitudeType \
328 norm_1( const MultiVectorBase<SCALAR >& V ); \
329 \
330 template void scale( SCALAR alpha, const Ptr<MultiVectorBase<SCALAR > > &V ); \
331 \
332 template void scaleUpdate( const VectorBase<SCALAR >& a, \
333 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
334 \
335 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, SCALAR alpha ); \
336 \
337 template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, \
338 const MultiVectorBase<SCALAR >& U ); \
339 \
340 template void update( SCALAR alpha, const MultiVectorBase<SCALAR >& U, \
341 const Ptr<MultiVectorBase<SCALAR > > &V ); \
342 \
343 template void update( const ArrayView<const SCALAR > &alpha, SCALAR beta, \
344 const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
345 \
346 template void update( const MultiVectorBase<SCALAR >& U, \
347 const ArrayView<const SCALAR > &alpha, SCALAR beta, \
348 const Ptr<MultiVectorBase<SCALAR > > &V ); \
349 \
350 template void linear_combination( \
351 const ArrayView<const SCALAR > &alpha, \
352 const ArrayView<const Ptr<const MultiVectorBase<SCALAR > > > &X, \
353 const SCALAR &beta, \
354 const Ptr<MultiVectorBase<SCALAR > > &Y \
355 ); \
356 \
357 template void randomize( SCALAR l, SCALAR u, \
358 const Ptr<MultiVectorBase<SCALAR > > &V ); \
359 \
360 template void Vt_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
361 const SCALAR & alpha ); \
362 \
363 template void Vp_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
364 const SCALAR & alpha ); \
365 \
366 template void Vp_V( const Ptr<MultiVectorBase<SCALAR > > &Z, \
367 const MultiVectorBase<SCALAR >& X ); \
368 \
369 template void V_VpV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
370 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
371 \
372 template void V_VmV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
373 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
374 \
375 template void V_StVpV( \
376 const Ptr<MultiVectorBase<SCALAR > > &Z, const SCALAR &alpha, \
377 const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y \
378 ); \
379
380
381#endif // THYRA_MULTI_VECTOR_STD_OPS_HPP
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TypeTo as(const TypeFrom &t)
T_To & dyn_cast(T_From &from)