Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_Sparse3TensorPartition.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) 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 Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef STOKHOS_SPARSE_3_TENSOR_PARTITION_HPP
43#define STOKHOS_SPARSE_3_TENSOR_PARTITION_HPP
44
46
47#include "Teuchos_ArrayRCP.hpp"
48#include "Teuchos_Array.hpp"
49
50namespace Stokhos {
51
52 template <typename TupleType>
53 class RCB {
54 public:
55 typedef typename TupleType::value_type coord_type;
56 typedef typename Teuchos::ArrayView<TupleType>::size_type size_type;
57 typedef typename TupleType::id_type id_type;
58
59 struct CoordCompare {
61 CoordCompare(const size_type& d_) : d(d_) {}
62 bool operator() (const TupleType& a, const TupleType& b) const {
63 return a(d) < b(d);
64 }
65 };
66
67 struct Box {
71 Teuchos::RCP<Box> left, right;
72 Teuchos::Array<TupleType> coords;
73
74 Box(const Teuchos::ArrayView<TupleType>& c) : coords(c.begin(), c.end()) {
77 }
78
79 // Compute bounding box around points
81 xmin = coords[0](0); xmax = coords[0](0);
82 ymin = coords[0](1); ymax = coords[0](1);
83 zmin = coords[0](2); zmax = coords[0](2);
84 for (size_type i=0; i<coords.size(); ++i) {
85 coord_type x = coords[i](0);
86 coord_type y = coords[i](1);
87 coord_type z = coords[i](2);
88 if (x < xmin) xmin = x;
89 if (y < ymin) ymin = y;
90 if (z < zmin) zmin = z;
91 if (x > xmax) xmax = x;
92 if (y > ymax) ymax = y;
93 if (z > zmax) zmax = z;
94 }
95 delta_x = xmax - xmin + 1;
96 delta_y = ymax - ymin + 1;
97 delta_z = zmax - zmin + 1;
98
99 // std::cout << "delta_x = " << delta_x
100 // << " delta_y = " << delta_y
101 // << " delta_z = " << delta_z << std::endl;
102 }
103
104 // Compute dimension to split
106 split_dim = 0;
107 if (delta_y >= delta_x && delta_y >= delta_z) split_dim = 1;
108 if (delta_z >= delta_x && delta_z >= delta_y) split_dim = 2;
109
110 //std::cout << "splitting dimension = " << split_dim << std::endl;
111 }
112
113 // Split box into two pieces with roughly equal numbers of points
114 void split() {
115 // Sort points based on splitting dimension
117 std::sort(coords.begin(), coords.end(), cmp);
118
119 // Divide coords into two bins of roughly equal size, keeping
120 // coords with equal values for split dimension together
121 size_type n = coords.size();
122 size_type s = n / 2;
123
124 while (s < n-1 && coords[s](split_dim) == coords[s+1](split_dim)) ++s;
125 //std::cout << "n = " << n << " s = " << s << std::endl;
126
127 if (s > 0)
128 left = Teuchos::rcp(new Box(coords.view(0, s)));
129 if (s < n)
130 right = Teuchos::rcp(new Box(coords.view(s, n-s)));
131
132 // Clear my coordinate array since we aren't a leaf
133 //Teuchos::Array<TupleType>().swap(coords);
134 //coords.resize(0);
135 }
136
137 };
138
140 RCB(const coord_type& max_length_,
141 const size_type& max_parts_,
142 const Teuchos::ArrayView<TupleType>& coords_) :
143 max_length(max_length_),
144 max_parts(max_parts_),
145 coords(coords_.begin(), coords_.end()) {
146 partition();
147 }
148
150 ~RCB() {}
151
154
156 Teuchos::RCP<Box> get_partition_root() const { return root; }
157
159 Teuchos::RCP< Teuchos::Array< Teuchos::RCP<Box> > > get_parts() const {
160 return parts; }
161
162 // Create list of part IDs for each tuple
163 Teuchos::ArrayRCP<id_type> get_part_IDs() const {
164 Teuchos::ArrayRCP<id_type> part_ids(coords.size());
165 for (size_type part=0; part<num_parts; ++part) {
166 Teuchos::RCP<Box> box = (*parts)[part];
167 size_type n = box->coords.size();
168 for (size_type i=0; i<n; ++i)
169 part_ids[ box->coords[i].ID() ] = part;
170 }
171 return part_ids;
172 }
173
174 private:
175
178 Teuchos::Array<TupleType> coords;
179 Teuchos::RCP<Box> root;
180 Teuchos::RCP< Teuchos::Array< Teuchos::RCP<Box> > > parts;
181
183 void partition() {
184
185 // Create root bounding box
186 root = Teuchos::rcp(new Box(coords()));
187 num_parts = 1;
188 parts = Teuchos::rcp(new Teuchos::Array< Teuchos::RCP<Box> >);
189
190 // Create array of boxes that are too big
191 Teuchos::Array< Teuchos::RCP<Box> > boxes_to_split;
192 if (root->delta_x > max_length ||
193 root->delta_y > max_length ||
194 root->delta_z > max_length)
195 boxes_to_split.push_back(root);
196 else
197 parts->push_back(root);
198
199 // Split each box until all boxes are less than tolerance
200 while(boxes_to_split.size() > 0 && num_parts < max_parts) {
201 Teuchos::RCP<Box> box = boxes_to_split.back();
202 boxes_to_split.pop_back();
203 box->split();
204 ++num_parts;
205 if (box->left != Teuchos::null) {
206 if (box->left->delta_x > max_length ||
207 box->left->delta_y > max_length ||
208 box->left->delta_z > max_length)
209 boxes_to_split.push_back(box->left);
210 else
211 parts->push_back(box->left);
212 }
213 if (box->right != Teuchos::null) {
214 if (box->right->delta_x > max_length ||
215 box->right->delta_y > max_length ||
216 box->right->delta_z > max_length)
217 boxes_to_split.push_back(box->right);
218 else
219 parts->push_back(box->right);
220 }
221 }
222
223 TEUCHOS_ASSERT(parts->size() == num_parts);
224 }
225
226 };
227
228 template <typename ordinal_type, typename scalar_type>
229 struct CijkData {
230 typedef ordinal_type value_type;
231 typedef ordinal_type id_type;
232 ordinal_type gid;
233 ordinal_type i, j, k;
234 scalar_type c;
235
236 ordinal_type operator() (ordinal_type d) const {
237 if (d == 0) return i;
238 if (d == 1) return j;
239 if (d == 2) return k;
240 return -1;
241 }
242
243 ordinal_type ID() const { return gid; }
244 };
245
250 };
251
252 template <typename ordinal_type, typename scalar_type>
253 Teuchos::ArrayRCP< CijkData<ordinal_type,scalar_type> >
256 CijkSymmetryType symmetry_type) {
258 typedef typename Cijk_type::k_iterator k_iterator;
259 typedef typename Cijk_type::kj_iterator kj_iterator;
260 typedef typename Cijk_type::kji_iterator kji_iterator;
261
262 ordinal_type num_cijk = Cijk.num_entries();
263 Teuchos::ArrayRCP< CijkData<ordinal_type,scalar_type> > coordinate_list(
264 num_cijk);
265 ordinal_type idx = 0;
266 k_iterator k_begin = Cijk.k_begin();
267 k_iterator k_end = Cijk.k_end();
268 for (k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
269 ordinal_type k = index(k_it);
270 kj_iterator j_begin = Cijk.j_begin(k_it);
271 kj_iterator j_end = Cijk.j_end(k_it);
272 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
273 ordinal_type j = index(j_it);
274 kji_iterator i_begin = Cijk.i_begin(j_it);
275 kji_iterator i_end = Cijk.i_end(j_it);
276 for (kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
277 ordinal_type i = index(i_it);
278 if (symmetry_type == CIJK_NO_SYMMETRY) {
279 coordinate_list[idx].i = i;
280 coordinate_list[idx].j = j;
281 coordinate_list[idx].k = k;
282 coordinate_list[idx].c = value(i_it);
283 coordinate_list[idx].gid = idx;
284 ++idx;
285 }
286 else if (symmetry_type == CIJK_TWO_WAY_SYMMETRY && j >= k) {
287 coordinate_list[idx].i = i;
288 coordinate_list[idx].j = j;
289 coordinate_list[idx].k = k;
290 if (j == k)
291 coordinate_list[idx].c = 0.5*value(i_it);
292 else
293 coordinate_list[idx].c = value(i_it);
294 coordinate_list[idx].gid = idx;
295 ++idx;
296 }
297 else if (symmetry_type == CIJK_SIX_WAY_SYMMETRY && i >= j && j >= k) {
298 coordinate_list[idx].i = i;
299 coordinate_list[idx].j = j;
300 coordinate_list[idx].k = k;
301 if (i == j && j == k)
302 coordinate_list[idx].c = (1.0/6.0)*value(i_it);
303 else
304 coordinate_list[idx].c = value(i_it);
305 coordinate_list[idx].gid = idx;
306 ++idx;
307 }
308 }
309 }
310 }
311 coordinate_list.resize(idx);
312
313 return coordinate_list;
314 }
315
316} // namespace Stokhos
317
318#endif // STOKHOS_SPARSE_3_TENSOR_PARTITION_HPP
Teuchos::Array< TupleType > coords
Teuchos::ArrayRCP< id_type > get_part_IDs() const
size_type get_num_parts() const
Get number of parts.
Teuchos::RCP< Box > get_partition_root() const
Get root of partition.
Teuchos::RCP< Teuchos::Array< Teuchos::RCP< Box > > > get_parts() const
Get parts array.
TupleType::value_type coord_type
Teuchos::ArrayView< TupleType >::size_type size_type
RCB(const coord_type &max_length_, const size_type &max_parts_, const Teuchos::ArrayView< TupleType > &coords_)
Constructor.
Teuchos::RCP< Teuchos::Array< Teuchos::RCP< Box > > > parts
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
k_iterator k_begin() const
Iterator pointing to first k entry.
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
ordinal_type num_entries() const
Return number of non-zero entries.
k_iterator k_end() const
Iterator pointing to last k entry.
Top-level namespace for Stokhos classes and functions.
Teuchos::ArrayRCP< CijkData< ordinal_type, scalar_type > > build_cijk_coordinate_list(const Sparse3Tensor< ordinal_type, scalar_type > &Cijk, CijkSymmetryType symmetry_type)
ordinal_type operator()(ordinal_type d) const
Teuchos::Array< TupleType > coords
Box(const Teuchos::ArrayView< TupleType > &c)
bool operator()(const TupleType &a, const TupleType &b) const