Thyra  Version of the Day
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_ROpDotProd.hpp"
53 #include "RTOpPack_ROpNorm1.hpp"
54 #include "RTOpPack_ROpNormInf.hpp"
55 #include "RTOpPack_TOpAssignVectors.hpp"
56 #include "RTOpPack_TOpAXPY.hpp"
57 #include "RTOpPack_TOpLinearCombination.hpp"
58 #include "RTOpPack_TOpScaleVector.hpp"
59 #include "Teuchos_Assert.hpp"
60 #include "Teuchos_Assert.hpp"
61 #include "Teuchos_as.hpp"
62 
63 
64 template<class Scalar>
65 void Thyra::norms( const MultiVectorBase<Scalar>& V,
66  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType> &norms )
67 {
68  typedef Teuchos::ScalarTraits<Scalar> ST;
69  const int m = V.domain()->dim();
70  Array<Scalar> prods(m);
71  V.range()->scalarProds(V, V, prods());
72  for ( int j = 0; j < m; ++j )
73  norms[j] = ST::magnitude(ST::squareroot(prods[j]));
74 }
75 
76 
77 template<class Scalar>
78 void Thyra::dots( const MultiVectorBase<Scalar>& V1, const MultiVectorBase<Scalar>& V2,
79  const ArrayView<Scalar> &dots )
80 {
81  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
82  const int m = V1.domain()->dim();
83  RTOpPack::ROpDotProd<Scalar> dot_op;
84  Array<RCP<RTOpPack::ReductTarget> > rcp_dot_targs(m);
85  Array<Ptr<RTOpPack::ReductTarget> > dot_targs(m);
86  for( int kc = 0; kc < m; ++kc ) {
87  rcp_dot_targs[kc] = dot_op.reduct_obj_create();
88  dot_targs[kc] = rcp_dot_targs[kc].ptr();
89  }
90  applyOp<Scalar>( dot_op, tuple(ptrInArg(V1), ptrInArg(V2)),
91  ArrayView<Ptr<MultiVectorBase<Scalar> > >(null),
92  dot_targs );
93  for( int kc = 0; kc < m; ++kc ) {
94  dots[kc] = dot_op(*dot_targs[kc]);
95  }
96 }
97 
98 
99 template<class Scalar>
100 void Thyra::sums( const MultiVectorBase<Scalar>& V, const ArrayView<Scalar> &sums )
101 {
102  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
103  const int m = V.domain()->dim();
104  RTOpPack::ROpSum<Scalar> sum_op;
105  Array<RCP<RTOpPack::ReductTarget> > rcp_op_targs(m);
106  Array<Ptr<RTOpPack::ReductTarget> > op_targs(m);
107  for( int kc = 0; kc < m; ++kc ) {
108  rcp_op_targs[kc] = sum_op.reduct_obj_create();
109  op_targs[kc] = rcp_op_targs[kc].ptr();
110  }
111  applyOp<Scalar>(sum_op, tuple(ptrInArg(V)),
112  ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null), op_targs);
113  for( int kc = 0; kc < m; ++kc ) {
114  sums[kc] = sum_op(*op_targs[kc]);
115  }
116 }
117 
118 
119 template<class Scalar>
120 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
121 Thyra::norm_1( const MultiVectorBase<Scalar>& V )
122 {
123  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
124  // Primary column-wise reduction (sum of absolute values)
125  RTOpPack::ROpNorm1<Scalar> sum_abs_op;
126  // Secondary reduction (max over all columns = induced norm_1 matrix norm)
127  RTOpPack::ROpNormInf<Scalar> max_op;
128  // Reduction object (must be same for both sum_abs and max_targ objects)
129  RCP<RTOpPack::ReductTarget>
130  max_targ = max_op.reduct_obj_create();
131  // Perform the reductions
132  Thyra::applyOp<Scalar>(sum_abs_op, max_op, tuple(ptrInArg(V))(),
133  ArrayView<const Ptr<MultiVectorBase<Scalar> > >(null),
134  max_targ.ptr());
135  // Return the final value
136  return max_op(*max_targ);
137 }
138 
139 
140 template<class Scalar>
141 void Thyra::scale( Scalar alpha, const Ptr<MultiVectorBase<Scalar> > &V )
142 {
143  using Teuchos::tuple; using Teuchos::null;
144  typedef ScalarTraits<Scalar> ST;
145  if (alpha==ST::zero()) {
146  assign( V, ST::zero() );
147  return;
148  }
149  if (alpha==ST::one()) {
150  return;
151  }
152  RTOpPack::TOpScaleVector<Scalar> scale_vector_op(alpha);
153  applyOp<Scalar>(scale_vector_op,
154  ArrayView<Ptr<const MultiVectorBase<Scalar> > >(null),
155  tuple(V), null );
156 }
157 
158 
159 template<class Scalar>
160 void Thyra::scaleUpdate( const VectorBase<Scalar>& a,
161  const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
162 {
163 #ifdef TEUCHOS_DEBUG
164  bool is_compatible = U.range()->isCompatible(*a.space());
165  TEUCHOS_TEST_FOR_EXCEPTION(
166  !is_compatible, Exceptions::IncompatibleVectorSpaces,
167  "update(...), Error, U.range()->isCompatible(*a.space())==false" );
168  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  ele_wise_prod<Scalar>( 1.0, a, *U.col(j), V->col(j).ptr() );
180  }
181 }
182 
183 
184 template<class Scalar>
185 void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V, Scalar alpha )
186 {
187  V->assign(alpha);
188 }
189 
190 
191 template<class Scalar>
192 void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V,
193  const MultiVectorBase<Scalar>& U )
194 {
195  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
196  RTOpPack::TOpAssignVectors<Scalar> assign_vectors_op;
197  applyOp<Scalar>( assign_vectors_op, tuple(ptrInArg(U)), tuple(V), null );
198 }
199 
200 
201 template<class Scalar>
202 void Thyra::update( Scalar alpha, const MultiVectorBase<Scalar>& U,
203  const Ptr<MultiVectorBase<Scalar> > &V )
204 {
205  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::null;
206  RTOpPack::TOpAXPY<Scalar> axpy_op(alpha);
207  applyOp<Scalar>( axpy_op, tuple(ptrInArg(U)), tuple(V), null );
208 }
209 
210 
211 template<class Scalar>
212 void Thyra::update( const ArrayView<const Scalar> &alpha, Scalar beta,
213  const MultiVectorBase<Scalar>& U, const Ptr<MultiVectorBase<Scalar> > &V )
214 {
215 #ifdef TEUCHOS_DEBUG
216  bool is_compatible = U.range()->isCompatible(*V->range());
217  TEUCHOS_TEST_FOR_EXCEPTION(
218  !is_compatible, Exceptions::IncompatibleVectorSpaces,
219  "update(...), Error, U.range()->isCompatible((V->range())==false");
220  is_compatible = U.domain()->isCompatible(*V->domain());
221  TEUCHOS_TEST_FOR_EXCEPTION(
222  !is_compatible, Exceptions::IncompatibleVectorSpaces,
223  "update(...), Error, U.domain().isCompatible(V->domain())==false");
224 #endif
225  const int m = U.domain()->dim();
226  for( int j = 0; j < m; ++j )
227  Vp_StV<Scalar>( V->col(j).ptr(), alpha[j]*beta, *U.col(j) );
228 }
229 
230 
231 template<class Scalar>
232 void Thyra::update( const MultiVectorBase<Scalar>& U,
233  const ArrayView<const Scalar> &alpha, Scalar beta,
234  const Ptr<MultiVectorBase<Scalar> > &V )
235 {
236 #ifdef TEUCHOS_DEBUG
237  bool is_compatible = U.range()->isCompatible(*V->range());
238  TEUCHOS_TEST_FOR_EXCEPTION(
239  !is_compatible, Exceptions::IncompatibleVectorSpaces,
240  "update(...), Error, U.range()->isCompatible((V->range())==false");
241  is_compatible = U.domain()->isCompatible(*V->domain());
242  TEUCHOS_TEST_FOR_EXCEPTION(
243  !is_compatible, Exceptions::IncompatibleVectorSpaces,
244  "update(...), Error, U.domain().isCompatible(V->domain())==false");
245 #endif
246  const int m = U.domain()->dim();
247  for( int j = 0; j < m; ++j ) {
248  Vt_S<Scalar>( V->col(j).ptr(), alpha[j]*beta );
249  Vp_StV<Scalar>( V->col(j).ptr(), 1.0, *U.col(j) );
250  }
251 }
252 
253 
254 template<class Scalar>
255 void Thyra::linear_combination(
256  const ArrayView<const Scalar> &alpha,
257  const ArrayView<const Ptr<const MultiVectorBase<Scalar> > > &X,
258  const Scalar &beta,
259  const Ptr<MultiVectorBase<Scalar> > &Y
260  )
261 {
262  using Teuchos::tuple; using Teuchos::null;
263 #ifdef TEUCHOS_DEBUG
264  TEUCHOS_ASSERT_EQUALITY(alpha.size(),X.size());
265 #endif
266  const int m = alpha.size();
267  if ( beta == ScalarTraits<Scalar>::one() && m == 1 ) {
268  update( alpha[0], *X[0], Y );
269  return;
270  }
271  else if (m == 0) {
272  scale( beta, Y );
273  return;
274  }
275  RTOpPack::TOpLinearCombination<Scalar> lin_comb_op(alpha, beta);
276  Thyra::applyOp<Scalar>( lin_comb_op, X, tuple(Y), null );
277 }
278 
279 
280 template<class Scalar>
281 void Thyra::randomize( Scalar l, Scalar u,
282  const Ptr<MultiVectorBase<Scalar> > &V )
283 {
284  const int m = V->domain()->dim();
285  for( int j = 0; j < m; ++j )
286  randomize( l, u, V->col(j).ptr() );
287  // Todo: call applyOp(...) directly!
288 }
289 
290 
291 template<class Scalar>
292 void Thyra::Vt_S( const Ptr<MultiVectorBase<Scalar> > &Z,
293  const Scalar& alpha )
294 {
295  const int m = Z->domain()->dim();
296  for( int j = 0; j < m; ++j )
297  Vt_S( Z->col(j).ptr(), alpha );
298  // Todo: call applyOp(...) directly!
299 }
300 
301 
302 template<class Scalar>
303 void Thyra::Vp_S( const Ptr<MultiVectorBase<Scalar> > &Z,
304  const Scalar& alpha )
305 {
306  const int m = Z->domain()->dim();
307  for( int j = 0; j < m; ++j )
308  Vp_S( Z->col(j).ptr(), alpha );
309  // Todo: call applyOp(...) directly!
310 }
311 
312 
313 template<class Scalar>
314 void Thyra::Vp_V( const Ptr<MultiVectorBase<Scalar> > &Z,
315  const MultiVectorBase<Scalar>& X )
316 {
317  using Teuchos::tuple; using Teuchos::ptrInArg;
318  typedef Teuchos::ScalarTraits<Scalar> ST;
319  linear_combination<Scalar>( tuple(ST::one()), tuple(ptrInArg(X)),
320  ST::one(), Z );
321 }
322 
323 
324 template<class Scalar>
325 void Thyra::V_VpV( const Ptr<MultiVectorBase<Scalar> > &Z,
326  const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
327 {
328  using Teuchos::tuple; using Teuchos::ptrInArg;
329  typedef Teuchos::ScalarTraits<Scalar> ST;
330  linear_combination<Scalar>(
331  tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
332  ST::zero(), Z
333  );
334 }
335 
336 
337 template<class Scalar>
338 void Thyra::V_VmV( const Ptr<MultiVectorBase<Scalar> > &Z,
339  const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y )
340 {
341  using Teuchos::tuple; using Teuchos::ptrInArg; using Teuchos::as;
342  typedef Teuchos::ScalarTraits<Scalar> ST;
343  linear_combination<Scalar>(
344  tuple(ST::one(), as<Scalar>(-ST::one())), tuple(ptrInArg(X), ptrInArg(Y)),
345  ST::zero(), Z
346  );
347 }
348 
349 
350 template<class Scalar>
351 void Thyra::V_StVpV(
352  const Ptr<MultiVectorBase<Scalar> > &Z, const Scalar &alpha,
353  const MultiVectorBase<Scalar>& X, const MultiVectorBase<Scalar>& Y
354  )
355 {
356  using Teuchos::tuple; using Teuchos::ptrInArg;
357  typedef Teuchos::ScalarTraits<Scalar> ST;
358  linear_combination<Scalar>(
359  tuple(alpha, ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
360  ST::zero(), Z
361  );
362 }
363 
364 
365 //
366 // Explicit instant macro
367 //
368 
369 #define THYRA_MULTI_VECTOR_STD_OPS_INSTANT(SCALAR) \
370  \
371  template void norms( const MultiVectorBase<SCALAR >& V, \
372  const ArrayView<ScalarTraits<SCALAR >::magnitudeType> &norms ); \
373  \
374  template void dots( const MultiVectorBase<SCALAR >& V1, const MultiVectorBase<SCALAR >& V2, \
375  const ArrayView<SCALAR > &dots ); \
376  \
377  template void sums( const MultiVectorBase<SCALAR >& V, const ArrayView<SCALAR > &sums ); \
378  \
379  template Teuchos::ScalarTraits<SCALAR >::magnitudeType \
380  norm_1( const MultiVectorBase<SCALAR >& V ); \
381  \
382  template void scale( SCALAR alpha, const Ptr<MultiVectorBase<SCALAR > > &V ); \
383  \
384  template void scaleUpdate( const VectorBase<SCALAR >& a, \
385  const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
386  \
387  template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, SCALAR alpha ); \
388  \
389  template void assign( const Ptr<MultiVectorBase<SCALAR > > &V, \
390  const MultiVectorBase<SCALAR >& U ); \
391  \
392  template void update( SCALAR alpha, const MultiVectorBase<SCALAR >& U, \
393  const Ptr<MultiVectorBase<SCALAR > > &V ); \
394  \
395  template void update( const ArrayView<const SCALAR > &alpha, SCALAR beta, \
396  const MultiVectorBase<SCALAR >& U, const Ptr<MultiVectorBase<SCALAR > > &V ); \
397  \
398  template void update( const MultiVectorBase<SCALAR >& U, \
399  const ArrayView<const SCALAR > &alpha, SCALAR beta, \
400  const Ptr<MultiVectorBase<SCALAR > > &V ); \
401  \
402  template void linear_combination( \
403  const ArrayView<const SCALAR > &alpha, \
404  const ArrayView<const Ptr<const MultiVectorBase<SCALAR > > > &X, \
405  const SCALAR &beta, \
406  const Ptr<MultiVectorBase<SCALAR > > &Y \
407  ); \
408  \
409  template void randomize( SCALAR l, SCALAR u, \
410  const Ptr<MultiVectorBase<SCALAR > > &V ); \
411  \
412  template void Vt_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
413  const SCALAR & alpha ); \
414  \
415  template void Vp_S( const Ptr<MultiVectorBase<SCALAR > > &Z, \
416  const SCALAR & alpha ); \
417  \
418  template void Vp_V( const Ptr<MultiVectorBase<SCALAR > > &Z, \
419  const MultiVectorBase<SCALAR >& X ); \
420  \
421  template void V_VpV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
422  const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
423  \
424  template void V_VmV( const Ptr<MultiVectorBase<SCALAR > > &Z, \
425  const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y ); \
426  \
427  template void V_StVpV( \
428  const Ptr<MultiVectorBase<SCALAR > > &Z, const SCALAR &alpha, \
429  const MultiVectorBase<SCALAR >& X, const MultiVectorBase<SCALAR >& Y \
430  ); \
431 
432 
433 #endif // THYRA_MULTI_VECTOR_STD_OPS_HPP
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(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 LinearOpBase< Scalar > > scale(const Scalar &scalar, const RCP< const LinearOpBase< Scalar > > &Op, const std::string &label="")
Build an implicit const scaled linear operator.
void norms(const MultiVectorBase< Scalar > &V, const ArrayView< typename ScalarTraits< Scalar >::magnitudeType > &norms)
Column-wise multi-vector natural norm.
void randomize(Scalar l, Scalar u, const Ptr< MultiVectorBase< Scalar > > &V)
Generate a random multi-vector with elements uniformly distributed elements.
void dots(const MultiVectorBase< Scalar > &V1, const MultiVectorBase< Scalar > &V2, const ArrayView< Scalar > &dots)
Multi-vector dot product.
void update(Scalar alpha, const MultiVectorBase< Scalar > &U, const Ptr< MultiVectorBase< Scalar > > &V)
alpha*U + V -> V.
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.