CppAD: A C++ Algorithmic Differentiation Package  20130918
color_general.hpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines