RTOp Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RTOpPack_TOpSetSubVector_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // RTOp: Interfaces and Support Software for Vector Reduction Transformation
5 // Operations
6 // Copyright (2006) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39 //
40 // ***********************************************************************
41 // @HEADER
42 
43 #ifndef RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
44 #define RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
45 
46 
47 namespace RTOpPack {
48 
49 
50 template<class Scalar>
52  :RTOpT<Scalar>("TOpSetSubVector")
53 {}
54 
55 
56 template<class Scalar>
58  :RTOpT<Scalar>("TOpSetSubVector")
59 {
60  set_sub_vec(sub_vec);
61 }
62 
63 
64 template<class Scalar>
66 {
67  sub_vec_ = sub_vec;
68 }
69 
70 
71 // Overridden from RTOpT
72 
73 
74 template<class Scalar>
76 {
77  return false;
78 }
79 
80 
81 template<class Scalar>
83  const ArrayView<const ConstSubVectorView<Scalar> > &sub_vecs,
84  const ArrayView<const SubVectorView<Scalar> > &targ_sub_vecs,
85  const Ptr<ReductTarget> &reduct_obj
86  ) const
87 {
88 
89  typedef typename Teuchos::ArrayRCP<Scalar>::iterator iter_t;
90  typedef typename Teuchos::ArrayRCP<const Scalar>::iterator const_iter_t;
91  typedef typename Teuchos::ArrayRCP<const Teuchos_Ordinal>::iterator const_indices_iter_t;
92 
93  validate_apply_op( *this, 0, 1, false, sub_vecs, targ_sub_vecs, reduct_obj.getConst() );
94 
95  // Get the dense subvector data that we will set
96  const SubVectorView<Scalar> &z = targ_sub_vecs[0];
97  const index_type z_global_offset = z.globalOffset();
98  const index_type z_sub_dim = z.subDim();
99  iter_t z_val = z.values().begin();
100  const ptrdiff_t z_val_s = z.stride();
101 
102  // Get the sparse subvector data
103  const index_type v_global_offset = sub_vec_.globalOffset();
104  const index_type v_sub_dim = sub_vec_.subDim();
105  const index_type v_sub_nz = sub_vec_.subNz();
106  const_iter_t v_val = sub_vec_.values().begin();
107  const ptrdiff_t v_val_s = sub_vec_.valuesStride();
108  const bool has_v_ind = !is_null(sub_vec_.indices());
109  const_indices_iter_t v_ind = sub_vec_.indices().begin();
110  const ptrdiff_t v_ind_s = sub_vec_.indicesStride();
111  const ptrdiff_t v_l_off = sub_vec_.localOffset();
112  //const bool v_sorted = sub_vec_.isSorted();
113 
114  //
115  // Set part of the sub-vector v for this chunk z.
116  //
117 
118  if( v_global_offset + v_sub_dim < z_global_offset + 1
119  || z_global_offset + z_sub_dim < v_global_offset + 1 )
120  {
121  // The sub-vector that we are setting does not overlap with this vector
122  // chunk!
123  return;
124  }
125 
126  if( v_sub_nz == 0 )
127  return; // The sparse sub-vector we are reading from is empty?
128 
129  // Get the number of elements that overlap
130  index_type num_overlap;
131  if( v_global_offset <= z_global_offset ) {
132  if( v_global_offset + v_sub_dim >= z_global_offset + z_sub_dim )
133  num_overlap = z_sub_dim;
134  else
135  num_overlap = (v_global_offset + v_sub_dim) - z_global_offset;
136  }
137  else {
138  if( z_global_offset + z_sub_dim >= v_global_offset + v_sub_dim )
139  num_overlap = v_sub_dim;
140  else
141  num_overlap = (z_global_offset + z_sub_dim) - v_global_offset;
142  }
143 
144  // Set the part of the sub-vector that overlaps
145  if (has_v_ind) {
146  // Sparse elements
147  // Set the overlapping elements to zero first.
148  if( v_global_offset >= z_global_offset )
149  z_val += (v_global_offset - z_global_offset) * z_val_s;
150  for( index_type k = 0; k < num_overlap; ++k, z_val += z_val_s )
151  *z_val = 0.0;
152  // Now set the sparse entries
153  z_val = targ_sub_vecs[0].values().begin();
154  for( index_type k = 0; k < v_sub_nz; ++k, v_val += v_val_s, v_ind += v_ind_s ) {
155  const index_type i = v_global_offset + v_l_off + (*v_ind);
156  if( z_global_offset < i && i <= z_global_offset + z_sub_dim )
157  z_val[ z_val_s * (i - z_global_offset - 1) ] = *v_val;
158  }
159  // ToDo: Implement a faster version for v sorted and eliminate the
160  // if statement in the above loop.
161  }
162  else {
163  // Dense elements
164  if( v_global_offset <= z_global_offset )
165  v_val += (z_global_offset - v_global_offset) * v_val_s;
166  else
167  z_val += (v_global_offset - z_global_offset) * z_val_s;
168  for( index_type k = 0; k < num_overlap; ++k, v_val += v_val_s, z_val += z_val_s )
169  *z_val = *v_val;
170  }
171 }
172 
173 
174 } // namespace RTOpPack
175 
176 
177 #endif // RTOPPACK_TOP_SET_SUB_VECTOR_DEF_HPP
bool is_null(const boost::shared_ptr< T > &p)
void set_sub_vec(const SparseSubVectorT< Scalar > &sub_vec)
Class for a changeable sub-vector.
Class for a (sparse or dense) sub-vector.
Teuchos_Ordinal index_type
Class for a non-changeable sub-vector.
Templated interface to vector reduction/transformation operators {abstract}.
void validate_apply_op(const RTOpT< Scalar > &op, const int allowed_num_sub_vecs, const int allowed_num_targ_sub_vecs, const bool expect_reduct_obj, const ArrayView< const ConstSubVectorView< Scalar > > &sub_vecs, const ArrayView< const SubVectorView< Scalar > > &targ_sub_vecs, const Ptr< const ReductTarget > &reduct_obj)
Validate the input to an apply_op(...) function.
void apply_op_impl(const ArrayView< const ConstSubVectorView< Scalar > > &sub_vecs, const ArrayView< const SubVectorView< Scalar > > &targ_sub_vecs, const Ptr< ReductTarget > &reduct_obj) const
Ptr< const T > getConst() const