epoc32/include/stdapis/boost/graph/minimum_degree_ordering.hpp
branchSymbian2
changeset 2 2fe1408b6811
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/epoc32/include/stdapis/boost/graph/minimum_degree_ordering.hpp	Tue Mar 16 16:12:26 2010 +0000
     1.3 @@ -0,0 +1,655 @@
     1.4 +//-*-c++-*-
     1.5 +//=======================================================================
     1.6 +// Copyright 1997-2001 University of Notre Dame.
     1.7 +// Authors: Lie-Quan Lee, Jeremy Siek
     1.8 +//
     1.9 +// Distributed under the Boost Software License, Version 1.0. (See
    1.10 +// accompanying file LICENSE_1_0.txt or copy at
    1.11 +// http://www.boost.org/LICENSE_1_0.txt)
    1.12 +//=======================================================================
    1.13 +//
    1.14 +#ifndef MINIMUM_DEGREE_ORDERING_HPP
    1.15 +#define MINIMUM_DEGREE_ORDERING_HPP
    1.16 +
    1.17 +#include <vector>
    1.18 +#include <cassert>
    1.19 +#include <boost/config.hpp>
    1.20 +#include <boost/pending/bucket_sorter.hpp>
    1.21 +#include <boost/detail/numeric_traits.hpp> // for integer_traits
    1.22 +#include <boost/graph/graph_traits.hpp>
    1.23 +#include <boost/property_map.hpp>
    1.24 +
    1.25 +namespace boost {
    1.26 +
    1.27 +  namespace detail {
    1.28 +
    1.29 +    // 
    1.30 +    // Given a set of n integers (where the integer values range from
    1.31 +    // zero to n-1), we want to keep track of a collection of stacks
    1.32 +    // of integers. It so happens that an integer will appear in at
    1.33 +    // most one stack at a time, so the stacks form disjoint sets.
    1.34 +    // Because of these restrictions, we can use one big array to
    1.35 +    // store all the stacks, intertwined with one another.
    1.36 +    // No allocation/deallocation happens in the push()/pop() methods
    1.37 +    // so this is faster than using std::stack's.
    1.38 +    //
    1.39 +    template <class SignedInteger>
    1.40 +    class Stacks {
    1.41 +      typedef SignedInteger value_type;
    1.42 +      typedef typename std::vector<value_type>::size_type size_type;
    1.43 +    public:
    1.44 +      Stacks(size_type n) : data(n) {}
    1.45 +      
    1.46 +      //: stack 
    1.47 +      class stack {
    1.48 +        typedef typename std::vector<value_type>::iterator Iterator;
    1.49 +      public:
    1.50 +        stack(Iterator _data, const value_type& head)
    1.51 +          :  data(_data), current(head) {}
    1.52 +
    1.53 +        // did not use default argument here to avoid internal compiler error
    1.54 +        // in g++.
    1.55 +        stack(Iterator _data)
    1.56 +          : data(_data), current(-(std::numeric_limits<value_type>::max)()) {}
    1.57 +        
    1.58 +        void pop() {
    1.59 +          assert(! empty());
    1.60 +          current = data[current];
    1.61 +        }
    1.62 +        void push(value_type v) {
    1.63 +          data[v] = current; 
    1.64 +          current = v;
    1.65 +        }
    1.66 +        bool empty() {
    1.67 +          return current == -(std::numeric_limits<value_type>::max)(); 
    1.68 +        }
    1.69 +        value_type& top() { return current; }
    1.70 +      private:
    1.71 +        Iterator data;
    1.72 +        value_type current;
    1.73 +      };
    1.74 +
    1.75 +      // To return a stack object 
    1.76 +      stack make_stack()
    1.77 +        { return stack(data.begin()); }
    1.78 +    protected:
    1.79 +      std::vector<value_type> data;
    1.80 +    };
    1.81 +
    1.82 +
    1.83 +    // marker class, a generalization of coloring. 
    1.84 +    //
    1.85 +    // This class is to provide a generalization of coloring which has
    1.86 +    // complexity of amortized constant time to set all vertices' color
    1.87 +    // back to be untagged. It implemented by increasing a tag.
    1.88 +    //
    1.89 +    // The colors are:
    1.90 +    //   not tagged 
    1.91 +    //   tagged
    1.92 +    //   multiple_tagged
    1.93 +    //   done
    1.94 +    //
    1.95 +    template <class SignedInteger, class Vertex, class VertexIndexMap>
    1.96 +    class Marker {
    1.97 +      typedef SignedInteger value_type;
    1.98 +      typedef typename std::vector<value_type>::size_type size_type;
    1.99 +      
   1.100 +      static value_type done() 
   1.101 +      { return (std::numeric_limits<value_type>::max)()/2; }
   1.102 +    public:
   1.103 +      Marker(size_type _num, VertexIndexMap index_map) 
   1.104 +        : tag(1 - (std::numeric_limits<value_type>::max)()),
   1.105 +          data(_num, - (std::numeric_limits<value_type>::max)()),
   1.106 +          id(index_map) {}
   1.107 +      
   1.108 +      void mark_done(Vertex node) { data[get(id, node)] = done(); }
   1.109 +      
   1.110 +      bool is_done(Vertex node) { return data[get(id, node)] == done(); }
   1.111 +      
   1.112 +      void mark_tagged(Vertex node) { data[get(id, node)] = tag; }
   1.113 +      
   1.114 +      void mark_multiple_tagged(Vertex node) { data[get(id, node)] = multiple_tag; }
   1.115 +  
   1.116 +      bool is_tagged(Vertex node) const { return data[get(id, node)] >= tag; }
   1.117 +
   1.118 +      bool is_not_tagged(Vertex node) const { return data[get(id, node)] < tag; }
   1.119 +
   1.120 +      bool is_multiple_tagged(Vertex node) const 
   1.121 +        { return data[get(id, node)] >= multiple_tag; }
   1.122 +
   1.123 +      void increment_tag() {
   1.124 +        const size_type num = data.size();
   1.125 +        ++tag;
   1.126 +        if ( tag >= done() ) {
   1.127 +          tag = 1 - (std::numeric_limits<value_type>::max)();
   1.128 +          for (size_type i = 0; i < num; ++i)
   1.129 +            if ( data[i] < done() ) 
   1.130 +              data[i] = - (std::numeric_limits<value_type>::max)();
   1.131 +        }
   1.132 +      }
   1.133 +      
   1.134 +      void set_multiple_tag(value_type mdeg0) 
   1.135 +      { 
   1.136 +        const size_type num = data.size();
   1.137 +        multiple_tag = tag + mdeg0; 
   1.138 +        
   1.139 +        if ( multiple_tag >= done() ) {
   1.140 +          tag = 1-(std::numeric_limits<value_type>::max)();
   1.141 +          
   1.142 +          for (size_type i=0; i<num; i++)
   1.143 +            if ( data[i] < done() ) 
   1.144 +              data[i] = -(std::numeric_limits<value_type>::max)();
   1.145 +          
   1.146 +          multiple_tag = tag + mdeg0; 
   1.147 +        }
   1.148 +      }
   1.149 +      
   1.150 +      void set_tag_as_multiple_tag() { tag = multiple_tag; }
   1.151 +      
   1.152 +    protected:
   1.153 +      value_type tag;
   1.154 +      value_type multiple_tag;
   1.155 +      std::vector<value_type> data;
   1.156 +      VertexIndexMap id;
   1.157 +    };
   1.158 +    
   1.159 +    template< class Iterator, class SignedInteger, 
   1.160 +       class Vertex, class VertexIndexMap, int offset = 1 >
   1.161 +    class Numbering {
   1.162 +      typedef SignedInteger number_type;
   1.163 +      number_type num; //start from 1 instead of zero
   1.164 +      Iterator   data;
   1.165 +      number_type max_num;
   1.166 +      VertexIndexMap id;
   1.167 +    public:
   1.168 +      Numbering(Iterator _data, number_type _max_num, VertexIndexMap id) 
   1.169 +        : num(1), data(_data), max_num(_max_num), id(id) {}
   1.170 +      void operator()(Vertex node) { data[get(id, node)] = -num; }
   1.171 +      bool all_done(number_type i = 0) const { return num + i > max_num; }
   1.172 +      void increment(number_type i = 1) { num += i; }
   1.173 +      bool is_numbered(Vertex node) const {
   1.174 +        return data[get(id, node)] < 0;
   1.175 +      }
   1.176 +      void indistinguishable(Vertex i, Vertex j) {
   1.177 +        data[get(id, i)] = - (get(id, j) + offset);
   1.178 +      }
   1.179 +    };
   1.180 +
   1.181 +    template <class SignedInteger, class Vertex, class VertexIndexMap>
   1.182 +    class degreelists_marker {
   1.183 +    public:
   1.184 +      typedef SignedInteger value_type;
   1.185 +      typedef typename std::vector<value_type>::size_type size_type;
   1.186 +      degreelists_marker(size_type n, VertexIndexMap id)
   1.187 +        : marks(n, 0), id(id) {}
   1.188 +      void mark_need_update(Vertex i) { marks[get(id, i)] = 1;  }
   1.189 +      bool need_update(Vertex i) { return marks[get(id, i)] == 1; }
   1.190 +      bool outmatched_or_done (Vertex i) { return marks[get(id, i)] == -1; }
   1.191 +      void mark(Vertex i) { marks[get(id, i)] = -1; }
   1.192 +      void unmark(Vertex i) { marks[get(id, i)] = 0; }
   1.193 +    private:
   1.194 +      std::vector<value_type> marks;
   1.195 +      VertexIndexMap id;
   1.196 +    };
   1.197 +
   1.198 +    // Helper function object for edge removal
   1.199 +    template <class Graph, class MarkerP, class NumberD, class Stack,
   1.200 +      class VertexIndexMap>
   1.201 +    class predicateRemoveEdge1 {
   1.202 +      typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
   1.203 +      typedef typename graph_traits<Graph>::edge_descriptor edge_t;
   1.204 +    public:
   1.205 +      predicateRemoveEdge1(Graph& _g, MarkerP& _marker, 
   1.206 +                           NumberD _numbering, Stack& n_e, VertexIndexMap id)
   1.207 +        : g(&_g), marker(&_marker), numbering(_numbering),
   1.208 +          neighbor_elements(&n_e), id(id) {}
   1.209 +
   1.210 +      bool operator()(edge_t e) {
   1.211 +        vertex_t dist = target(e, *g);
   1.212 +        if ( marker->is_tagged(dist) )
   1.213 +          return true;
   1.214 +        marker->mark_tagged(dist);
   1.215 +        if (numbering.is_numbered(dist)) {
   1.216 +          neighbor_elements->push(get(id, dist));
   1.217 +          return true;
   1.218 +        }
   1.219 +        return false;
   1.220 +      }
   1.221 +    private:
   1.222 +      Graph*   g;
   1.223 +      MarkerP* marker;
   1.224 +      NumberD  numbering;
   1.225 +      Stack*   neighbor_elements;
   1.226 +      VertexIndexMap id;
   1.227 +    };
   1.228 +
   1.229 +    // Helper function object for edge removal
   1.230 +    template <class Graph, class MarkerP>
   1.231 +    class predicate_remove_tagged_edges
   1.232 +    {
   1.233 +      typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
   1.234 +      typedef typename graph_traits<Graph>::edge_descriptor edge_t;
   1.235 +    public:
   1.236 +      predicate_remove_tagged_edges(Graph& _g, MarkerP& _marker)
   1.237 +        : g(&_g), marker(&_marker) {}
   1.238 +
   1.239 +      bool operator()(edge_t e) {
   1.240 +        vertex_t dist = target(e, *g);
   1.241 +        if ( marker->is_tagged(dist) )
   1.242 +          return true;
   1.243 +        return false;
   1.244 +      }
   1.245 +    private:
   1.246 +      Graph*   g;
   1.247 +      MarkerP* marker;
   1.248 +    };
   1.249 +
   1.250 +    template<class Graph, class DegreeMap, 
   1.251 +             class InversePermutationMap, 
   1.252 +             class PermutationMap,
   1.253 +             class SuperNodeMap, 
   1.254 +             class VertexIndexMap>
   1.255 +    class mmd_impl
   1.256 +    {
   1.257 +      // Typedefs
   1.258 +      typedef graph_traits<Graph> Traits;
   1.259 +      typedef typename Traits::vertices_size_type size_type;
   1.260 +      typedef typename detail::integer_traits<size_type>::difference_type 
   1.261 +        diff_t;
   1.262 +      typedef typename Traits::vertex_descriptor vertex_t;
   1.263 +      typedef typename Traits::adjacency_iterator adj_iter;
   1.264 +      typedef iterator_property_map<vertex_t*, 
   1.265 +        identity_property_map, vertex_t, vertex_t&> IndexVertexMap;
   1.266 +      typedef detail::Stacks<diff_t> Workspace;
   1.267 +      typedef bucket_sorter<size_type, vertex_t, DegreeMap, VertexIndexMap> 
   1.268 +        DegreeLists;
   1.269 +      typedef Numbering<InversePermutationMap, diff_t, vertex_t,VertexIndexMap>
   1.270 +        NumberingD;
   1.271 +      typedef degreelists_marker<diff_t, vertex_t, VertexIndexMap> 
   1.272 +        DegreeListsMarker;
   1.273 +      typedef Marker<diff_t, vertex_t, VertexIndexMap> MarkerP;
   1.274 +
   1.275 +      // Data Members
   1.276 +
   1.277 +      // input parameters
   1.278 +      Graph& G;
   1.279 +      int delta;
   1.280 +      DegreeMap degree;
   1.281 +      InversePermutationMap inverse_perm;
   1.282 +      PermutationMap perm;
   1.283 +      SuperNodeMap supernode_size;
   1.284 +      VertexIndexMap vertex_index_map;
   1.285 +
   1.286 +      // internal data-structures
   1.287 +      std::vector<vertex_t> index_vertex_vec;
   1.288 +      size_type n;
   1.289 +      IndexVertexMap index_vertex_map;
   1.290 +      DegreeLists degreelists;
   1.291 +      NumberingD numbering;
   1.292 +      DegreeListsMarker degree_lists_marker;
   1.293 +      MarkerP marker;
   1.294 +      Workspace work_space;
   1.295 +    public:
   1.296 +      mmd_impl(Graph& g, size_type n_, int delta, DegreeMap degree, 
   1.297 +               InversePermutationMap inverse_perm, 
   1.298 +               PermutationMap perm,
   1.299 +               SuperNodeMap supernode_size, 
   1.300 +               VertexIndexMap id) 
   1.301 +        : G(g), delta(delta), degree(degree), 
   1.302 +        inverse_perm(inverse_perm), 
   1.303 +        perm(perm), 
   1.304 +        supernode_size(supernode_size), 
   1.305 +        vertex_index_map(id),
   1.306 +        index_vertex_vec(n_), 
   1.307 +        n(n_),
   1.308 +        degreelists(n_ + 1, n_, degree, id),
   1.309 +        numbering(inverse_perm, n_, vertex_index_map),
   1.310 +        degree_lists_marker(n_, vertex_index_map), 
   1.311 +        marker(n_, vertex_index_map),
   1.312 +        work_space(n_)
   1.313 +      {
   1.314 +        typename graph_traits<Graph>::vertex_iterator v, vend;
   1.315 +        size_type vid = 0;
   1.316 +        for (tie(v, vend) = vertices(G); v != vend; ++v, ++vid)
   1.317 +          index_vertex_vec[vid] = *v;
   1.318 +        index_vertex_map = IndexVertexMap(&index_vertex_vec[0]);
   1.319 +
   1.320 +        // Initialize degreelists.  Degreelists organizes the nodes
   1.321 +        // according to their degree.
   1.322 +        for (tie(v, vend) = vertices(G); v != vend; ++v) {
   1.323 +          put(degree, *v, out_degree(*v, G));
   1.324 +          degreelists.push(*v);
   1.325 +        }
   1.326 +      }
   1.327 +
   1.328 +      void do_mmd()
   1.329 +      {
   1.330 +        // Eliminate the isolated nodes -- these are simply the nodes
   1.331 +        // with no neighbors, which are accessible as a list (really, a
   1.332 +        // stack) at location 0.  Since these don't affect any other
   1.333 +        // nodes, we can eliminate them without doing degree updates.
   1.334 +        typename DegreeLists::stack list_isolated = degreelists[0];
   1.335 +        while (!list_isolated.empty()) {
   1.336 +          vertex_t node = list_isolated.top();
   1.337 +          marker.mark_done(node);
   1.338 +          numbering(node);
   1.339 +          numbering.increment();
   1.340 +          list_isolated.pop();
   1.341 +        }
   1.342 +        size_type min_degree = 1;
   1.343 +        typename DegreeLists::stack list_min_degree = degreelists[min_degree];
   1.344 +
   1.345 +        while (list_min_degree.empty()) {
   1.346 +          ++min_degree;
   1.347 +          list_min_degree = degreelists[min_degree];
   1.348 +        }
   1.349 +
   1.350 +        // check if the whole eliminating process is done
   1.351 +        while (!numbering.all_done()) {
   1.352 +
   1.353 +          size_type min_degree_limit = min_degree + delta; // WARNING
   1.354 +          typename Workspace::stack llist = work_space.make_stack();
   1.355 +
   1.356 +          // multiple elimination
   1.357 +          while (delta >= 0) {
   1.358 +
   1.359 +            // Find the next non-empty degree
   1.360 +            for (list_min_degree = degreelists[min_degree]; 
   1.361 +                 list_min_degree.empty() && min_degree <= min_degree_limit; 
   1.362 +                 ++min_degree, list_min_degree = degreelists[min_degree])
   1.363 +              ;
   1.364 +            if (min_degree > min_degree_limit)
   1.365 +              break;
   1.366 +
   1.367 +            const vertex_t node = list_min_degree.top();
   1.368 +            const size_type node_id = get(vertex_index_map, node);
   1.369 +            list_min_degree.pop();
   1.370 +            numbering(node);
   1.371 +
   1.372 +            // check if node is the last one
   1.373 +            if (numbering.all_done(supernode_size[node])) {
   1.374 +              numbering.increment(supernode_size[node]);
   1.375 +              break;
   1.376 +            }
   1.377 +            marker.increment_tag();
   1.378 +            marker.mark_tagged(node);
   1.379 +
   1.380 +            this->eliminate(node);
   1.381 +
   1.382 +            numbering.increment(supernode_size[node]);
   1.383 +            llist.push(node_id);
   1.384 +          } // multiple elimination
   1.385 +
   1.386 +          if (numbering.all_done()) 
   1.387 +            break;
   1.388 +
   1.389 +          this->update( llist, min_degree);
   1.390 +        }
   1.391 +
   1.392 +      } // do_mmd()
   1.393 +
   1.394 +      void eliminate(vertex_t node)
   1.395 +      {
   1.396 +        typename Workspace::stack element_neighbor = work_space.make_stack();
   1.397 +
   1.398 +        // Create two function objects for edge removal
   1.399 +        typedef typename Workspace::stack WorkStack;
   1.400 +        predicateRemoveEdge1<Graph, MarkerP, NumberingD, 
   1.401 +                             WorkStack, VertexIndexMap>
   1.402 +          p(G, marker, numbering, element_neighbor, vertex_index_map);
   1.403 +
   1.404 +        predicate_remove_tagged_edges<Graph, MarkerP> p2(G, marker);
   1.405 +
   1.406 +        // Reconstruct the adjacent node list, push element neighbor in a List.
   1.407 +        remove_out_edge_if(node, p, G);
   1.408 +        //during removal element neighbors are collected.
   1.409 +
   1.410 +        while (!element_neighbor.empty()) {
   1.411 +          // element absorb
   1.412 +          size_type e_id = element_neighbor.top();
   1.413 +          vertex_t element = get(index_vertex_map, e_id);
   1.414 +          adj_iter i, i_end;
   1.415 +          for (tie(i, i_end) = adjacent_vertices(element, G); i != i_end; ++i){
   1.416 +            vertex_t i_node = *i;
   1.417 +            if (!marker.is_tagged(i_node) && !numbering.is_numbered(i_node)) {
   1.418 +              marker.mark_tagged(i_node);
   1.419 +              add_edge(node, i_node, G);
   1.420 +            }
   1.421 +          }
   1.422 +          element_neighbor.pop();
   1.423 +        }
   1.424 +        adj_iter v, ve;
   1.425 +        for (tie(v, ve) = adjacent_vertices(node, G); v != ve; ++v) {
   1.426 +          vertex_t v_node = *v;
   1.427 +          if (!degree_lists_marker.need_update(v_node) 
   1.428 +              && !degree_lists_marker.outmatched_or_done(v_node)) {
   1.429 +            degreelists.remove(v_node);
   1.430 +          }
   1.431 +          //update out edges of v_node
   1.432 +          remove_out_edge_if(v_node, p2, G);
   1.433 +
   1.434 +          if ( out_degree(v_node, G) == 0 ) { // indistinguishable nodes
   1.435 +            supernode_size[node] += supernode_size[v_node];
   1.436 +            supernode_size[v_node] = 0;
   1.437 +            numbering.indistinguishable(v_node, node);
   1.438 +            marker.mark_done(v_node);
   1.439 +            degree_lists_marker.mark(v_node);
   1.440 +          } else {                           // not indistinguishable nodes
   1.441 +            add_edge(v_node, node, G);
   1.442 +            degree_lists_marker.mark_need_update(v_node);
   1.443 +          }
   1.444 +        }
   1.445 +      } // eliminate()
   1.446 +
   1.447 +
   1.448 +      template <class Stack>
   1.449 +      void update(Stack llist, size_type& min_degree)
   1.450 +      {
   1.451 +        size_type min_degree0 = min_degree + delta + 1;
   1.452 +
   1.453 +        while (! llist.empty()) {
   1.454 +          size_type deg, deg0 = 0;
   1.455 +
   1.456 +          marker.set_multiple_tag(min_degree0);
   1.457 +          typename Workspace::stack q2list = work_space.make_stack();
   1.458 +          typename Workspace::stack qxlist = work_space.make_stack();
   1.459 +
   1.460 +          vertex_t current = get(index_vertex_map, llist.top());
   1.461 +          adj_iter i, ie;
   1.462 +          for (tie(i,ie) = adjacent_vertices(current, G); i != ie; ++i) {
   1.463 +            vertex_t i_node = *i;
   1.464 +            const size_type i_id = get(vertex_index_map, i_node);
   1.465 +            if (supernode_size[i_node] != 0) {
   1.466 +              deg0 += supernode_size[i_node];
   1.467 +              marker.mark_multiple_tagged(i_node);
   1.468 +              if (degree_lists_marker.need_update(i_node)) {
   1.469 +                if (out_degree(i_node, G) == 2)
   1.470 +                  q2list.push(i_id);
   1.471 +                else
   1.472 +                  qxlist.push(i_id);
   1.473 +              }
   1.474 +            }
   1.475 +          }
   1.476 +
   1.477 +          while (!q2list.empty()) {
   1.478 +            const size_type u_id = q2list.top();
   1.479 +            vertex_t u_node = get(index_vertex_map, u_id);
   1.480 +            // if u_id is outmatched by others, no need to update degree
   1.481 +            if (degree_lists_marker.outmatched_or_done(u_node)) {
   1.482 +              q2list.pop();
   1.483 +              continue;
   1.484 +            }
   1.485 +            marker.increment_tag();
   1.486 +            deg = deg0;
   1.487 +
   1.488 +            adj_iter nu = adjacent_vertices(u_node, G).first;
   1.489 +            vertex_t neighbor = *nu;
   1.490 +            if (neighbor == u_node) {
   1.491 +              ++nu;
   1.492 +              neighbor = *nu;
   1.493 +            }
   1.494 +            if (numbering.is_numbered(neighbor)) {
   1.495 +              adj_iter i, ie;
   1.496 +              for (tie(i,ie) = adjacent_vertices(neighbor, G);
   1.497 +                   i != ie; ++i) {
   1.498 +                const vertex_t i_node = *i;
   1.499 +                if (i_node == u_node || supernode_size[i_node] == 0)
   1.500 +                  continue;
   1.501 +                if (marker.is_tagged(i_node)) {
   1.502 +                  if (degree_lists_marker.need_update(i_node)) {
   1.503 +                    if ( out_degree(i_node, G) == 2 ) { // is indistinguishable
   1.504 +                      supernode_size[u_node] += supernode_size[i_node];
   1.505 +                      supernode_size[i_node] = 0;
   1.506 +                      numbering.indistinguishable(i_node, u_node);
   1.507 +                      marker.mark_done(i_node);
   1.508 +                      degree_lists_marker.mark(i_node);
   1.509 +                    } else                        // is outmatched
   1.510 +                      degree_lists_marker.mark(i_node);
   1.511 +                  }
   1.512 +                } else {
   1.513 +                  marker.mark_tagged(i_node);
   1.514 +                  deg += supernode_size[i_node];
   1.515 +                }
   1.516 +              }
   1.517 +            } else
   1.518 +              deg += supernode_size[neighbor];
   1.519 +
   1.520 +            deg -= supernode_size[u_node];
   1.521 +            degree[u_node] = deg; //update degree
   1.522 +            degreelists[deg].push(u_node);
   1.523 +            //u_id has been pushed back into degreelists
   1.524 +            degree_lists_marker.unmark(u_node);
   1.525 +            if (min_degree > deg) 
   1.526 +              min_degree = deg;
   1.527 +            q2list.pop();
   1.528 +          } // while (!q2list.empty())
   1.529 +
   1.530 +          while (!qxlist.empty()) {
   1.531 +            const size_type u_id = qxlist.top();
   1.532 +            const vertex_t u_node = get(index_vertex_map, u_id);
   1.533 +
   1.534 +            // if u_id is outmatched by others, no need to update degree
   1.535 +            if (degree_lists_marker.outmatched_or_done(u_node)) {
   1.536 +              qxlist.pop();
   1.537 +              continue;
   1.538 +            }
   1.539 +            marker.increment_tag();
   1.540 +            deg = deg0;
   1.541 +            adj_iter i, ie;
   1.542 +            for (tie(i, ie) = adjacent_vertices(u_node, G); i != ie; ++i) {
   1.543 +              vertex_t i_node = *i;
   1.544 +              if (marker.is_tagged(i_node)) 
   1.545 +                continue;
   1.546 +              marker.mark_tagged(i_node);
   1.547 +
   1.548 +              if (numbering.is_numbered(i_node)) {
   1.549 +                adj_iter j, je;
   1.550 +                for (tie(j, je) = adjacent_vertices(i_node, G); j != je; ++j) {
   1.551 +                  const vertex_t j_node = *j;
   1.552 +                  if (marker.is_not_tagged(j_node)) {
   1.553 +                    marker.mark_tagged(j_node);
   1.554 +                    deg += supernode_size[j_node];
   1.555 +                  }
   1.556 +                }
   1.557 +              } else
   1.558 +                deg += supernode_size[i_node];
   1.559 +            } // for adjacent vertices of u_node
   1.560 +            deg -= supernode_size[u_node];
   1.561 +            degree[u_node] = deg;
   1.562 +            degreelists[deg].push(u_node);
   1.563 +            // u_id has been pushed back into degreelists
   1.564 +            degree_lists_marker.unmark(u_node); 
   1.565 +            if (min_degree > deg)
   1.566 +              min_degree = deg;
   1.567 +            qxlist.pop();
   1.568 +          } // while (!qxlist.empty()) {
   1.569 +
   1.570 +          marker.set_tag_as_multiple_tag();
   1.571 +          llist.pop();
   1.572 +        } // while (! llist.empty())
   1.573 +
   1.574 +      } // update()
   1.575 +
   1.576 +
   1.577 +      void build_permutation(InversePermutationMap next,
   1.578 +                             PermutationMap prev) 
   1.579 +      {
   1.580 +        // collect the permutation info
   1.581 +        size_type i;
   1.582 +        for (i = 0; i < n; ++i) {
   1.583 +          diff_t size = supernode_size[get(index_vertex_map, i)];
   1.584 +          if ( size <= 0 ) {
   1.585 +            prev[i] = next[i];
   1.586 +            supernode_size[get(index_vertex_map, i)]
   1.587 +              = next[i] + 1;  // record the supernode info
   1.588 +          } else
   1.589 +            prev[i] = - next[i];
   1.590 +        }
   1.591 +        for (i = 1; i < n + 1; ++i) {
   1.592 +          if ( prev[i-1] > 0 )
   1.593 +            continue;
   1.594 +          diff_t parent = i;
   1.595 +          while ( prev[parent - 1] < 0 ) {
   1.596 +            parent = - prev[parent - 1];
   1.597 +          }
   1.598 +
   1.599 +          diff_t root = parent;
   1.600 +          diff_t num = prev[root - 1] + 1;
   1.601 +          next[i-1] = - num;
   1.602 +          prev[root-1] = num;
   1.603 +
   1.604 +          parent = i;
   1.605 +          diff_t next_node = - prev[parent - 1];
   1.606 +          while (next_node > 0) {
   1.607 +            prev[parent-1] = - root;
   1.608 +            parent = next_node;
   1.609 +            next_node = - prev[parent - 1];
   1.610 +          }
   1.611 +        }
   1.612 +        for (i = 0; i < n; i++) {
   1.613 +          diff_t num = - next[i] - 1;
   1.614 +          next[i] = num;
   1.615 +          prev[num] = i;
   1.616 +        }
   1.617 +      } // build_permutation()
   1.618 +    };
   1.619 +
   1.620 +  } //namespace detail
   1.621 +
   1.622 +
   1.623 +  // MMD algorithm
   1.624 +  // 
   1.625 +  //The implementation presently includes the enhancements for mass
   1.626 +  //elimination, incomplete degree update, multiple elimination, and
   1.627 +  //external degree.
   1.628 +  //
   1.629 +  //Important Note: This implementation requires the BGL graph to be
   1.630 +  //directed.  Therefore, nonzero entry (i, j) in a symmetrical matrix
   1.631 +  //A coresponds to two directed edges (i->j and j->i).
   1.632 +  //
   1.633 +  //see Alan George and Joseph W. H. Liu, The Evolution of the Minimum
   1.634 +  //Degree Ordering Algorithm, SIAM Review, 31, 1989, Page 1-19
   1.635 +  template<class Graph, class DegreeMap, 
   1.636 +           class InversePermutationMap, 
   1.637 +           class PermutationMap,
   1.638 +           class SuperNodeMap, class VertexIndexMap>
   1.639 +  void minimum_degree_ordering
   1.640 +    (Graph& G, 
   1.641 +     DegreeMap degree, 
   1.642 +     InversePermutationMap inverse_perm, 
   1.643 +     PermutationMap perm, 
   1.644 +     SuperNodeMap supernode_size, 
   1.645 +     int delta, 
   1.646 +     VertexIndexMap vertex_index_map)
   1.647 +  {
   1.648 +    detail::mmd_impl<Graph,DegreeMap,InversePermutationMap,
   1.649 +      PermutationMap, SuperNodeMap, VertexIndexMap> 
   1.650 +      impl(G, num_vertices(G), delta, degree, inverse_perm, 
   1.651 +           perm, supernode_size, vertex_index_map);
   1.652 +    impl.do_mmd();
   1.653 +    impl.build_permutation(inverse_perm, perm);
   1.654 +  }
   1.655 +
   1.656 +} // namespace boost
   1.657 +
   1.658 +#endif // MINIMUM_DEGREE_ORDERING_HPP