ceph/src/boost/libs/graph/example/minimum_degree_ordering.cpp

   1 //-*-c++-*-
   2 //=======================================================================
   3 // Copyright 1997-2001 University of Notre Dame.
   4 // Authors: Lie-Quan Lee
   5 //
   6 // Distributed under the Boost Software License, Version 1.0. (See
   7 // accompanying file LICENSE_1_0.txt or copy at
   8 // http://www.boost.org/LICENSE_1_0.txt)
   9 //=======================================================================
  10
  11 /*
  12   This file is to demo how to use minimum_degree_ordering algorithm.
  13
  14   Important Note: This implementation requires the BGL graph to be
  15   directed.  Therefore, nonzero entry (i, j) in a symmetrical matrix
  16   A coresponds to two directed edges (i->j and j->i).
  17
  18   The bcsstk01.rsa is an example graph in Harwell-Boeing format,
  19   and bcsstk01 is the ordering produced by Liu's MMD implementation.
  20   Link this file with iohb.c to get the harwell-boeing I/O functions.
  21   To run this example, type:
  22
  23   ./minimum_degree_ordering bcsstk01.rsa bcsstk01
  24
  25 */
  26
  27 #include <boost/config.hpp>
  28 #include <fstream>
  29 #include <iostream>
  30 #include "boost/graph/adjacency_list.hpp"
  31 #include "boost/graph/graph_utility.hpp"
  32 #include "boost/graph/minimum_degree_ordering.hpp"
  33
  34 void terminate(const char* msg)
  35 {
  36     std::cerr << msg << std::endl;
  37     abort();
  38 }
  39
  40 // copy and modify from mtl harwell boeing stream
  41 struct harwell_boeing
  42 {
  43     harwell_boeing(char* filename)
  44     {
  45         int Nrhs;
  46         char* Type;
  47         Type = new char[4];
  48         isComplex = false;
  49         // Never called:
  50         // readHB_info(filename, &M, &N, &nonzeros, &Type, &Nrhs);
  51         colptr = (int*)malloc((N + 1) * sizeof(int));
  52         if (colptr == NULL)
  53             terminate("Insufficient memory for colptr.\n");
  54         rowind = (int*)malloc(nonzeros * sizeof(int));
  55         if (rowind == NULL)
  56             terminate("Insufficient memory for rowind.\n");
  57
  58         if (Type[0] == 'C')
  59         {
  60             isComplex = true;
  61             val = (double*)malloc(nonzeros * sizeof(double) * 2);
  62             if (val == NULL)
  63                 terminate("Insufficient memory for val.\n");
  64         }
  65         else
  66         {
  67             if (Type[0] != 'P')
  68             {
  69                 val = (double*)malloc(nonzeros * sizeof(double));
  70                 if (val == NULL)
  71                     terminate("Insufficient memory for val.\n");
  72             }
  73         }
  74         // Never called:
  75         // readHB_mat_double(filename, colptr, rowind, val);
  76
  77         cnt = 0;
  78         col = 0;
  79         delete[] Type;
  80     }
  81
  82     ~harwell_boeing()
  83     {
  84         free(colptr);
  85         free(rowind);
  86         free(val);
  87     }
  88
  89     inline int nrows() const { return M; }
  90
  91     int cnt;
  92     int col;
  93     int* colptr;
  94     bool isComplex;
  95     int M;
  96     int N;
  97     int nonzeros;
  98     int* rowind;
  99     double* val;
 100 };
 101
 102 int main(int argc, char* argv[])
 103 {
 104     using namespace std;
 105     using namespace boost;
 106
 107     if (argc < 2)
 108     {
 109         cout << argv[0] << " HB file" << endl;
 110         return -1;
 111     }
 112
 113     int delta = 0;
 114
 115     if (argc >= 4)
 116         delta = atoi(argv[3]);
 117
 118     harwell_boeing hbs(argv[1]);
 119
 120     // must be BGL directed graph now
 121     typedef adjacency_list< vecS, vecS, directedS > Graph;
 122
 123     int n = hbs.nrows();
 124
 125     cout << "n is " << n << endl;
 126
 127     Graph G(n);
 128
 129     int num_edge = 0;
 130
 131     for (int i = 0; i < n; ++i)
 132         for (int j = hbs.colptr[i]; j < hbs.colptr[i + 1]; ++j)
 133             if ((hbs.rowind[j - 1] - 1) > i)
 134             {
 135                 add_edge(hbs.rowind[j - 1] - 1, i, G);
 136                 add_edge(i, hbs.rowind[j - 1] - 1, G);
 137                 num_edge++;
 138             }
 139
 140     cout << "number of off-diagnal elements: " << num_edge << endl;
 141
 142     typedef std::vector< int > Vector;
 143
 144     Vector inverse_perm(n, 0);
 145     Vector perm(n, 0);
 146
 147     Vector supernode_sizes(n, 1); // init has to be 1
 148
 149     boost::property_map< Graph, vertex_index_t >::type id
 150         = get(vertex_index, G);
 151
 152     Vector degree(n, 0);
 153
 154     minimum_degree_ordering(G,
 155         make_iterator_property_map(&degree[0], id, degree[0]), &inverse_perm[0],
 156         &perm[0],
 157         make_iterator_property_map(&supernode_sizes[0], id, supernode_sizes[0]),
 158         delta, id);
 159
 160     if (argc >= 3)
 161     {
 162         ifstream input(argv[2]);
 163         if (input.fail())
 164         {
 165             cout << argv[3] << " is failed to open!. " << endl;
 166             return -1;
 167         }
 168         int comp;
 169         bool is_correct = true;
 170         int i;
 171         for (i = 0; i < n; i++)
 172         {
 173             input >> comp;
 174             if (comp != inverse_perm[i] + 1)
 175             {
 176                 cout << "at i= " << i << ": " << comp
 177                      << " ***is NOT EQUAL to*** " << inverse_perm[i] + 1
 178                      << endl;
 179                 is_correct = false;
 180             }
 181         }
 182         for (i = 0; i < n; i++)
 183         {
 184             input >> comp;
 185             if (comp != perm[i] + 1)
 186             {
 187                 cout << "at i= " << i << ": " << comp
 188                      << " ***is NOT EQUAL to*** " << perm[i] + 1 << endl;
 189                 is_correct = false;
 190             }
 191         }
 192         if (is_correct)
 193             cout << "Permutation and inverse permutation are correct. " << endl;
 194         else
 195             cout << "WARNING -- Permutation or inverse permutation is not the "
 196                  << "same ones generated by Liu's " << endl;
 197     }
 198     return 0;
 199 }