Collection of Concrete Vector Reduction/Transformation Operator Implementations  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
RTOpPack_TOpLinearCombination_def.hpp
1 // @HEADER
2 // *****************************************************************************
3 // RTOp: Interfaces and Support Software for Vector Reduction Transformation
4 // Operations
5 //
6 // Copyright 2006 NTESS and the RTOp contributors.
7 // SPDX-License-Identifier: BSD-3-Clause
8 // *****************************************************************************
9 // @HEADER
10 
11 #ifndef RTOPPACK_TOP_LINEAR_COMBINATION_DEF_HPP
12 #define RTOPPACK_TOP_LINEAR_COMBINATION_DEF_HPP
13 
14 
15 #include "Teuchos_Workspace.hpp"
16 
17 
18 namespace RTOpPack {
19 
20 
21 template<class Scalar>
22 TOpLinearCombination<Scalar>::TOpLinearCombination(
23  const ArrayView<const Scalar> &alpha_in,
24  const Scalar &beta_in
25  )
26  :beta_(beta_in)
27 {
28  if (alpha_in.size())
29  this->alpha(alpha_in);
30  this->setOpNameBase("TOpLinearCombination");
31 }
32 
33 
34 
35 template<class Scalar>
36 void TOpLinearCombination<Scalar>::alpha(
37  const ArrayView<const Scalar> &alpha_in )
38 {
39  TEUCHOS_TEST_FOR_EXCEPT( alpha_in.size() == 0 );
40  alpha_ = alpha_in;
41 }
42 
43 
44 template<class Scalar>
45 const ArrayView<const Scalar>
46 TOpLinearCombination<Scalar>::alpha() const
47 { return alpha_; }
48 
49 
50 template<class Scalar>
51 void TOpLinearCombination<Scalar>::beta( const Scalar& beta_in ) { beta_ = beta_in; }
52 
53 
54 template<class Scalar>
55 Scalar TOpLinearCombination<Scalar>::beta() const { return beta_; }
56 
57 
58 template<class Scalar>
59 int TOpLinearCombination<Scalar>::num_vecs() const { return alpha_.size(); }
60 
61 
62 // Overridden from RTOpT
63 
64 
65 template<class Scalar>
66 void TOpLinearCombination<Scalar>::apply_op_impl(
67  const ArrayView<const ConstSubVectorView<Scalar> > &sub_vecs,
68  const ArrayView<const SubVectorView<Scalar> > &targ_sub_vecs,
69  const Ptr<ReductTarget> &reduct_obj_inout
70  ) const
71 {
72 
73  using Teuchos::as;
74  using Teuchos::Workspace;
75  typedef Teuchos::ScalarTraits<Scalar> ST;
76  typedef typename Teuchos::ArrayRCP<Scalar>::iterator iter_t;
77  typedef typename Teuchos::ArrayRCP<const Scalar>::iterator const_iter_t;
78  Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
79 
80 #ifdef TEUCHOS_DEBUG
81  validate_apply_op<Scalar>(*this, as<int>(alpha_.size()), 1, false,
82  sub_vecs, targ_sub_vecs, reduct_obj_inout.getConst());
83 #else
84  (void)reduct_obj_inout;
85 #endif
86 
87  const int l_num_vecs = alpha_.size();
88 
89  // Get iterators to local data
90  const RTOpPack::index_type subDim = targ_sub_vecs[0].subDim();
91  iter_t z0_val = targ_sub_vecs[0].values().begin();
92  const ptrdiff_t z0_s = targ_sub_vecs[0].stride();
93  Workspace<const_iter_t> v_val(wss,l_num_vecs);
94  Workspace<ptrdiff_t> v_s(wss,l_num_vecs,false);
95  for( int k = 0; k < l_num_vecs; ++k ) {
96 #ifdef TEUCHOS_DEBUG
97  TEUCHOS_TEST_FOR_EXCEPT( sub_vecs[k].subDim() != subDim );
98  TEUCHOS_TEST_FOR_EXCEPT( sub_vecs[k].globalOffset() != targ_sub_vecs[0].globalOffset() );
99 #endif
100  v_val[k] = sub_vecs[k].values().begin();
101  v_s[k] = sub_vecs[k].stride();
102  }
103 
104  //
105  // Perform the operation and specialize the cases for l_num_vecs = 1 and 2
106  // in order to get good performance.
107  //
108  if( l_num_vecs == 1 ) {
109  //
110  // z0 = alpha*v0 + beta*z0
111  //
112  const Scalar l_alpha = alpha_[0], l_beta = beta_;
113  const_iter_t v0_val = v_val[0];
114  const ptrdiff_t v0_s = v_s[0];
115  if( l_beta==ST::zero() ) {
116  // z0 = alpha*v0
117  if( z0_s==1 && v0_s==1 ) {
118  for( int j = 0; j < subDim; ++j )
119  (*z0_val++) = l_alpha * (*v0_val++);
120  }
121  else {
122  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s )
123  (*z0_val) = l_alpha * (*v0_val);
124  }
125  }
126  else if( l_beta==ST::one() ) {
127  //
128  // z0 = alpha*v0 + z0
129  //
130  if( z0_s==1 && v0_s==1 ) {
131  for( int j = 0; j < subDim; ++j )
132  (*z0_val++) += l_alpha * (*v0_val++);
133  }
134  else {
135  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s )
136  (*z0_val) += l_alpha * (*v0_val);
137  }
138  }
139  else {
140  // z0 = alpha*v0 + beta*z0
141  if( z0_s==1 && v0_s==1 ) {
142  for( int j = 0; j < subDim; ++j, ++z0_val )
143  (*z0_val) = l_alpha * (*v0_val++) + l_beta*(*z0_val);
144  }
145  else {
146  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s )
147  (*z0_val) = l_alpha * (*v0_val) + l_beta*(*z0_val);
148  }
149  }
150  }
151  else if( l_num_vecs == 2 ) {
152  //
153  // z0 = alpha0*v0 + alpha1*v1 + beta*z0
154  //
155  const Scalar alpha0 = alpha_[0], alpha1=alpha_[1], l_beta = beta_;
156  const_iter_t v0_val = v_val[0];
157  const ptrdiff_t v0_s = v_s[0];
158  const_iter_t v1_val = v_val[1];
159  const ptrdiff_t v1_s = v_s[1];
160  if( l_beta==ST::zero() ) {
161  if( alpha0 == ST::one() ) {
162  if( alpha1 == ST::one() ) {
163  // z0 = v0 + v1
164  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
165  for( int j = 0; j < subDim; ++j )
166  (*z0_val++) = (*v0_val++) + (*v1_val++);
167  }
168  else {
169  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
170  (*z0_val) = (*v0_val) + (*v1_val);
171  }
172  }
173  else {
174  // z0 = v0 + alpha1*v1
175  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
176  for( int j = 0; j < subDim; ++j )
177  (*z0_val++) = (*v0_val++) + alpha1*(*v1_val++);
178  }
179  else {
180  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
181  (*z0_val) = (*v0_val) + alpha1*(*v1_val);
182  }
183  }
184  }
185  else {
186  if( alpha1 == ST::one() ) {
187  // z0 = alpha0*v0 + v1
188  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
189  for( int j = 0; j < subDim; ++j )
190  (*z0_val++) = alpha0*(*v0_val++) + (*v1_val++);
191  }
192  else {
193  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
194  (*z0_val) = alpha0*(*v0_val) + (*v1_val);
195  }
196  }
197  else {
198  // z0 = alpha0*v0 + alpha1*v1
199  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
200  for( int j = 0; j < subDim; ++j )
201  (*z0_val++) = alpha0*(*v0_val++) + alpha1*(*v1_val++);
202  }
203  else {
204  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
205  (*z0_val) = alpha0*(*v0_val) + alpha1*(*v1_val);
206  }
207  }
208  }
209  }
210  else if( l_beta==ST::one() ) {
211  if( alpha0 == ST::one() ) {
212  if( alpha1 == ST::one() ) {
213  // z0 = v0 + v1 + z0
214  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
215  for( int j = 0; j < subDim; ++j, ++z0_val )
216  (*z0_val) += (*v0_val++) + (*v1_val++);
217  }
218  else {
219  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
220  (*z0_val) += (*v0_val) + (*v1_val);
221  }
222  }
223  else {
224  // z0 = v0 + alpha1*v1 + z0
225  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
226  for( int j = 0; j < subDim; ++j, ++z0_val )
227  (*z0_val) += (*v0_val++) + alpha1*(*v1_val++);
228  }
229  else {
230  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
231  (*z0_val) += (*v0_val) + alpha1*(*v1_val);
232  }
233  }
234  }
235  else {
236  if( alpha1 == ST::one() ) {
237  // z0 = alpha0*v0 + v1 + z0
238  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
239  for( int j = 0; j < subDim; ++j, ++z0_val )
240  (*z0_val) += alpha0*(*v0_val++) + (*v1_val++);
241  }
242  else {
243  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
244  (*z0_val) += alpha0*(*v0_val) + (*v1_val);
245  }
246  }
247  else {
248  // z0 = alpha0*v0 + alpha1*v1 + z0
249  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
250  for( int j = 0; j < subDim; ++j, ++z0_val )
251  (*z0_val) += alpha0*(*v0_val++) + alpha1*(*v1_val++);
252  }
253  else {
254  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
255  (*z0_val) += alpha0*(*v0_val) + alpha1*(*v1_val);
256  }
257  }
258  }
259  }
260  else {
261  if( alpha0 == ST::one() ) {
262  if( alpha1 == ST::one() ) {
263  // z0 = v0 + v1 + beta*z0
264  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
265  for( int j = 0; j < subDim; ++j, ++z0_val )
266  (*z0_val) = (*v0_val++) + (*v1_val++) + l_beta*(*z0_val);
267  }
268  else {
269  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
270  (*z0_val) = (*v0_val) + (*v1_val) + l_beta*(*z0_val);
271  }
272  }
273  else {
274  // z0 = v0 + alpha1*v1 + beta*z0
275  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
276  for( int j = 0; j < subDim; ++j, ++z0_val )
277  (*z0_val) = (*v0_val++) + alpha1*(*v1_val++) + l_beta*(*z0_val);
278  }
279  else {
280  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
281  (*z0_val) = (*v0_val) + alpha1*(*v1_val) + l_beta*(*z0_val);
282  }
283  }
284  }
285  else {
286  if( alpha1 == ST::one() ) {
287  // z0 = alpha0*v0 + v1 + beta*z0
288  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
289  for( int j = 0; j < subDim; ++j, ++z0_val )
290  (*z0_val) = alpha0*(*v0_val++) + (*v1_val++) + l_beta*(*z0_val);
291  }
292  else {
293  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
294  (*z0_val) = alpha0*(*v0_val) + (*v1_val) + l_beta*(*z0_val);
295  }
296  }
297  else {
298  // z0 = alpha0*v0 + alpha1*v1 + beta*z0
299  if( z0_s==1 && v0_s==1 && v1_s==1 ) {
300  for( int j = 0; j < subDim; ++j, ++z0_val )
301  (*z0_val) = alpha0*(*v0_val++) + alpha1*(*v1_val++) + l_beta*(*z0_val);
302  }
303  else {
304  for( int j = 0; j < subDim; ++j, z0_val+=z0_s, v0_val+=v0_s, v1_val+=v1_s )
305  (*z0_val) = alpha0*(*v0_val) + alpha1*(*v1_val) + l_beta*(*z0_val);
306  }
307  }
308  }
309  }
310  }
311  else {
312  //
313  // Totally general implementation (but least efficient)
314  //
315  // z0 *= beta
316  if( beta_ == ST::zero() ) {
317  for( int j = 0; j < subDim; ++j, z0_val += z0_s )
318  (*z0_val) = ST::zero();
319  }
320  else if( beta_ != ST::one() ) {
321  for( int j = 0; j < subDim; ++j, z0_val += z0_s )
322  (*z0_val) *= beta_;
323  }
324  // z0 += sum( alpha[k]*v[k], k=0...l_num_vecs-1)
325  z0_val = targ_sub_vecs[0].values().begin();
326  for( int j = 0; j < subDim; ++j, z0_val += z0_s ) {
327  for( int k = 0; k < l_num_vecs; ++k ) {
328  const Scalar
329  &alpha_k = alpha_[k],
330  &v_k_val = *v_val[k];
331  (*z0_val) += alpha_k * v_k_val;
332  v_val[k] += v_s[k];
333  }
334  }
335  }
336 }
337 
338 
339 } // namespace RTOpPack
340 
341 
342 #endif // RTOPPACK_TOP_LINEAR_COMBINATION_DEF_HPP
RTOpPack::TOpLinearCombination::apply_op_impl
void apply_op_impl(const ArrayView< const ConstSubVectorView< Scalar > > &sub_vecs, const ArrayView< const SubVectorView< Scalar > > &targ_sub_vecs, const Ptr< ReductTarget > &reduct_obj_inout) const
Definition: RTOpPack_TOpLinearCombination_def.hpp:66
RTOpPack::TOpLinearCombination::alpha
const ArrayView< const Scalar > alpha() const
Definition: RTOpPack_TOpLinearCombination_def.hpp:46
RTOpPack::RTOpT< Scalar >::setOpNameBase
void setOpNameBase(const std::string &op_name_base)
Teuchos::as
TypeTo as(const TypeFrom &t)
RTOpPack::TOpLinearCombination::TOpLinearCombination
TOpLinearCombination(const ArrayView< const Scalar > &alpha_in=Teuchos::null, const Scalar &beta=Teuchos::ScalarTraits< Scalar >::zero())
Definition: RTOpPack_TOpLinearCombination_def.hpp:22
TEUCHOS_TEST_FOR_EXCEPT
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
Teuchos::get_default_workspace_store
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()
Generated by   doxygen 1.8.5