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 //
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 
59 template<class Scalar>
60 void Thyra::norms( const MultiVectorBase<Scalar>& V,
61  const ArrayView<typename ScalarTraits<Scalar>::magnitudeType> &norms )
62 {
63  typedef Teuchos::ScalarTraits<Scalar> ST;
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 
72 template<class Scalar>
73 void Thyra::dots( const MultiVectorBase<Scalar>& V1, const MultiVectorBase<Scalar>& V2,
74  const ArrayView<Scalar> &dots )
75 {
76  V2.dots(V1, dots);
77 }
78 
79 
80 template<class Scalar>
81 void 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 
100 template<class Scalar>
101 typename Teuchos::ScalarTraits<Scalar>::magnitudeType
102 Thyra::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 
121 template<class Scalar>
122 void Thyra::scale( Scalar alpha, const Ptr<MultiVectorBase<Scalar> > &V )
123 {
124  V->scale(alpha);
125 }
126 
127 
128 template<class Scalar>
129 void 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());
134  TEUCHOS_TEST_FOR_EXCEPTION(
135  !is_compatible, Exceptions::IncompatibleVectorSpaces,
136  "update(...), Error, U.range()->isCompatible(*a.space())==false" );
137  is_compatible = U.range()->isCompatible(*V->range());
138  TEUCHOS_TEST_FOR_EXCEPTION(
139  !is_compatible, Exceptions::IncompatibleVectorSpaces,
140  "update(...), Error, U.range()->isCompatible((V->range())==false" );
141  is_compatible = U.domain()->isCompatible(*V->domain());
142  TEUCHOS_TEST_FOR_EXCEPTION(
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 
153 template<class Scalar>
154 void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V, Scalar alpha )
155 {
156  V->assign(alpha);
157 }
158 
159 
160 template<class Scalar>
161 void Thyra::assign( const Ptr<MultiVectorBase<Scalar> > &V,
162  const MultiVectorBase<Scalar>& U )
163 {
164  V->assign(U);
165 }
166 
167 
168 template<class Scalar>
169 void Thyra::update( Scalar alpha, const MultiVectorBase<Scalar>& U,
170  const Ptr<MultiVectorBase<Scalar> > &V )
171 {
172  V->update(alpha, U);
173 }
174 
175 
176 template<class Scalar>
177 void 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());
182  TEUCHOS_TEST_FOR_EXCEPTION(
183  !is_compatible, Exceptions::IncompatibleVectorSpaces,
184  "update(...), Error, U.range()->isCompatible((V->range())==false");
185  is_compatible = U.domain()->isCompatible(*V->domain());
186  TEUCHOS_TEST_FOR_EXCEPTION(
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 
196 template<class Scalar>
197 void 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());
203  TEUCHOS_TEST_FOR_EXCEPTION(
204  !is_compatible, Exceptions::IncompatibleVectorSpaces,
205  "update(...), Error, U.range()->isCompatible((V->range())==false");
206  is_compatible = U.domain()->isCompatible(*V->domain());
207  TEUCHOS_TEST_FOR_EXCEPTION(
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 
219 template<class Scalar>
220 void 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 
231 template<class Scalar>
232 void 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 
242 template<class Scalar>
243 void Thyra::Vt_S( const Ptr<MultiVectorBase<Scalar> > &Z,
244  const Scalar& alpha )
245 {
246  Z->scale(alpha);
247 }
248 
249 
250 template<class Scalar>
251 void 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 
261 template<class Scalar>
262 void Thyra::Vp_V( const Ptr<MultiVectorBase<Scalar> > &Z,
263  const MultiVectorBase<Scalar>& X )
264 {
265  using Teuchos::tuple; using Teuchos::ptrInArg;
266  typedef Teuchos::ScalarTraits<Scalar> ST;
267  linear_combination<Scalar>( tuple(ST::one()), tuple(ptrInArg(X)),
268  ST::one(), Z );
269 }
270 
271 
272 template<class Scalar>
273 void 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;
277  typedef Teuchos::ScalarTraits<Scalar> ST;
278  linear_combination<Scalar>(
279  tuple(ST::one(), ST::one()), tuple(ptrInArg(X), ptrInArg(Y)),
280  ST::zero(), Z
281  );
282 }
283 
284 
285 template<class Scalar>
286 void 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;
290  typedef Teuchos::ScalarTraits<Scalar> ST;
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 
298 template<class Scalar>
299 void 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;
305  typedef Teuchos::ScalarTraits<Scalar> ST;
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
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 by   doxygen 1.8.5