hess_color.cc

Go to the documentation of this file.

// Hessian coloring algorithms -*- C++ -*-

// Copyright (C) 2008 Mihaly Csaba Markot
//
// This file is a COCONUT inference engine.  This library
// is free software; you can redistribute it and/or modify it under the
// terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2, or (at your option)
// any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// Library GNU General Public License for more details.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the Library GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the Library GNU General Public License.

#include <hess_color.h>
#include <print_seq.h>

namespace coco
{

#define HESS_COLOR_DEBUG 0

void hess_star_color1(const unsigned int n,
                      const linalg::matrix<int>& H,
                      std::vector<std::vector<unsigned int> >& parts,
                      std::vector<linalg::sparse_vector<unsigned int> >& provides)
{
// initialize working vectors
  std::vector<unsigned int> color(n,n);  // colors are indexed from 0..n-1,
                                         // color[i] == n means that H[i] is
                                         // not colored yet

  std::vector<unsigned int> forb_colors(n,n); // forb_colors[j] = i means that
                                              // color j is forbidden to the
                                              // vertex H[i]
  unsigned int v, w, x, y, c;
  linalg::matrix<int>::OneD::const_iterator itv, itw, itx;
  int maxcolor(-1);  // highest index of used colors (number of required
                     // colors == maxcolor+1)

  for (v=0; v<n; ++v)
  {
#if HESS_COLOR_DEBUG > 1
    std::cout << "hess_star_color1: processing vertex " << v << std::endl;
#endif
    for (itv=H[v].begin(); itv<H[v].end(); ++itv) // iterate through N_1[v]
    {
      w = itv.column();
      if (color[w] != n)  // w is colored
        forb_colors[color[w]] = v;
      for (itw=H[w].begin(); itw<H[w].end(); ++itw) // iterate through N_1[w]
      {
        x = itw.column();
        if (color[x] == n) continue;
        if (color[w] == n)  // ww is not colored
          forb_colors[color[x]] = v;
        else
        {
          for (itx=H[x].begin(); itx<H[x].end(); ++itx) // iterate through N_1[x]
          {
            y = itx.column();
            if (y == w || color[y] == n) continue;
            if (color[y] == color[w])
            {
              forb_colors[color[x]] = v;
              break;
            }
          } // for itx
        }
      } // for itw
    } // for itv

    // assign to H[v] the smallest non-forbidden color
    for (c=0; c<n; ++c)
    {
      if (forb_colors[c] != v)
      {
        color[v] = c;
#if HESS_COLOR_DEBUG > 1
        std::cout << "hess_star_color1: setting color[" << v << "] = "
                  << c << std::endl;
#endif
        if ((int)c > maxcolor)
          maxcolor = c;
        break;
      }
    }
    if (c==n) // no free color has been found => error
      throw api_exception(apiee_internal, "internal error in hess_star_color1");
  } // for v

  // create parts by listing all the same colored rows/columns of the Hessian,
  // i.e. parts[i] will consist of all the row/column indices which has color i
  parts.clear();
  for (int ii=0; ii<=maxcolor; ++ii)
  {
    std::vector<unsigned int> vec;
    parts.push_back(vec);
  }
  for (v=0; v<n; ++v)
    parts[color[v]].push_back(v);

  unsigned int i,j,k;
  unsigned int cnt, ind;

#if HESS_COLOR_DEBUG
  std::cout << "hess_star_color1: maxcolor = " << maxcolor << std::endl;
  std::cout << "hess_star_color1: parts:" << std::endl;
  for (i=0; i<=(unsigned int)maxcolor; ++i)
    std::cout << parts[i] << std::endl;
#endif
  
  for (i=0; i<parts.size(); ++i)  // check partition i (ith binary vector)
  {
    provides.push_back(linalg::sparse_vector<unsigned int>());
    for (j=0; j<n; ++j)  // for all rows of H
    {
      cnt = 0; ind = 0;
      for (k=0; k<parts[i].size(); k++)
      {
        if (H(j,parts[i][k]) != 0)
        {
          cnt++;
          ind = parts[i][k];
          if (cnt > 1)  // not a unique nonzero
            break;
        }
      }
      if (cnt == 1)
        provides[i][j] = ind+1;
    }
  }
#if HESS_COLOR_DEBUG
  std::cout << "hess_color_eval_provides: provides = " << std::endl;
  std::cout << provides << std::endl;
  // mtl::print_all_matrix(provides);
#endif
}

} // namespace ipopt

---