CppAD: A C++ Algorithmic Differentiation Package
20130918
|
00001 /* $Id$ */ 00002 # ifndef CPPAD_COLOR_GENERAL_INCLUDED 00003 # define CPPAD_COLOR_GENERAL_INCLUDED 00004 00005 /* -------------------------------------------------------------------------- 00006 CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-14 Bradley M. Bell 00007 00008 CppAD is distributed under multiple licenses. This distribution is under 00009 the terms of the 00010 Eclipse Public License Version 1.0. 00011 00012 A copy of this license is included in the COPYING file of this distribution. 00013 Please visit http://www.coin-or.org/CppAD/ for information on other licenses. 00014 -------------------------------------------------------------------------- */ 00015 # include <cppad/local/cppad_colpack.hpp> 00016 00017 namespace CppAD { // BEGIN_CPPAD_NAMESPACE 00018 /*! 00019 \file color_general.hpp 00020 Coloring algorithm for a general sparse matrix. 00021 */ 00022 // -------------------------------------------------------------------------- 00023 /*! 00024 Determine which rows of a general sparse matrix can be computed together; 00025 i.e., do not have non-zero overlapping values. 00026 00027 \tparam VectorSize 00028 is a simple vector class with elements of type \c size_t. 00029 00030 \tparam VectorSet 00031 is an unspecified type with the exception that it must support the 00032 operations under \c pattern and the following operations where 00033 \c p is a \c VectorSet object: 00034 \n 00035 <code>VectorSet p</code> 00036 Constructs a new vector of sets object. 00037 \n 00038 <code>p.resize(ns, ne)</code> 00039 resizes \p to \c ns sets with elements between zero \c ne. 00040 All of the \c ns sets are initially empty. 00041 \n 00042 <code>p.add_element(s, e)</code> 00043 add element \c e to set with index \c s. 00044 00045 \param pattern [in] 00046 Is a representation of the sparsity pattern for the matrix. 00047 Note that \c color_general does not change the values in \c pattern, 00048 but it is not \c const because its iterator facility modifies some of its 00049 internal data. 00050 \n 00051 <code>m = pattern.n_set()</code> 00052 \n 00053 sets \c m to the number of rows in the sparse matrix. 00054 All of the row indices are less than this value. 00055 \n 00056 <code>n = pattern.end()</code> 00057 \n 00058 sets \c n to the number of columns in the sparse matrix. 00059 All of the column indices are less than this value. 00060 \n 00061 <code>pattern.begin(i)</code> 00062 instructs the iterator facility to start iterating over 00063 columns in the i-th row of the sparsity pattern. 00064 \n 00065 <code>j = pattern.next_element()</code> 00066 Sets \c j to the next possibly non-zero column 00067 in the row specified by the previous call to <code>pattern.begin</code>. 00068 If there are no more such columns, the value 00069 <code>pattern.end()</code> is returned. 00070 00071 \param row [in] 00072 is a vector specifying which row indices to compute. 00073 00074 \param col [in] 00075 is a vector, with the same size as \c row, 00076 that specifies which column indices to compute. 00077 For each valid index \c k, the index pair 00078 <code>(row[k], col[k])</code> must be present in the sparsity pattern. 00079 It may be that some entries in the sparsity pattern do not need to be computed; 00080 i.e, do not appear in the set of 00081 <code>(row[k], col[k])</code> entries. 00082 00083 \param color [out] 00084 is a vector with size \c m. 00085 The input value of its elements does not matter. 00086 Upon return, it is a coloring for the rows of the sparse matrix. 00087 \n 00088 \n 00089 If for come \c i, <code>color[i] == m</code>, then 00090 the i-th row does not appear in the vector \c row. 00091 Otherwise, <code>color[i] < m</code>. 00092 \n 00093 \n 00094 Suppose two differen rows, <code>i != r</code> have the same color and 00095 column index \c j is such that both of the pairs 00096 <code>(i, j)</code> and <code>(r, j)</code> appear in the sparsity pattern. 00097 It follows that neither of these pairs appear in the set of 00098 <code>(row[k], col[k])</code> entries. 00099 \n 00100 \n 00101 This routine tries to minimize, with respect to the choice of colors, 00102 the maximum, with respct to \c k, of <code>color[ row[k] ]</code>. 00103 */ 00104 template <class VectorSet, class VectorSize> 00105 void color_general_cppad( 00106 VectorSet& pattern , 00107 const VectorSize& row , 00108 const VectorSize& col , 00109 CppAD::vector<size_t>& color ) 00110 { size_t i, j, k, ell, r; 00111 00112 size_t K = row.size(); 00113 size_t m = pattern.n_set(); 00114 size_t n = pattern.end(); 00115 00116 CPPAD_ASSERT_UNKNOWN( col.size() == K ); 00117 CPPAD_ASSERT_UNKNOWN( color.size() == m ); 00118 00119 // We define the set of rows, columns, and pairs that appear 00120 // by the set ( row[k], col[k] ) for k = 0, ... , K-1. 00121 00122 // initialize rows that appear 00123 CppAD::vector<bool> row_appear(m); 00124 for(i = 0; i < m; i++) 00125 row_appear[i] = false; 00126 00127 // rows and columns that appear 00128 VectorSet c2r_appear, r2c_appear; 00129 c2r_appear.resize(n, m); 00130 r2c_appear.resize(m, n); 00131 for(k = 0; k < K; k++) 00132 { CPPAD_ASSERT_UNKNOWN( pattern.is_element(row[k], col[k]) ); 00133 row_appear[ row[k] ] = true; 00134 c2r_appear.add_element(col[k], row[k]); 00135 r2c_appear.add_element(row[k], col[k]); 00136 } 00137 00138 // for each column, which rows are non-zero and do not appear 00139 VectorSet not_appear; 00140 not_appear.resize(n, m); 00141 for(i = 0; i < m; i++) 00142 { pattern.begin(i); 00143 j = pattern.next_element(); 00144 while( j != pattern.end() ) 00145 { if( ! c2r_appear.is_element(j , i) ) 00146 not_appear.add_element(j, i); 00147 j = pattern.next_element(); 00148 } 00149 } 00150 00151 // initial coloring 00152 color.resize(m); 00153 ell = 0; 00154 for(i = 0; i < m; i++) 00155 { if( row_appear[i] ) 00156 color[i] = ell++; 00157 else color[i] = m; 00158 } 00159 /* 00160 See GreedyPartialD2Coloring Algorithm Section 3.6.2 of 00161 Graph Coloring in Optimization Revisited by 00162 Assefaw Gebremedhin, Fredrik Maane, Alex Pothen 00163 00164 The algorithm above was modified (by Brad Bell) to take advantage of the 00165 fact that only the entries (subset of the sparsity pattern) specified by 00166 row and col need to be computed. 00167 */ 00168 CppAD::vector<bool> forbidden(m); 00169 for(i = 1; i < m; i++) // for each row that appears 00170 if( color[i] < m ) 00171 { 00172 // initial all colors as ok for this row 00173 // (value of forbidden for ell > initial color[i] does not matter) 00174 for(ell = 0; ell <= color[i]; ell++) 00175 forbidden[ell] = false; 00176 00177 // ----------------------------------------------------- 00178 // Forbid colors for which this row would destroy results: 00179 // 00180 // for each column that is non-zero for this row 00181 pattern.begin(i); 00182 j = pattern.next_element(); 00183 while( j != pattern.end() ) 00184 { // for each row that appears with this column 00185 c2r_appear.begin(j); 00186 r = c2r_appear.next_element(); 00187 while( r != c2r_appear.end() ) 00188 { // if this is not the same row, forbid its color 00189 if( (r < i) & (color[r] < m) ) 00190 forbidden[ color[r] ] = true; 00191 r = c2r_appear.next_element(); 00192 } 00193 j = pattern.next_element(); 00194 } 00195 00196 00197 // ----------------------------------------------------- 00198 // Forbid colors that destroy results needed for this row. 00199 // 00200 // for each column that appears with this row 00201 r2c_appear.begin(i); 00202 j = r2c_appear.next_element(); 00203 while( j != r2c_appear.end() ) 00204 { // For each row that is non-zero for this column 00205 // (the appear rows have already been checked above). 00206 not_appear.begin(j); 00207 r = not_appear.next_element(); 00208 while( r != not_appear.end() ) 00209 { // if this is not the same row, forbid its color 00210 if( (r < i) & (color[r] < m) ) 00211 forbidden[ color[r] ] = true; 00212 r = not_appear.next_element(); 00213 } 00214 j = r2c_appear.next_element(); 00215 } 00216 00217 // pick the color with smallest index 00218 ell = 0; 00219 while( forbidden[ell] ) 00220 { ell++; 00221 CPPAD_ASSERT_UNKNOWN( ell <= color[i] ); 00222 } 00223 color[i] = ell; 00224 } 00225 return; 00226 } 00227 00228 # if CPPAD_HAS_COLPACK 00229 /*! 00230 Colpack version of determining which rows of a sparse matrix 00231 can be computed together. 00232 00233 \copydetails color_general 00234 */ 00235 template <class VectorSet, class VectorSize> 00236 void color_general_colpack( 00237 VectorSet& pattern , 00238 const VectorSize& row , 00239 const VectorSize& col , 00240 CppAD::vector<size_t>& color ) 00241 { size_t i, j, k; 00242 size_t m = pattern.n_set(); 00243 size_t n = pattern.end(); 00244 00245 // Determine number of non-zero entries in each row 00246 CppAD::vector<size_t> n_nonzero(m); 00247 size_t n_nonzero_total = 0; 00248 for(i = 0; i < m; i++) 00249 { n_nonzero[i] = 0; 00250 pattern.begin(i); 00251 j = pattern.next_element(); 00252 while( j != pattern.end() ) 00253 { n_nonzero[i]++; 00254 j = pattern.next_element(); 00255 } 00256 n_nonzero_total += n_nonzero[i]; 00257 } 00258 00259 // Allocate memory and fill in Adolc sparsity pattern 00260 CppAD::vector<unsigned int*> adolc_pattern(m); 00261 CppAD::vector<unsigned int> adolc_memory(m + n_nonzero_total); 00262 size_t i_memory = 0; 00263 for(i = 0; i < m; i++) 00264 { adolc_pattern[i] = adolc_memory.data() + i_memory; 00265 adolc_pattern[i][0] = n_nonzero[i]; 00266 pattern.begin(i); 00267 j = pattern.next_element(); 00268 k = 1; 00269 while(j != pattern.end() ) 00270 { adolc_pattern[i][k++] = j; 00271 j = pattern.next_element(); 00272 } 00273 CPPAD_ASSERT_UNKNOWN( k == 1 + n_nonzero[i] ); 00274 i_memory += k; 00275 } 00276 CPPAD_ASSERT_UNKNOWN( i_memory == m + n_nonzero_total ); 00277 00278 // Must use an external routine for this part of the calculation because 00279 // ColPack/ColPackHeaders.h has as 'using namespace std' at global level. 00280 cppad_colpack_general(color, m, n, adolc_pattern); 00281 00282 return; 00283 } 00284 # endif // CPPAD_HAS_COLPACK 00285 00286 } // END_CPPAD_NAMESPACE 00287 # endif