Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
27 template<class Scalar>
28 void 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 
40 template<class Scalar>
41 void Thyra::dots( const MultiVectorBase<Scalar>& V1, const MultiVectorBase<Scalar>& V2,
42  const ArrayView<Scalar> &dots )
43 {
44  V2.dots(V1, dots);
45 }
46 
47 
48 template<class Scalar>
49 void 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 
68 template<class Scalar>
69 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
70 Thyra::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 
89 template<class Scalar>
90 void Thyra::scale( Scalar alpha, const Ptr<MultiVectorBase<Scalar> > &V )
91 {
92  V->scale(alpha);
93 }
94 
95 
96 template<class Scalar>
97 void 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 
121 template<class Scalar>
122 void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V, Scalar alpha )
123 {
124  V->assign(alpha);
125 }
126 
127 
128 template<class Scalar>
129 void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V,
130  const MultiVectorBase<Scalar>& U )
131 {
132  V->assign(U);
133 }
134 
135 
136 template<class Scalar>
137 void Thyra::update( Scalar alpha, const MultiVectorBase<Scalar>& U,
138  const Ptr<MultiVectorBase<Scalar> > &V )
139 {
140  V->update(alpha, U);
141 }
142 
143 
144 template<class Scalar>
145 void 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 
164 template<class Scalar>
165 void 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 
187 template<class Scalar>
188 void 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 
199 template<class Scalar>
200 void 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 
210 template<class Scalar>
211 void Thyra::Vt_S( const Ptr<MultiVectorBase<Scalar> > &Z,
212  const Scalar& alpha )
213 {
214  Z->scale(alpha);
215 }
216 
217 
218 template<class Scalar>
219 void 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 
229 template<class Scalar>
230 void 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 
240 template<class Scalar>
241 void 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 
253 template<class Scalar>
254 void 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 
266 template<class Scalar>
267 void 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
Thyra::MultiVectorBase::Vp_S
void Vp_S(const Ptr< MultiVectorBase< Scalar > > &Z, const Scalar &alpha)
Z(i,j) += alpha, i = 0...Z-&gt;range()-&gt;dim()-1, j = 0...Z-&gt;domain()-&gt;dim()-1.
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Thyra::MultiVectorBase::sums
void sums(const MultiVectorBase< Scalar > &V, const ArrayView< Scalar > &sums)
Multi-vector column sum.
Thyra::MultiVectorBase::norms
void norms(const MultiVectorBase< Scalar > &V, const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms)
Column-wise multi-vector natural norm.
Thyra::MultiVectorBase::randomize
void randomize(Scalar l, Scalar u, const Ptr< MultiVectorBase< Scalar > > &V)
Generate a random multi-vector with elements uniformly distributed elements.
Teuchos::as
TypeTo as(const TypeFrom &t)
Thyra::MultiVectorBase::dots
void dots(const MultiVectorBase< Scalar > &V1, const MultiVectorBase< Scalar > &V2, const ArrayView< Scalar > &dots)
Multi-vector dot product.
Generated on Thu Nov 21 2024 09:20:52 for Thyra by   doxygen 1.8.5