CppAD: A C++ Algorithmic Differentiation Package  20171217
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
sparse_internal.hpp
Go to the documentation of this file.
1 # ifndef CPPAD_LOCAL_SPARSE_INTERNAL_HPP
2 # define CPPAD_LOCAL_SPARSE_INTERNAL_HPP
3 
4 /* --------------------------------------------------------------------------
5 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
6 
7 CppAD is distributed under multiple licenses. This distribution is under
8 the terms of the
9  Eclipse Public License Version 1.0.
10 
11 A copy of this license is included in the COPYING file of this distribution.
12 Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
13 -------------------------------------------------------------------------- */
14 
15 // necessary definitions
16 # include <cppad/core/define.hpp>
20 
21 namespace CppAD { namespace local { // BEGIN_CPPAD_LOCAL_NAMESPACE
22 /*!
23 \file sparse_internal.hpp
24 Routines that enable code to be independent of which internal spasity pattern
25 is used.
26 */
27 // ---------------------------------------------------------------------------
28 /*!
29 Template structure used obtain the internal sparsity pattern type
30 form the corresponding element type.
31 The general form is not valid, must use a specialization.
32 
33 \tparam Element_type
34 type of an element in the sparsity structrue.
35 
36 \par <code>internal_sparsity<Element_type>::pattern_type</code>
37 is the type of the corresponding internal sparsity pattern.
38 */
39 template <class Element_type> struct internal_sparsity;
40 /// Specilization for \c bool elements.
41 template <>
42 struct internal_sparsity<bool>
43 {
45 };
46 /// Specilization for <code>std::set<size_t></code> elements.
47 template <>
48 struct internal_sparsity< std::set<size_t> >
49 {
51 };
52 // ---------------------------------------------------------------------------
53 /*!
54 Update the internal sparsity pattern for a sub-set of rows
55 
56 \tparam SizeVector
57 The type used for index sparsity patterns. This is a simple vector
58 with elements of type size_t.
59 
60 \tparam InternalSparsitiy
61 The type used for intenal sparsity patterns. This can be either
62 sparse_pack or sparse_list.
63 
64 \param zero_empty
65 If this is true, the internal sparstity pattern corresponds to row zero
66 must be empty on input and will be emtpy output; i.e., any corresponding
67 values in pattern_in will be ignored.
68 
69 \param input_empty
70 If this is true, the initial sparsity pattern for row
71 internal_index[i] is empty for all i.
72 In this case, one is setting the sparsity patterns; i.e.,
73 the output pattern in row internal_index[i] is the corresponding
74 entries in pattern.
75 
76 \param transpose
77 If this is true, pattern_in is transposed.
78 
79 \param internal_index
80 This specifies the sub-set of rows in internal_sparsity that we are updating.
81 If traspose is false (true),
82 this is the mapping from row (column) index in pattern_in to the corresponding
83 row index in the internal_pattern.
84 
85 \param internal_pattern
86 On input, the number of sets internal_pattern.n_set(),
87 and possible elements internal_pattern.end(), have been set.
88 If input_empty is true, and all of the sets
89 in internal_index are empty on input.
90 On output, the entries in pattern_in are added to internal_pattern.
91 To be specific, suppose transpose is false, and (i, j) is a possibly
92 non-zero entry in pattern_in, the entry (internal_index[i], j) is added
93 to internal_pattern.
94 On the other hand, if transpose is true,
95 the entry (internal_index[j], i) is added to internal_pattern.
96 
97 \param pattern_in
98 This is the sparsity pattern for variables,
99 or its transpose, depending on the value of transpose.
100 */
101 template <class SizeVector, class InternalSparsity>
103  bool zero_empty ,
104  bool input_empty ,
105  bool transpose ,
106  const vector<size_t>& internal_index ,
107  InternalSparsity& internal_pattern ,
108  const sparse_rc<SizeVector>& pattern_in )
109 {
110  size_t nr = internal_index.size();
111 # ifndef NDEBUG
112  size_t nc = internal_pattern.end();
113  if( transpose )
114  { CPPAD_ASSERT_UNKNOWN( pattern_in.nr() == nc );
115  CPPAD_ASSERT_UNKNOWN( pattern_in.nc() == nr );
116  }
117  else
118  { CPPAD_ASSERT_UNKNOWN( pattern_in.nr() == nr );
119  CPPAD_ASSERT_UNKNOWN( pattern_in.nc() == nc );
120  }
121  if( input_empty ) for(size_t i = 0; i < nr; i++)
122  { size_t i_var = internal_index[i];
123  CPPAD_ASSERT_UNKNOWN( internal_pattern.number_elements(i_var) == 0 );
124  }
125 # endif
126  const SizeVector& row( pattern_in.row() );
127  const SizeVector& col( pattern_in.col() );
128  size_t nnz = row.size();
129  for(size_t k = 0; k < nnz; k++)
130  { size_t r = row[k];
131  size_t c = col[k];
132  if( transpose )
133  std::swap(r, c);
134  //
135  size_t i_var = internal_index[r];
136  CPPAD_ASSERT_UNKNOWN( i_var < internal_pattern.n_set() );
137  CPPAD_ASSERT_UNKNOWN( c < nc );
138  bool ignore = zero_empty && i_var == 0;
139  if( ! ignore )
140  internal_pattern.post_element( internal_index[r], c );
141  }
142  // process posts
143  for(size_t i = 0; i < nr; ++i)
144  internal_pattern.process_post( internal_index[i] );
145 }
146 template <class InternalSparsity>
148  bool zero_empty ,
149  bool input_empty ,
150  bool transpose ,
151  const vector<size_t>& internal_index ,
152  InternalSparsity& internal_pattern ,
153  const vectorBool& pattern_in )
154 { size_t nr = internal_index.size();
155  size_t nc = internal_pattern.end();
156 # ifndef NDEBUG
157  CPPAD_ASSERT_UNKNOWN( pattern_in.size() == nr * nc );
158  if( input_empty ) for(size_t i = 0; i < nr; i++)
159  { size_t i_var = internal_index[i];
160  CPPAD_ASSERT_UNKNOWN( internal_pattern.number_elements(i_var) == 0 );
161  }
162 # endif
163  for(size_t i = 0; i < nr; i++)
164  { for(size_t j = 0; j < nc; j++)
165  { bool flag = pattern_in[i * nc + j];
166  if( transpose )
167  flag = pattern_in[j * nr + i];
168  if( flag )
169  { size_t i_var = internal_index[i];
170  CPPAD_ASSERT_UNKNOWN( i_var < internal_pattern.n_set() );
171  CPPAD_ASSERT_UNKNOWN( j < nc );
172  bool ignore = zero_empty && i_var == 0;
173  if( ! ignore )
174  internal_pattern.post_element( i_var, j);
175  }
176  }
177  }
178  // process posts
179  for(size_t i = 0; i < nr; ++i)
180  internal_pattern.process_post( internal_index[i] );
181  return;
182 }
183 template <class InternalSparsity>
185  bool zero_empty ,
186  bool input_empty ,
187  bool transpose ,
188  const vector<size_t>& internal_index ,
189  InternalSparsity& internal_pattern ,
190  const vector<bool>& pattern_in )
191 { size_t nr = internal_index.size();
192  size_t nc = internal_pattern.end();
193 # ifndef NDEBUG
194  CPPAD_ASSERT_UNKNOWN( pattern_in.size() == nr * nc );
195  if( input_empty ) for(size_t i = 0; i < nr; i++)
196  { size_t i_var = internal_index[i];
197  CPPAD_ASSERT_UNKNOWN( internal_pattern.number_elements(i_var) == 0 );
198  }
199 # endif
200  for(size_t i = 0; i < nr; i++)
201  { for(size_t j = 0; j < nc; j++)
202  { bool flag = pattern_in[i * nc + j];
203  if( transpose )
204  flag = pattern_in[j * nr + i];
205  if( flag )
206  { size_t i_var = internal_index[i];
207  CPPAD_ASSERT_UNKNOWN( i_var < internal_pattern.n_set() );
208  CPPAD_ASSERT_UNKNOWN( j < nc );
209  bool ignore = zero_empty && i_var == 0;
210  if( ! ignore )
211  internal_pattern.post_element( i_var, j);
212  }
213  }
214  }
215  // process posts
216  for(size_t i = 0; i < nr; ++i)
217  internal_pattern.process_post( internal_index[i] );
218  return;
219 }
220 template <class InternalSparsity>
222  bool zero_empty ,
223  bool input_empty ,
224  bool transpose ,
225  const vector<size_t>& internal_index ,
226  InternalSparsity& internal_pattern ,
227  const vector< std::set<size_t> >& pattern_in )
228 { size_t nr = internal_index.size();
229  size_t nc = internal_pattern.end();
230 # ifndef NDEBUG
231  if( input_empty ) for(size_t i = 0; i < nr; i++)
232  { size_t i_var = internal_index[i];
233  CPPAD_ASSERT_UNKNOWN( internal_pattern.number_elements(i_var) == 0 );
234  }
235 # endif
236  if( transpose )
237  { CPPAD_ASSERT_UNKNOWN( pattern_in.size() == nc );
238  for(size_t j = 0; j < nc; j++)
239  { std::set<size_t>::const_iterator itr( pattern_in[j].begin() );
240  while( itr != pattern_in[j].end() )
241  { size_t i = *itr;
242  size_t i_var = internal_index[i];
243  CPPAD_ASSERT_UNKNOWN( i_var < internal_pattern.n_set() );
244  CPPAD_ASSERT_UNKNOWN( j < nc );
245  bool ignore = zero_empty && i_var == 0;
246  if( ! ignore )
247  internal_pattern.post_element( i_var, j);
248  ++itr;
249  }
250  }
251  }
252  else
253  { CPPAD_ASSERT_UNKNOWN( pattern_in.size() == nr );
254  for(size_t i = 0; i < nr; i++)
255  { std::set<size_t>::const_iterator itr( pattern_in[i].begin() );
256  while( itr != pattern_in[i].end() )
257  { size_t j = *itr;
258  size_t i_var = internal_index[i];
259  CPPAD_ASSERT_UNKNOWN( i_var < internal_pattern.n_set() );
260  CPPAD_ASSERT_UNKNOWN( j < nc );
261  bool ignore = zero_empty && i_var == 0;
262  if( ! ignore )
263  internal_pattern.post_element( i_var, j);
264  ++itr;
265  }
266  }
267  }
268  // process posts
269  for(size_t i = 0; i < nr; ++i)
270  internal_pattern.process_post( internal_index[i] );
271  return;
272 }
273 // ---------------------------------------------------------------------------
274 /*!
275 Get sparsity pattern for a sub-set of variables
276 
277 \tparam SizeVector
278 The type used for index sparsity patterns. This is a simple vector
279 with elements of type size_t.
280 
281 \tparam InternalSparsitiy
282 The type used for intenal sparsity patterns. This can be either
283 sparse_pack or sparse_list.
284 
285 \param transpose
286 If this is true, pattern_out is transposed.
287 
288 \param internal_index
289 If transpose is false (true)
290 this is the mapping from row (column) an index in pattern_out
291 to the corresponding row index in internal_pattern.
292 
293 \param internal_pattern
294 This is the internal sparsity pattern.
295 
296 \param pattern_out
297 The input value of pattern_out does not matter.
298 Upon return it is an index sparsity pattern for each of the variables
299 in internal_index, or its transpose, depending on the value of transpose.
300 */
301 template <class SizeVector, class InternalSparsity>
303  bool transpose ,
304  const vector<size_t>& internal_index ,
305  const InternalSparsity& internal_pattern ,
306  sparse_rc<SizeVector>& pattern_out )
307 { typedef typename InternalSparsity::const_iterator iterator;
308  // number variables
309  size_t nr = internal_index.size();
310  // column size of interanl sparstiy pattern
311  size_t nc = internal_pattern.end();
312  // determine nnz, the number of possibly non-zero index pairs
313  size_t nnz = 0;
314  for(size_t i = 0; i < nr; i++)
315  { CPPAD_ASSERT_UNKNOWN( internal_index[i] < internal_pattern.n_set() );
316  iterator itr(internal_pattern, internal_index[i]);
317  size_t j = *itr;
318  while( j < nc )
319  { ++nnz;
320  j = *(++itr);
321  }
322  }
323  // transposed
324  if( transpose )
325  { pattern_out.resize(nc, nr, nnz);
326  //
327  size_t k = 0;
328  for(size_t i = 0; i < nr; i++)
329  { iterator itr(internal_pattern, internal_index[i]);
330  size_t j = *itr;
331  while( j < nc )
332  { pattern_out.set(k++, j, i);
333  j = *(++itr);
334  }
335  }
336  return;
337  }
338  // not transposed
339  pattern_out.resize(nr, nc, nnz);
340  //
341  size_t k = 0;
342  for(size_t i = 0; i < nr; i++)
343  { iterator itr(internal_pattern, internal_index[i]);
344  size_t j = *itr;
345  while( j < nc )
346  { pattern_out.set(k++, i, j);
347  j = *(++itr);
348  }
349  }
350  return;
351 }
352 template <class InternalSparsity>
354  bool transpose ,
355  const vector<size_t>& internal_index ,
356  const InternalSparsity& internal_pattern ,
357  vectorBool& pattern_out )
358 { typedef typename InternalSparsity::const_iterator iterator;
359  // number variables
360  size_t nr = internal_index.size();
361  //
362  // column size of interanl sparstiy pattern
363  size_t nc = internal_pattern.end();
364  //
365  pattern_out.resize(nr * nc);
366  for(size_t ij = 0; ij < nr * nc; ij++)
367  pattern_out[ij] = false;
368  //
369  for(size_t i = 0; i < nr; i++)
370  { CPPAD_ASSERT_UNKNOWN( internal_index[i] < internal_pattern.n_set() );
371  iterator itr(internal_pattern, internal_index[i]);
372  size_t j = *itr;
373  while( j < nc )
374  { if( transpose )
375  pattern_out[j * nr + i] = true;
376  else
377  pattern_out[i * nc + j] = true;
378  j = *(++itr);
379  }
380  }
381  return;
382 }
383 template <class InternalSparsity>
385  bool transpose ,
386  const vector<size_t>& internal_index ,
387  const InternalSparsity& internal_pattern ,
388  vector<bool>& pattern_out )
389 { typedef typename InternalSparsity::const_iterator iterator;
390  // number variables
391  size_t nr = internal_index.size();
392  //
393  // column size of interanl sparstiy pattern
394  size_t nc = internal_pattern.end();
395  //
396  pattern_out.resize(nr * nc);
397  for(size_t ij = 0; ij < nr * nc; ij++)
398  pattern_out[ij] = false;
399  //
400  for(size_t i = 0; i < nr; i++)
401  { CPPAD_ASSERT_UNKNOWN( internal_index[i] < internal_pattern.n_set() );
402  iterator itr(internal_pattern, internal_index[i]);
403  size_t j = *itr;
404  while( j < nc )
405  { if( transpose )
406  pattern_out[j * nr + i] = true;
407  else
408  pattern_out[i * nc + j] = true;
409  j = *(++itr);
410  }
411  }
412  return;
413 }
414 template <class InternalSparsity>
416  bool transpose ,
417  const vector<size_t>& internal_index ,
418  const InternalSparsity& internal_pattern ,
419  vector< std::set<size_t> >& pattern_out )
420 { typedef typename InternalSparsity::const_iterator iterator;
421  // number variables
422  size_t nr = internal_index.size();
423  //
424  // column size of interanl sparstiy pattern
425  size_t nc = internal_pattern.end();
426  //
427  if( transpose )
428  pattern_out.resize(nc);
429  else
430  pattern_out.resize(nr);
431  for(size_t k = 0; k < pattern_out.size(); k++)
432  pattern_out[k].clear();
433  //
434  for(size_t i = 0; i < nr; i++)
435  { CPPAD_ASSERT_UNKNOWN( internal_index[i] < internal_pattern.n_set() );
436  iterator itr(internal_pattern, internal_index[i]);
437  size_t j = *itr;
438  while( j < nc )
439  { if( transpose )
440  pattern_out[j].insert(i);
441  else
442  pattern_out[i].insert(j);
443  j = *(++itr);
444  }
445  }
446  return;
447 }
448 
449 
450 
451 } } // END_CPPAD_LOCAL_NAMESPACE
452 
453 # endif
Vector of sets of positive integers, each set stored as a singly linked list.
Definition: sparse_list.hpp:35
const SizeVector & col(void) const
column indices
Definition: sparse_rc.hpp:296
Define processor symbols and macros that are used by CppAD.
Vector of sets of positive integers stored as singly linked lists with the element values strictly in...
void get_internal_sparsity(bool transpose, const vector< size_t > &internal_index, const InternalSparsity &internal_pattern, sparse_rc< SizeVector > &pattern_out)
Get sparsity pattern for a sub-set of variables.
Vector of sets of positive integers stored as a packed array of bools.
void resize(size_t n)
change the number of elements in this vector.
Definition: vector.hpp:399
size_t nr(void) const
number of rows in matrix
Definition: sparse_rc.hpp:284
size_t size(void) const
number of elements currently in this vector.
Definition: vector.hpp:387
Vector of sets of postivie integers, each set stored as a packed boolean array.
Definition: sparse_pack.hpp:33
sparsity pattern for a matrix with indices of type size_t
Definition: sparse_rc.hpp:212
void resize(size_t nr, size_t nc, size_t nnz)
resize
Definition: sparse_rc.hpp:258
#define CPPAD_ASSERT_UNKNOWN(exp)
Check that exp is true, if not terminate execution.
size_t nc(void) const
number of columns in matrix
Definition: sparse_rc.hpp:287
void set_internal_sparsity(bool zero_empty, bool input_empty, bool transpose, const vector< size_t > &internal_index, InternalSparsity &internal_pattern, const sparse_rc< SizeVector > &pattern_in)
Update the internal sparsity pattern for a sub-set of rows.
void set(size_t k, size_t r, size_t c)
set row and column for a possibly non-zero element
Definition: sparse_rc.hpp:266
void resize(size_t n)
change number of elements in this vector
Definition: vector.hpp:721
Vector of sets of positive integers stored as size_t vector with the element values strictly increasi...
const SizeVector & row(void) const
row indices
Definition: sparse_rc.hpp:293
size_t size(void) const
number of elements in this vector
Definition: vector.hpp:713
Template structure used obtain the internal sparsity pattern type form the corresponding element type...