williamr@2: //-*-c++-*-
williamr@2: //=======================================================================
williamr@2: // Copyright 1997-2001 University of Notre Dame.
williamr@2: // Authors: Lie-Quan Lee, Jeremy Siek
williamr@2: //
williamr@2: // Distributed under the Boost Software License, Version 1.0. (See
williamr@2: // accompanying file LICENSE_1_0.txt or copy at
williamr@2: // http://www.boost.org/LICENSE_1_0.txt)
williamr@2: //=======================================================================
williamr@2: //
williamr@2: #ifndef MINIMUM_DEGREE_ORDERING_HPP
williamr@2: #define MINIMUM_DEGREE_ORDERING_HPP
williamr@2: 
williamr@2: #include <vector>
williamr@2: #include <cassert>
williamr@2: #include <boost/config.hpp>
williamr@2: #include <boost/pending/bucket_sorter.hpp>
williamr@2: #include <boost/detail/numeric_traits.hpp> // for integer_traits
williamr@2: #include <boost/graph/graph_traits.hpp>
williamr@2: #include <boost/property_map.hpp>
williamr@2: 
williamr@2: namespace boost {
williamr@2: 
williamr@2:   namespace detail {
williamr@2: 
williamr@2:     // 
williamr@2:     // Given a set of n integers (where the integer values range from
williamr@2:     // zero to n-1), we want to keep track of a collection of stacks
williamr@2:     // of integers. It so happens that an integer will appear in at
williamr@2:     // most one stack at a time, so the stacks form disjoint sets.
williamr@2:     // Because of these restrictions, we can use one big array to
williamr@2:     // store all the stacks, intertwined with one another.
williamr@2:     // No allocation/deallocation happens in the push()/pop() methods
williamr@2:     // so this is faster than using std::stack's.
williamr@2:     //
williamr@2:     template <class SignedInteger>
williamr@2:     class Stacks {
williamr@2:       typedef SignedInteger value_type;
williamr@2:       typedef typename std::vector<value_type>::size_type size_type;
williamr@2:     public:
williamr@2:       Stacks(size_type n) : data(n) {}
williamr@2:       
williamr@2:       //: stack 
williamr@2:       class stack {
williamr@2:         typedef typename std::vector<value_type>::iterator Iterator;
williamr@2:       public:
williamr@2:         stack(Iterator _data, const value_type& head)
williamr@2:           :  data(_data), current(head) {}
williamr@2: 
williamr@2:         // did not use default argument here to avoid internal compiler error
williamr@2:         // in g++.
williamr@2:         stack(Iterator _data)
williamr@2:           : data(_data), current(-(std::numeric_limits<value_type>::max)()) {}
williamr@2:         
williamr@2:         void pop() {
williamr@2:           assert(! empty());
williamr@2:           current = data[current];
williamr@2:         }
williamr@2:         void push(value_type v) {
williamr@2:           data[v] = current; 
williamr@2:           current = v;
williamr@2:         }
williamr@2:         bool empty() {
williamr@2:           return current == -(std::numeric_limits<value_type>::max)(); 
williamr@2:         }
williamr@2:         value_type& top() { return current; }
williamr@2:       private:
williamr@2:         Iterator data;
williamr@2:         value_type current;
williamr@2:       };
williamr@2: 
williamr@2:       // To return a stack object 
williamr@2:       stack make_stack()
williamr@2:         { return stack(data.begin()); }
williamr@2:     protected:
williamr@2:       std::vector<value_type> data;
williamr@2:     };
williamr@2: 
williamr@2: 
williamr@2:     // marker class, a generalization of coloring. 
williamr@2:     //
williamr@2:     // This class is to provide a generalization of coloring which has
williamr@2:     // complexity of amortized constant time to set all vertices' color
williamr@2:     // back to be untagged. It implemented by increasing a tag.
williamr@2:     //
williamr@2:     // The colors are:
williamr@2:     //   not tagged 
williamr@2:     //   tagged
williamr@2:     //   multiple_tagged
williamr@2:     //   done
williamr@2:     //
williamr@2:     template <class SignedInteger, class Vertex, class VertexIndexMap>
williamr@2:     class Marker {
williamr@2:       typedef SignedInteger value_type;
williamr@2:       typedef typename std::vector<value_type>::size_type size_type;
williamr@2:       
williamr@2:       static value_type done() 
williamr@2:       { return (std::numeric_limits<value_type>::max)()/2; }
williamr@2:     public:
williamr@2:       Marker(size_type _num, VertexIndexMap index_map) 
williamr@2:         : tag(1 - (std::numeric_limits<value_type>::max)()),
williamr@2:           data(_num, - (std::numeric_limits<value_type>::max)()),
williamr@2:           id(index_map) {}
williamr@2:       
williamr@2:       void mark_done(Vertex node) { data[get(id, node)] = done(); }
williamr@2:       
williamr@2:       bool is_done(Vertex node) { return data[get(id, node)] == done(); }
williamr@2:       
williamr@2:       void mark_tagged(Vertex node) { data[get(id, node)] = tag; }
williamr@2:       
williamr@2:       void mark_multiple_tagged(Vertex node) { data[get(id, node)] = multiple_tag; }
williamr@2:   
williamr@2:       bool is_tagged(Vertex node) const { return data[get(id, node)] >= tag; }
williamr@2: 
williamr@2:       bool is_not_tagged(Vertex node) const { return data[get(id, node)] < tag; }
williamr@2: 
williamr@2:       bool is_multiple_tagged(Vertex node) const 
williamr@2:         { return data[get(id, node)] >= multiple_tag; }
williamr@2: 
williamr@2:       void increment_tag() {
williamr@2:         const size_type num = data.size();
williamr@2:         ++tag;
williamr@2:         if ( tag >= done() ) {
williamr@2:           tag = 1 - (std::numeric_limits<value_type>::max)();
williamr@2:           for (size_type i = 0; i < num; ++i)
williamr@2:             if ( data[i] < done() ) 
williamr@2:               data[i] = - (std::numeric_limits<value_type>::max)();
williamr@2:         }
williamr@2:       }
williamr@2:       
williamr@2:       void set_multiple_tag(value_type mdeg0) 
williamr@2:       { 
williamr@2:         const size_type num = data.size();
williamr@2:         multiple_tag = tag + mdeg0; 
williamr@2:         
williamr@2:         if ( multiple_tag >= done() ) {
williamr@2:           tag = 1-(std::numeric_limits<value_type>::max)();
williamr@2:           
williamr@2:           for (size_type i=0; i<num; i++)
williamr@2:             if ( data[i] < done() ) 
williamr@2:               data[i] = -(std::numeric_limits<value_type>::max)();
williamr@2:           
williamr@2:           multiple_tag = tag + mdeg0; 
williamr@2:         }
williamr@2:       }
williamr@2:       
williamr@2:       void set_tag_as_multiple_tag() { tag = multiple_tag; }
williamr@2:       
williamr@2:     protected:
williamr@2:       value_type tag;
williamr@2:       value_type multiple_tag;
williamr@2:       std::vector<value_type> data;
williamr@2:       VertexIndexMap id;
williamr@2:     };
williamr@2:     
williamr@2:     template< class Iterator, class SignedInteger, 
williamr@2:        class Vertex, class VertexIndexMap, int offset = 1 >
williamr@2:     class Numbering {
williamr@2:       typedef SignedInteger number_type;
williamr@2:       number_type num; //start from 1 instead of zero
williamr@2:       Iterator   data;
williamr@2:       number_type max_num;
williamr@2:       VertexIndexMap id;
williamr@2:     public:
williamr@2:       Numbering(Iterator _data, number_type _max_num, VertexIndexMap id) 
williamr@2:         : num(1), data(_data), max_num(_max_num), id(id) {}
williamr@2:       void operator()(Vertex node) { data[get(id, node)] = -num; }
williamr@2:       bool all_done(number_type i = 0) const { return num + i > max_num; }
williamr@2:       void increment(number_type i = 1) { num += i; }
williamr@2:       bool is_numbered(Vertex node) const {
williamr@2:         return data[get(id, node)] < 0;
williamr@2:       }
williamr@2:       void indistinguishable(Vertex i, Vertex j) {
williamr@2:         data[get(id, i)] = - (get(id, j) + offset);
williamr@2:       }
williamr@2:     };
williamr@2: 
williamr@2:     template <class SignedInteger, class Vertex, class VertexIndexMap>
williamr@2:     class degreelists_marker {
williamr@2:     public:
williamr@2:       typedef SignedInteger value_type;
williamr@2:       typedef typename std::vector<value_type>::size_type size_type;
williamr@2:       degreelists_marker(size_type n, VertexIndexMap id)
williamr@2:         : marks(n, 0), id(id) {}
williamr@2:       void mark_need_update(Vertex i) { marks[get(id, i)] = 1;  }
williamr@2:       bool need_update(Vertex i) { return marks[get(id, i)] == 1; }
williamr@2:       bool outmatched_or_done (Vertex i) { return marks[get(id, i)] == -1; }
williamr@2:       void mark(Vertex i) { marks[get(id, i)] = -1; }
williamr@2:       void unmark(Vertex i) { marks[get(id, i)] = 0; }
williamr@2:     private:
williamr@2:       std::vector<value_type> marks;
williamr@2:       VertexIndexMap id;
williamr@2:     };
williamr@2: 
williamr@2:     // Helper function object for edge removal
williamr@2:     template <class Graph, class MarkerP, class NumberD, class Stack,
williamr@2:       class VertexIndexMap>
williamr@2:     class predicateRemoveEdge1 {
williamr@2:       typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
williamr@2:       typedef typename graph_traits<Graph>::edge_descriptor edge_t;
williamr@2:     public:
williamr@2:       predicateRemoveEdge1(Graph& _g, MarkerP& _marker, 
williamr@2:                            NumberD _numbering, Stack& n_e, VertexIndexMap id)
williamr@2:         : g(&_g), marker(&_marker), numbering(_numbering),
williamr@2:           neighbor_elements(&n_e), id(id) {}
williamr@2: 
williamr@2:       bool operator()(edge_t e) {
williamr@2:         vertex_t dist = target(e, *g);
williamr@2:         if ( marker->is_tagged(dist) )
williamr@2:           return true;
williamr@2:         marker->mark_tagged(dist);
williamr@2:         if (numbering.is_numbered(dist)) {
williamr@2:           neighbor_elements->push(get(id, dist));
williamr@2:           return true;
williamr@2:         }
williamr@2:         return false;
williamr@2:       }
williamr@2:     private:
williamr@2:       Graph*   g;
williamr@2:       MarkerP* marker;
williamr@2:       NumberD  numbering;
williamr@2:       Stack*   neighbor_elements;
williamr@2:       VertexIndexMap id;
williamr@2:     };
williamr@2: 
williamr@2:     // Helper function object for edge removal
williamr@2:     template <class Graph, class MarkerP>
williamr@2:     class predicate_remove_tagged_edges
williamr@2:     {
williamr@2:       typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
williamr@2:       typedef typename graph_traits<Graph>::edge_descriptor edge_t;
williamr@2:     public:
williamr@2:       predicate_remove_tagged_edges(Graph& _g, MarkerP& _marker)
williamr@2:         : g(&_g), marker(&_marker) {}
williamr@2: 
williamr@2:       bool operator()(edge_t e) {
williamr@2:         vertex_t dist = target(e, *g);
williamr@2:         if ( marker->is_tagged(dist) )
williamr@2:           return true;
williamr@2:         return false;
williamr@2:       }
williamr@2:     private:
williamr@2:       Graph*   g;
williamr@2:       MarkerP* marker;
williamr@2:     };
williamr@2: 
williamr@2:     template<class Graph, class DegreeMap, 
williamr@2:              class InversePermutationMap, 
williamr@2:              class PermutationMap,
williamr@2:              class SuperNodeMap, 
williamr@2:              class VertexIndexMap>
williamr@2:     class mmd_impl
williamr@2:     {
williamr@2:       // Typedefs
williamr@2:       typedef graph_traits<Graph> Traits;
williamr@2:       typedef typename Traits::vertices_size_type size_type;
williamr@2:       typedef typename detail::integer_traits<size_type>::difference_type 
williamr@2:         diff_t;
williamr@2:       typedef typename Traits::vertex_descriptor vertex_t;
williamr@2:       typedef typename Traits::adjacency_iterator adj_iter;
williamr@2:       typedef iterator_property_map<vertex_t*, 
williamr@2:         identity_property_map, vertex_t, vertex_t&> IndexVertexMap;
williamr@2:       typedef detail::Stacks<diff_t> Workspace;
williamr@2:       typedef bucket_sorter<size_type, vertex_t, DegreeMap, VertexIndexMap> 
williamr@2:         DegreeLists;
williamr@2:       typedef Numbering<InversePermutationMap, diff_t, vertex_t,VertexIndexMap>
williamr@2:         NumberingD;
williamr@2:       typedef degreelists_marker<diff_t, vertex_t, VertexIndexMap> 
williamr@2:         DegreeListsMarker;
williamr@2:       typedef Marker<diff_t, vertex_t, VertexIndexMap> MarkerP;
williamr@2: 
williamr@2:       // Data Members
williamr@2: 
williamr@2:       // input parameters
williamr@2:       Graph& G;
williamr@2:       int delta;
williamr@2:       DegreeMap degree;
williamr@2:       InversePermutationMap inverse_perm;
williamr@2:       PermutationMap perm;
williamr@2:       SuperNodeMap supernode_size;
williamr@2:       VertexIndexMap vertex_index_map;
williamr@2: 
williamr@2:       // internal data-structures
williamr@2:       std::vector<vertex_t> index_vertex_vec;
williamr@2:       size_type n;
williamr@2:       IndexVertexMap index_vertex_map;
williamr@2:       DegreeLists degreelists;
williamr@2:       NumberingD numbering;
williamr@2:       DegreeListsMarker degree_lists_marker;
williamr@2:       MarkerP marker;
williamr@2:       Workspace work_space;
williamr@2:     public:
williamr@2:       mmd_impl(Graph& g, size_type n_, int delta, DegreeMap degree, 
williamr@2:                InversePermutationMap inverse_perm, 
williamr@2:                PermutationMap perm,
williamr@2:                SuperNodeMap supernode_size, 
williamr@2:                VertexIndexMap id) 
williamr@2:         : G(g), delta(delta), degree(degree), 
williamr@2:         inverse_perm(inverse_perm), 
williamr@2:         perm(perm), 
williamr@2:         supernode_size(supernode_size), 
williamr@2:         vertex_index_map(id),
williamr@2:         index_vertex_vec(n_), 
williamr@2:         n(n_),
williamr@2:         degreelists(n_ + 1, n_, degree, id),
williamr@2:         numbering(inverse_perm, n_, vertex_index_map),
williamr@2:         degree_lists_marker(n_, vertex_index_map), 
williamr@2:         marker(n_, vertex_index_map),
williamr@2:         work_space(n_)
williamr@2:       {
williamr@2:         typename graph_traits<Graph>::vertex_iterator v, vend;
williamr@2:         size_type vid = 0;
williamr@2:         for (tie(v, vend) = vertices(G); v != vend; ++v, ++vid)
williamr@2:           index_vertex_vec[vid] = *v;
williamr@2:         index_vertex_map = IndexVertexMap(&index_vertex_vec[0]);
williamr@2: 
williamr@2:         // Initialize degreelists.  Degreelists organizes the nodes
williamr@2:         // according to their degree.
williamr@2:         for (tie(v, vend) = vertices(G); v != vend; ++v) {
williamr@2:           put(degree, *v, out_degree(*v, G));
williamr@2:           degreelists.push(*v);
williamr@2:         }
williamr@2:       }
williamr@2: 
williamr@2:       void do_mmd()
williamr@2:       {
williamr@2:         // Eliminate the isolated nodes -- these are simply the nodes
williamr@2:         // with no neighbors, which are accessible as a list (really, a
williamr@2:         // stack) at location 0.  Since these don't affect any other
williamr@2:         // nodes, we can eliminate them without doing degree updates.
williamr@2:         typename DegreeLists::stack list_isolated = degreelists[0];
williamr@2:         while (!list_isolated.empty()) {
williamr@2:           vertex_t node = list_isolated.top();
williamr@2:           marker.mark_done(node);
williamr@2:           numbering(node);
williamr@2:           numbering.increment();
williamr@2:           list_isolated.pop();
williamr@2:         }
williamr@2:         size_type min_degree = 1;
williamr@2:         typename DegreeLists::stack list_min_degree = degreelists[min_degree];
williamr@2: 
williamr@2:         while (list_min_degree.empty()) {
williamr@2:           ++min_degree;
williamr@2:           list_min_degree = degreelists[min_degree];
williamr@2:         }
williamr@2: 
williamr@2:         // check if the whole eliminating process is done
williamr@2:         while (!numbering.all_done()) {
williamr@2: 
williamr@2:           size_type min_degree_limit = min_degree + delta; // WARNING
williamr@2:           typename Workspace::stack llist = work_space.make_stack();
williamr@2: 
williamr@2:           // multiple elimination
williamr@2:           while (delta >= 0) {
williamr@2: 
williamr@2:             // Find the next non-empty degree
williamr@2:             for (list_min_degree = degreelists[min_degree]; 
williamr@2:                  list_min_degree.empty() && min_degree <= min_degree_limit; 
williamr@2:                  ++min_degree, list_min_degree = degreelists[min_degree])
williamr@2:               ;
williamr@2:             if (min_degree > min_degree_limit)
williamr@2:               break;
williamr@2: 
williamr@2:             const vertex_t node = list_min_degree.top();
williamr@2:             const size_type node_id = get(vertex_index_map, node);
williamr@2:             list_min_degree.pop();
williamr@2:             numbering(node);
williamr@2: 
williamr@2:             // check if node is the last one
williamr@2:             if (numbering.all_done(supernode_size[node])) {
williamr@2:               numbering.increment(supernode_size[node]);
williamr@2:               break;
williamr@2:             }
williamr@2:             marker.increment_tag();
williamr@2:             marker.mark_tagged(node);
williamr@2: 
williamr@2:             this->eliminate(node);
williamr@2: 
williamr@2:             numbering.increment(supernode_size[node]);
williamr@2:             llist.push(node_id);
williamr@2:           } // multiple elimination
williamr@2: 
williamr@2:           if (numbering.all_done()) 
williamr@2:             break;
williamr@2: 
williamr@2:           this->update( llist, min_degree);
williamr@2:         }
williamr@2: 
williamr@2:       } // do_mmd()
williamr@2: 
williamr@2:       void eliminate(vertex_t node)
williamr@2:       {
williamr@2:         typename Workspace::stack element_neighbor = work_space.make_stack();
williamr@2: 
williamr@2:         // Create two function objects for edge removal
williamr@2:         typedef typename Workspace::stack WorkStack;
williamr@2:         predicateRemoveEdge1<Graph, MarkerP, NumberingD, 
williamr@2:                              WorkStack, VertexIndexMap>
williamr@2:           p(G, marker, numbering, element_neighbor, vertex_index_map);
williamr@2: 
williamr@2:         predicate_remove_tagged_edges<Graph, MarkerP> p2(G, marker);
williamr@2: 
williamr@2:         // Reconstruct the adjacent node list, push element neighbor in a List.
williamr@2:         remove_out_edge_if(node, p, G);
williamr@2:         //during removal element neighbors are collected.
williamr@2: 
williamr@2:         while (!element_neighbor.empty()) {
williamr@2:           // element absorb
williamr@2:           size_type e_id = element_neighbor.top();
williamr@2:           vertex_t element = get(index_vertex_map, e_id);
williamr@2:           adj_iter i, i_end;
williamr@2:           for (tie(i, i_end) = adjacent_vertices(element, G); i != i_end; ++i){
williamr@2:             vertex_t i_node = *i;
williamr@2:             if (!marker.is_tagged(i_node) && !numbering.is_numbered(i_node)) {
williamr@2:               marker.mark_tagged(i_node);
williamr@2:               add_edge(node, i_node, G);
williamr@2:             }
williamr@2:           }
williamr@2:           element_neighbor.pop();
williamr@2:         }
williamr@2:         adj_iter v, ve;
williamr@2:         for (tie(v, ve) = adjacent_vertices(node, G); v != ve; ++v) {
williamr@2:           vertex_t v_node = *v;
williamr@2:           if (!degree_lists_marker.need_update(v_node) 
williamr@2:               && !degree_lists_marker.outmatched_or_done(v_node)) {
williamr@2:             degreelists.remove(v_node);
williamr@2:           }
williamr@2:           //update out edges of v_node
williamr@2:           remove_out_edge_if(v_node, p2, G);
williamr@2: 
williamr@2:           if ( out_degree(v_node, G) == 0 ) { // indistinguishable nodes
williamr@2:             supernode_size[node] += supernode_size[v_node];
williamr@2:             supernode_size[v_node] = 0;
williamr@2:             numbering.indistinguishable(v_node, node);
williamr@2:             marker.mark_done(v_node);
williamr@2:             degree_lists_marker.mark(v_node);
williamr@2:           } else {                           // not indistinguishable nodes
williamr@2:             add_edge(v_node, node, G);
williamr@2:             degree_lists_marker.mark_need_update(v_node);
williamr@2:           }
williamr@2:         }
williamr@2:       } // eliminate()
williamr@2: 
williamr@2: 
williamr@2:       template <class Stack>
williamr@2:       void update(Stack llist, size_type& min_degree)
williamr@2:       {
williamr@2:         size_type min_degree0 = min_degree + delta + 1;
williamr@2: 
williamr@2:         while (! llist.empty()) {
williamr@2:           size_type deg, deg0 = 0;
williamr@2: 
williamr@2:           marker.set_multiple_tag(min_degree0);
williamr@2:           typename Workspace::stack q2list = work_space.make_stack();
williamr@2:           typename Workspace::stack qxlist = work_space.make_stack();
williamr@2: 
williamr@2:           vertex_t current = get(index_vertex_map, llist.top());
williamr@2:           adj_iter i, ie;
williamr@2:           for (tie(i,ie) = adjacent_vertices(current, G); i != ie; ++i) {
williamr@2:             vertex_t i_node = *i;
williamr@2:             const size_type i_id = get(vertex_index_map, i_node);
williamr@2:             if (supernode_size[i_node] != 0) {
williamr@2:               deg0 += supernode_size[i_node];
williamr@2:               marker.mark_multiple_tagged(i_node);
williamr@2:               if (degree_lists_marker.need_update(i_node)) {
williamr@2:                 if (out_degree(i_node, G) == 2)
williamr@2:                   q2list.push(i_id);
williamr@2:                 else
williamr@2:                   qxlist.push(i_id);
williamr@2:               }
williamr@2:             }
williamr@2:           }
williamr@2: 
williamr@2:           while (!q2list.empty()) {
williamr@2:             const size_type u_id = q2list.top();
williamr@2:             vertex_t u_node = get(index_vertex_map, u_id);
williamr@2:             // if u_id is outmatched by others, no need to update degree
williamr@2:             if (degree_lists_marker.outmatched_or_done(u_node)) {
williamr@2:               q2list.pop();
williamr@2:               continue;
williamr@2:             }
williamr@2:             marker.increment_tag();
williamr@2:             deg = deg0;
williamr@2: 
williamr@2:             adj_iter nu = adjacent_vertices(u_node, G).first;
williamr@2:             vertex_t neighbor = *nu;
williamr@2:             if (neighbor == u_node) {
williamr@2:               ++nu;
williamr@2:               neighbor = *nu;
williamr@2:             }
williamr@2:             if (numbering.is_numbered(neighbor)) {
williamr@2:               adj_iter i, ie;
williamr@2:               for (tie(i,ie) = adjacent_vertices(neighbor, G);
williamr@2:                    i != ie; ++i) {
williamr@2:                 const vertex_t i_node = *i;
williamr@2:                 if (i_node == u_node || supernode_size[i_node] == 0)
williamr@2:                   continue;
williamr@2:                 if (marker.is_tagged(i_node)) {
williamr@2:                   if (degree_lists_marker.need_update(i_node)) {
williamr@2:                     if ( out_degree(i_node, G) == 2 ) { // is indistinguishable
williamr@2:                       supernode_size[u_node] += supernode_size[i_node];
williamr@2:                       supernode_size[i_node] = 0;
williamr@2:                       numbering.indistinguishable(i_node, u_node);
williamr@2:                       marker.mark_done(i_node);
williamr@2:                       degree_lists_marker.mark(i_node);
williamr@2:                     } else                        // is outmatched
williamr@2:                       degree_lists_marker.mark(i_node);
williamr@2:                   }
williamr@2:                 } else {
williamr@2:                   marker.mark_tagged(i_node);
williamr@2:                   deg += supernode_size[i_node];
williamr@2:                 }
williamr@2:               }
williamr@2:             } else
williamr@2:               deg += supernode_size[neighbor];
williamr@2: 
williamr@2:             deg -= supernode_size[u_node];
williamr@2:             degree[u_node] = deg; //update degree
williamr@2:             degreelists[deg].push(u_node);
williamr@2:             //u_id has been pushed back into degreelists
williamr@2:             degree_lists_marker.unmark(u_node);
williamr@2:             if (min_degree > deg) 
williamr@2:               min_degree = deg;
williamr@2:             q2list.pop();
williamr@2:           } // while (!q2list.empty())
williamr@2: 
williamr@2:           while (!qxlist.empty()) {
williamr@2:             const size_type u_id = qxlist.top();
williamr@2:             const vertex_t u_node = get(index_vertex_map, u_id);
williamr@2: 
williamr@2:             // if u_id is outmatched by others, no need to update degree
williamr@2:             if (degree_lists_marker.outmatched_or_done(u_node)) {
williamr@2:               qxlist.pop();
williamr@2:               continue;
williamr@2:             }
williamr@2:             marker.increment_tag();
williamr@2:             deg = deg0;
williamr@2:             adj_iter i, ie;
williamr@2:             for (tie(i, ie) = adjacent_vertices(u_node, G); i != ie; ++i) {
williamr@2:               vertex_t i_node = *i;
williamr@2:               if (marker.is_tagged(i_node)) 
williamr@2:                 continue;
williamr@2:               marker.mark_tagged(i_node);
williamr@2: 
williamr@2:               if (numbering.is_numbered(i_node)) {
williamr@2:                 adj_iter j, je;
williamr@2:                 for (tie(j, je) = adjacent_vertices(i_node, G); j != je; ++j) {
williamr@2:                   const vertex_t j_node = *j;
williamr@2:                   if (marker.is_not_tagged(j_node)) {
williamr@2:                     marker.mark_tagged(j_node);
williamr@2:                     deg += supernode_size[j_node];
williamr@2:                   }
williamr@2:                 }
williamr@2:               } else
williamr@2:                 deg += supernode_size[i_node];
williamr@2:             } // for adjacent vertices of u_node
williamr@2:             deg -= supernode_size[u_node];
williamr@2:             degree[u_node] = deg;
williamr@2:             degreelists[deg].push(u_node);
williamr@2:             // u_id has been pushed back into degreelists
williamr@2:             degree_lists_marker.unmark(u_node); 
williamr@2:             if (min_degree > deg)
williamr@2:               min_degree = deg;
williamr@2:             qxlist.pop();
williamr@2:           } // while (!qxlist.empty()) {
williamr@2: 
williamr@2:           marker.set_tag_as_multiple_tag();
williamr@2:           llist.pop();
williamr@2:         } // while (! llist.empty())
williamr@2: 
williamr@2:       } // update()
williamr@2: 
williamr@2: 
williamr@2:       void build_permutation(InversePermutationMap next,
williamr@2:                              PermutationMap prev) 
williamr@2:       {
williamr@2:         // collect the permutation info
williamr@2:         size_type i;
williamr@2:         for (i = 0; i < n; ++i) {
williamr@2:           diff_t size = supernode_size[get(index_vertex_map, i)];
williamr@2:           if ( size <= 0 ) {
williamr@2:             prev[i] = next[i];
williamr@2:             supernode_size[get(index_vertex_map, i)]
williamr@2:               = next[i] + 1;  // record the supernode info
williamr@2:           } else
williamr@2:             prev[i] = - next[i];
williamr@2:         }
williamr@2:         for (i = 1; i < n + 1; ++i) {
williamr@2:           if ( prev[i-1] > 0 )
williamr@2:             continue;
williamr@2:           diff_t parent = i;
williamr@2:           while ( prev[parent - 1] < 0 ) {
williamr@2:             parent = - prev[parent - 1];
williamr@2:           }
williamr@2: 
williamr@2:           diff_t root = parent;
williamr@2:           diff_t num = prev[root - 1] + 1;
williamr@2:           next[i-1] = - num;
williamr@2:           prev[root-1] = num;
williamr@2: 
williamr@2:           parent = i;
williamr@2:           diff_t next_node = - prev[parent - 1];
williamr@2:           while (next_node > 0) {
williamr@2:             prev[parent-1] = - root;
williamr@2:             parent = next_node;
williamr@2:             next_node = - prev[parent - 1];
williamr@2:           }
williamr@2:         }
williamr@2:         for (i = 0; i < n; i++) {
williamr@2:           diff_t num = - next[i] - 1;
williamr@2:           next[i] = num;
williamr@2:           prev[num] = i;
williamr@2:         }
williamr@2:       } // build_permutation()
williamr@2:     };
williamr@2: 
williamr@2:   } //namespace detail
williamr@2: 
williamr@2: 
williamr@2:   // MMD algorithm
williamr@2:   // 
williamr@2:   //The implementation presently includes the enhancements for mass
williamr@2:   //elimination, incomplete degree update, multiple elimination, and
williamr@2:   //external degree.
williamr@2:   //
williamr@2:   //Important Note: This implementation requires the BGL graph to be
williamr@2:   //directed.  Therefore, nonzero entry (i, j) in a symmetrical matrix
williamr@2:   //A coresponds to two directed edges (i->j and j->i).
williamr@2:   //
williamr@2:   //see Alan George and Joseph W. H. Liu, The Evolution of the Minimum
williamr@2:   //Degree Ordering Algorithm, SIAM Review, 31, 1989, Page 1-19
williamr@2:   template<class Graph, class DegreeMap, 
williamr@2:            class InversePermutationMap, 
williamr@2:            class PermutationMap,
williamr@2:            class SuperNodeMap, class VertexIndexMap>
williamr@2:   void minimum_degree_ordering
williamr@2:     (Graph& G, 
williamr@2:      DegreeMap degree, 
williamr@2:      InversePermutationMap inverse_perm, 
williamr@2:      PermutationMap perm, 
williamr@2:      SuperNodeMap supernode_size, 
williamr@2:      int delta, 
williamr@2:      VertexIndexMap vertex_index_map)
williamr@2:   {
williamr@2:     detail::mmd_impl<Graph,DegreeMap,InversePermutationMap,
williamr@2:       PermutationMap, SuperNodeMap, VertexIndexMap> 
williamr@2:       impl(G, num_vertices(G), delta, degree, inverse_perm, 
williamr@2:            perm, supernode_size, vertex_index_map);
williamr@2:     impl.do_mmd();
williamr@2:     impl.build_permutation(inverse_perm, perm);
williamr@2:   }
williamr@2: 
williamr@2: } // namespace boost
williamr@2: 
williamr@2: #endif // MINIMUM_DEGREE_ORDERING_HPP