color_symmetric.hpp

Go to the documentation of this file.

# ifndef CPPAD_LOCAL_COLOR_SYMMETRIC_HPP
# define CPPAD_LOCAL_COLOR_SYMMETRIC_HPP

# include <cppad/configure.hpp>
# include <cppad/local/cppad_colpack.hpp>

/* --------------------------------------------------------------------------
CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell

CppAD is distributed under multiple licenses. This distribution is under
the terms of the
                    Eclipse Public License Version 1.0.

A copy of this license is included in the COPYING file of this distribution.
Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
-------------------------------------------------------------------------- */

namespace CppAD { namespace local { // BEGIN_CPPAD_LOCAL_NAMESPACE
/*!
\file color_symmetric.hpp
Coloring algorithm for a symmetric sparse matrix.
*/
// --------------------------------------------------------------------------
/*!
CppAD algorithm for determining which rows of a symmetric sparse matrix can be
computed together.

\tparam VectorSize
is a simple vector class with elements of type size_t.

\tparam VectorSet
is a vector_of_sets class.

\param pattern [in]
Is a representation of the sparsity pattern for the matrix.

\param row [in/out]
is a vector specifying which row indices to compute.

\param col [in/out]
is a vector, with the same size as row,
that specifies which column indices to compute.
\n
\n
Input:
For each  valid index \c k, the index pair
<code>(row[k], col[k])</code> must be present in the sparsity pattern.
It may be that some entries in the sparsity pattern do not need to be computed;
i.e, do not appear in the set of
<code>(row[k], col[k])</code> entries.
\n
\n
Output:
On output, some of row and column indices may have been swapped
\code
     std::swap( row[k], col[k] )
\endcode
So the the the color for row[k] can be used to compute entry
(row[k], col[k]).

\param color [out]
is a vector with size m.
The input value of its elements does not matter.
Upon return, it is a coloring for the rows of the sparse matrix.
Note that if color[i] == m, then there is no index k for which
row[k] == i (for the return value of row).
\n
\n
Fix any (i, j) in the sparsity pattern.
Suppose that there is a row index i1 with
i1 != i, color[i1] == color[i] and (i1, j) is in the sparsity pattern.
If follows that for all j1 with
j1 != j and color[j1] == color[j],
(j1, i ) is not in the sparsity pattern.
\n
\n
This routine tries to minimize, with respect to the choice of colors,
the maximum, with respect to k, of <code>color[ row[k] ]</code>.
*/
template <class VectorSet>
void color_symmetric_cppad(
     const VectorSet&        pattern   ,
     CppAD::vector<size_t>&  row       ,
     CppAD::vector<size_t>&  col       ,
     CppAD::vector<size_t>&  color     )
{
     size_t K = row.size();
     size_t m = pattern.n_set();
     CPPAD_ASSERT_UNKNOWN( m == pattern.end() );
     CPPAD_ASSERT_UNKNOWN( color.size() == m );
     CPPAD_ASSERT_UNKNOWN( col.size()   == K );

     // row, column pairs that appear in ( row[k], col[k] )
     CppAD::vector< std::set<size_t> > pair_needed(m);
     std::set<size_t>::iterator itr1, itr2;
     for(size_t k1 = 0;  k1 < K; k1++)
     {    CPPAD_ASSERT_UNKNOWN( pattern.is_element(row[k1], col[k1]) );
          pair_needed[ row[k1] ].insert( col[k1] );
          pair_needed[ col[k1] ].insert( row[k1] );
     }

     // order the rows decending by number of pairs needed
     CppAD::vector<size_t> key(m), order2row(m);
     for(size_t i1 = 0; i1 < m; i1++)
     {    CPPAD_ASSERT_UNKNOWN( pair_needed[i1].size() <= m );
          key[i1] = m - pair_needed[i1].size();
     }
     CppAD::index_sort(key, order2row);

     // mapping from order index to row index
     CppAD::vector<size_t> row2order(m);
     for(size_t o1 = 0; o1 < m; o1++)
          row2order[ order2row[o1] ] = o1;

     // initial coloring
     color.resize(m);
     size_t c1 = 0;
     for(size_t o1 = 0; o1 < m; o1++)
     {    size_t i1 = order2row[o1];
          if( pair_needed[i1].empty() )
               color[i1] = m;
          else
               color[i1] = c1++;
     }

     // which colors are forbidden for this row
     CppAD::vector<bool> forbidden(m);

     // must start with row zero so that we remove results computed for it
     for(size_t o1 = 0; o1 < m; o1++) // for each row that appears (in order)
     if( color[ order2row[o1] ] < m )
     {    size_t i1 = order2row[o1];
          c1 = color[i1];

          // initial all colors as ok for this row
          // (value of forbidden for c > c1 does not matter)
          for(size_t c2 = 0; c2 <= c1; c2++)
               forbidden[c2] = false;

          // -----------------------------------------------------
          // Forbid grouping with rows that would destroy results that are
          // needed for this row.
          itr1 = pair_needed[i1].begin();
          while( itr1 != pair_needed[i1].end() )
          {    // entry (i1, j1) is needed for this row
               size_t j1 = *itr1;

               // Forbid rows i2 != i1 that have non-zero sparsity at (i2, j1).
               // Note that this is the same as non-zero sparsity at (j1, i2)
               typename VectorSet::const_iterator pattern_itr(pattern, j1);
               size_t i2 = *pattern_itr;
               while( i2 != pattern.end() )
               {    size_t c2 = color[i2];
                    if( c2 < c1 )
                         forbidden[c2] = true;
                    i2 = *(++pattern_itr);
               }
               itr1++;
          }
          // -----------------------------------------------------
          // Forbid grouping with rows that this row would destroy results for
          for(size_t o2 = 0; o2 < o1; o2++)
          {    size_t i2 = order2row[o2];
               size_t c2 = color[i2];
               itr2 = pair_needed[i2].begin();
               while( itr2 != pair_needed[i2].end() )
               {    size_t j2 = *itr2;
                    // row i2 needs pair (i2, j2).
                    // Forbid grouping with i1 if (i1, j2) has non-zero sparsity
                    if( pattern.is_element(i1, j2) )
                         forbidden[c2] = true;
                    itr2++;
               }
          }

          // pick the color with smallest index
          size_t c2 = 0;
          while( forbidden[c2] )
          {    c2++;
               CPPAD_ASSERT_UNKNOWN( c2 <= c1 );
          }
          color[i1] = c2;

          // no longer need results that are computed by this row
          itr1 = pair_needed[i1].begin();
          while( itr1 != pair_needed[i1].end() )
          {    size_t j1 = *itr1;
               if( row2order[j1] > o1 )
               {    itr2 = pair_needed[j1].find(i1);
                    if( itr2 != pair_needed[j1].end() )
                    {    pair_needed[j1].erase(itr2);
                         if( pair_needed[j1].empty() )
                              color[j1] = m;
                    }
               }
               itr1++;
          }
     }

     // determine which sparsity entries need to be reflected
     for(size_t k1 = 0; k1 < row.size(); k1++)
     {    size_t i1   = row[k1];
          size_t j1   = col[k1];
          itr1 = pair_needed[i1].find(j1);
          if( itr1 == pair_needed[i1].end() )
          {    row[k1] = j1;
               col[k1] = i1;
# ifndef NDEBUG
               itr1 = pair_needed[j1].find(i1);
               CPPAD_ASSERT_UNKNOWN( itr1 != pair_needed[j1].end() );
# endif
          }
     }
     return;
}

// --------------------------------------------------------------------------
/*!
Colpack algorithm for determining which rows of a symmetric sparse matrix
can be computed together.

\copydetails CppAD::local::color_symmetric_cppad
*/
template <class VectorSet>
void color_symmetric_colpack(
     const VectorSet&        pattern   ,
     CppAD::vector<size_t>&  row       ,
     CppAD::vector<size_t>&  col       ,
     CppAD::vector<size_t>&  color     )
{
# if ! CPPAD_HAS_COLPACK
     CPPAD_ASSERT_UNKNOWN(false);
     return;
# else
     size_t i, j, k;
     size_t m = pattern.n_set();
     CPPAD_ASSERT_UNKNOWN( m == pattern.end() );
     CPPAD_ASSERT_UNKNOWN( row.size() == col.size() );

     // Determine number of non-zero entries in each row
     CppAD::vector<size_t> n_nonzero(m);
     size_t n_nonzero_total = 0;
     for(i = 0; i < m; i++)
     {    n_nonzero[i] = 0;
          typename VectorSet::const_iterator pattern_itr(pattern, i);
          j = *pattern_itr;
          while( j != pattern.end() )
          {    n_nonzero[i]++;
               j = *(++pattern_itr);
          }
          n_nonzero_total += n_nonzero[i];
     }

     // Allocate memory and fill in Adolc sparsity pattern
     CppAD::vector<unsigned int*> adolc_pattern(m);
     CppAD::vector<unsigned int>  adolc_memory(m + n_nonzero_total);
     size_t i_memory = 0;
     for(i = 0; i < m; i++)
     {    adolc_pattern[i]    = adolc_memory.data() + i_memory;
          CPPAD_ASSERT_KNOWN(
               std::numeric_limits<unsigned int>::max() >= n_nonzero[i],
               "Matrix is too large for colpack"
          );
          adolc_pattern[i][0] = static_cast<unsigned int>( n_nonzero[i] );
          typename VectorSet::const_iterator pattern_itr(pattern, i);
          j = *pattern_itr;
          k = 1;
          while(j != pattern.end() )
          {
               CPPAD_ASSERT_KNOWN(
                    std::numeric_limits<unsigned int>::max() >= j,
                    "Matrix is too large for colpack"
               );
               adolc_pattern[i][k++] = static_cast<unsigned int>( j );
               j = *(++pattern_itr);
          }
          CPPAD_ASSERT_UNKNOWN( k == 1 + n_nonzero[i] );
          i_memory += k;
     }
     CPPAD_ASSERT_UNKNOWN( i_memory == m + n_nonzero_total );

     // Must use an external routine for this part of the calculation because
     // ColPack/ColPackHeaders.h has as 'using namespace std' at global level.
     cppad_colpack_symmetric(color, m, adolc_pattern);

     // determine which sparsity entries need to be reflected
     for(size_t k1 = 0; k1 < row.size(); k1++)
     {    size_t i1 = row[k1];
          size_t j1 = col[k1];
          bool reflect = false;
          for(size_t i2 = 0; i2 < m; i2++)
          if( (i1 != i2) & (color[i1]==color[i2]) )
          {    for(size_t k2 = 1; k2 <= adolc_pattern[i2][0]; k2++)
               {    size_t j2 = adolc_pattern[i2][k2];
                    reflect |= (j1 == j2);
               }
          }
          if( reflect )
          {    row[k1] = j1;
               col[k1] = i1;
          }
     }
     return;
# endif // CPPAD_HAS_COLPACK
}

} } // END_CPPAD_LOCAL_NAMESPACE

# endif