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