AbstractLinAlgPack: C++ Interfaces For Vectors, Matrices And Related Linear Algebra Objects  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AbstractLinAlgPack_COOMatrixPartitionedViewClassDef.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) 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 (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef COO_MATRIX_PARTITIONED_VIEW_CLASS_DEF_H
43 #define COO_MATRIX_PARTITIONED_VIEW_CLASS_DEF_H
44 
45 #include <sstream>
46 #include <algorithm>
47 #include <functional>
48 
49 #include "AbstractLinAlgPack_COOMatrixPartitionedViewClassDecl.hpp"
50 
51 namespace AbstractLinAlgPack {
52 
53 // /////////////////////////////////////////////////////////////////////////////
54 // Template function definitions for COOMatrixPartitionedView
55 
56 template <class T_Indice, class T_Value>
59  , size_type cols
60  , size_type nz
61  , value_type val[]
62  , const indice_type ivect[]
63  , const indice_type jvect[]
64  , const size_type inv_row_perm[]
65  , const size_type inv_col_perm[]
66  , const size_type num_row_part
67  , const size_type row_part[]
68  , const size_type num_col_part
69  , const size_type col_part[]
70  , const EPartitionOrder partition_order )
71 {
72  try {
73 
74  // 1) Check some preconditins before we start
75 
76  // Check the sparsisty density of COO matrix
77  if(nz > rows * cols)
78  throw std::out_of_range(
79  "COOMatrixPartitionedView<...>::create_view() : Input error, "
80  "nz can not be greater than rows * cols");
81 
82  // Check row and column partition information
83 
84  if(num_row_part > rows || num_col_part > rows)
85  throw std::out_of_range(
86  "COOMatrixPartitionedView<...>::create_view() : Input error, "
87  "num_rows_part and num_col_part can not be greater than"
88  " rows and cols respectivly");
89 
90  if(row_part[num_row_part] > rows + 1)
91  throw std::out_of_range(
92  "COOMatrixPartitionedView<...>::create_view() : Input error, "
93  "row_part[num_row_part] can not be greater than rows");
94 
95  if(col_part[num_col_part] > cols + 1)
96  throw std::out_of_range(
97  "COOMatrixPartitionedView<...>::create_view() : Input error, "
98  "col_part[num_col_part] can not be greater than cols");
99 
100  {for(size_type i = 1; i < num_row_part + 1; ++i)
101  if(row_part[i-1] >= row_part[i])
102  throw std::domain_error(
103  "COOMatrixPartitionedView<...>::create_view() : Input error, "
104  "row_part[i-1] < row_part[i] does not hold");}
105 
106  {for(size_type i = 1; i < num_col_part + 1; ++i)
107  if(col_part[i-1] >= col_part[i])
108  throw std::domain_error(
109  "COOMatrixPartitionedView<...>::create_view() : Input error, "
110  "col_part[i-1] < col_part[i] does not hold");}
111 
112  // Get references to referenced quantities
113  std::vector<size_type>
114  &_row_part = ref_row_part_.obj(),
115  &_col_part = ref_col_part_.obj(),
116  &_part_start = ref_part_start_.obj();
117  ele_type
118  &_ele = ref_ele_.obj();
119 
120  // 2) Initialize storage for data members and copy data
121  num_row_part_ = num_row_part;
122  num_col_part_ = num_col_part;
123  _row_part.resize(num_row_part_ + 1);
124  std::copy(row_part, row_part+ num_row_part_+1, _row_part.begin());
125  _col_part.resize(num_col_part_ + 1);
126  std::copy(col_part, col_part+ num_col_part_+1, _col_part.begin());
127  partition_order_ = partition_order;
128  _ele.resize(nz,element_type()); // hack to get around compiler error.
129  _part_start.resize(num_row_part_ * num_col_part_+1);
130  _part_start.assign(_part_start.size(),0); // set to 0
131 
132  // 3) Count the number of nonzero elements in each overall partition
133  {
134  // Use the storage locations _part_start[1] ... _part_start[num_part]
135  // to count the number of nonzero elements in each overall partition.
136  // In particular num_in_part[overall_p - 1] will give the number
137  // of nonzero elements in the overall partition number overall_p.
138  //
139  // _part_start = { 0, nz1, nz2,...,nzn }
140  //
141  size_type *num_in_part = &_part_start[0] + 1;
142 
143  // Loop (time = O(nz), space = O(1))
144 
145  // read iterators
146  const indice_type *ivect_itr = ivect,
147  *ivect_itr_end = ivect + nz,
148  *jvect_itr = jvect;
149 
150  for(;ivect_itr != ivect_itr_end; ++ivect_itr, ++jvect_itr) {
151  // Get row and column indices in the non-permuted matrix
152  indice_type i_org = *ivect_itr,
153  j_org = *jvect_itr;
154  // assert that they are in range
155  if(i_org < 1 || i_org > rows || j_org < 1 || j_org > cols) {
156  std::ostringstream omsg;
157  omsg << "COOMatrixPartitionedView<...>::create_view() : "
158  " Error, element k = " << ivect_itr - ivect
159  << " in the non-permuted matrix"
160  " is out of range with rows = " << rows
161  << ", cols = " << cols << ", i = " << i_org
162  << ", and j = " << j_org;
163  throw std::out_of_range(omsg.str());
164  }
165  // Get row and column indices in the permuted matrix
166  indice_type i = inv_row_perm[i_org - 1],
167  j = inv_col_perm[j_org - 1];
168  // assert that they are in range
169  if(i < 1 || i > rows || j < 1 || j > cols) {
170  std::ostringstream omsg;
171  omsg << "COOMatrixPartitionedView<...>::create_view() : "
172  " Error, element k = " << ivect_itr - ivect
173  << " in the permuted matrix"
174  " is out of range with rows = " << rows
175  << ", cols = " << cols << ", i = " << i
176  << ", and j = " << j;
177  throw std::out_of_range(omsg.str());
178  }
179  // get the overall partition number
180  size_type overall_p = overall_p_from_ij(_row_part,_col_part,i,j);
181  // Increment the number of nonzero elements.
182  num_in_part[overall_p - 1]++;
183  }
184  }
185 
186  // 4) Set part_start[ovarall_p] equal to the start in ele
187  // for the nonzero elements in that partition.
188  //
189  // _part_start = { start_1 = 0, start_2,..., start_n, ###}
190  //
191  {for(size_type i = 2; i < num_row_part_ * num_col_part_; ++i)
192  _part_start[i] += _part_start[i-1];}
193 
194  // 5) Shift the elements over
195  //
196  // _part_start = { 0, start_1 = 0, start_2,..., start_n}
197  //
198  {for(size_type i = num_row_part_ * num_col_part_; i > 0; --i)
199  _part_start[i] = _part_start[i-1];}
200 
201  // 6) Add the nonzero elements to each partition. When we
202  // are done we should have _part_start initialized properly.
203  //
204  // part_start = { start_1 = 0, start_2,..., start_n, total_nz }
205  //
206  {
207  // next_ele_insert[overall_p - 1] is the possition in ele
208  // for the next element to ensert
209  size_type *next_ele_insert = &_part_start[0] + 1;
210 
211  // Loop (time = O(nz), space = O(1))
212 
213  // read iterators
214  value_type *val_itr = val,
215  *val_itr_end = val + nz;
216  const indice_type *ivect_itr = ivect,
217  *jvect_itr = jvect;
218 
219  for(;val_itr != val_itr_end; ++val_itr, ++ivect_itr, ++jvect_itr) {
220  // Get row and column indices in the permuted matrix
221  indice_type i = inv_row_perm[*ivect_itr - 1],
222  j = inv_col_perm[*jvect_itr - 1];
223  // get the overall partition number
224  size_type overall_p = overall_p_from_ij(_row_part,_col_part,i,j);
225  // Add the element to the partition
226  _ele[next_ele_insert[overall_p - 1]++].initialize(val_itr,i,j);
227  }
228  }
229 
230  } // end try
231  catch(...) {
232  free();
233  throw; // rethrow the exception out of here
234  }
235 }
236 
237 template <class T_Indice, class T_Value>
239 {
240  num_row_part_ = coo_view.num_row_part_;
241  num_col_part_ = coo_view.num_col_part_;
242  ref_row_part_ = coo_view.ref_row_part_;
243  ref_col_part_ = coo_view.ref_col_part_;
244  partition_order_ = coo_view.partition_order_;
245  ref_ele_ = coo_view.ref_ele_;
246  ref_part_start_ = coo_view.ref_part_start_;
247 }
248 
249 template <class T_Indice, class T_Value>
251 {
252  // Reinitialize to uninitizlied
253  num_row_part_ = num_col_part_ = 0;
254  if(ref_row_part_.has_ref_set()) ref_row_part_.obj().resize(0);
255  if(ref_col_part_.has_ref_set()) ref_col_part_.obj().resize(0);
256  if(ref_ele_.has_ref_set()) ref_ele_.obj().resize(0,element_type());
257  if(ref_part_start_.has_ref_set()) ref_part_start_.obj().resize(0);
258 }
259 
260 template <class T_Indice, class T_Value>
262 {
263  assert_initialized();
264  const std::vector<size_type> &_row_part = ref_row_part_.const_obj();
265  std::copy(_row_part.begin(), _row_part.end(), row_part);
266 }
267 
268 template <class T_Indice, class T_Value>
270 {
271  assert_initialized();
272  const std::vector<size_type> &_col_part = ref_col_part_.const_obj();
273  std::copy(_col_part.begin(), _col_part.end(), col_part);
274 }
275 
276 template <class T_Indice, class T_Value>
278 COOMatrixPartitionedView<T_Indice,T_Value>::part_num(const vector_size_type& part
279  , size_type indice)
280 {
281  return std::upper_bound(part.begin(),part.end(),indice) - part.begin();
282 }
283 
284 template <class T_Indice, class T_Value>
285 COOMatrixPartitionedView<T_Indice,T_Value>::partition_type
286 COOMatrixPartitionedView<T_Indice,T_Value>::create_partition(Range1D rng_overall_p) const
287 {
288  assert_initialized();
289  // get reference to data structures
290  const std::vector<size_type>
291  &row_part = ref_row_part_.const_obj(),
292  &col_part = ref_col_part_.const_obj(),
293  &part_start = ref_part_start_.const_obj();
294  ele_type
295  &_ele = const_cast<ref_ele_type&>(ref_ele_).obj(); // This is ok
296  // Get upper and lower overall, row and column partition numbers
297  rng_overall_p = DenseLinAlgPack::full_range(rng_overall_p,1,num_row_part_*num_col_part_);
298  size_type l_p = rng_overall_p.lbound(),
299  u_p = rng_overall_p.ubound(),
300  l_r_p = imp_row_part_num(l_p),
301  l_c_p = imp_col_part_num(l_p),
302  u_r_p = imp_row_part_num(u_p),
303  u_c_p = imp_col_part_num(u_p);
304  // build argument list for creation of the partition.
305  size_type
306  rows = row_part[u_r_p] - row_part[l_r_p - 1],
307  cols = col_part[u_c_p] - col_part[l_c_p - 1],
308  nz = part_start[u_p] - part_start[l_p - 1];
309  element_type
310  *ele = &_ele[0] + part_start[l_p - 1];
311  difference_type
312  row_offset = - (row_part[l_r_p - 1] - 1),
313  col_offset = - (col_part[l_c_p - 1] - 1);
314 
315  return partition_type(rows,cols,nz,ele,row_offset,col_offset);
316 }
317 
318 } // end namespace AbstractLinAlgPack
319 
320 #endif // COO_MATRIX_PARTITIONED_VIEW_CLASS_DEF_H
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
void bind(const COOMatrixPartitionedView &coom_view)
Bind the view of another partitioned matrix.
void get_row_part(indice_type row_part[]) const
get the array of row partitions.
size_t size_type
Sparse pointer element type for a COO matrix (val, ivect, jvect).
void create_view(size_type rows, size_type cols, size_type nz, value_type val[], const indice_type ivect[], const indice_type jvect[], const size_type inv_row_perm[], const size_type inv_col_perm[], const size_type num_row_part, const size_type row_part[], const size_type num_col_part, const size_type col_part[], const EPartitionOrder partition_order)
Crete a view to a COO matrix.
void free()
Free the allocated memory and make uninitialized.
void get_col_part(indice_type col_part[]) const
get the array of column partitions.
size_type cols(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)