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
terminate(const char * msg)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 {
harwell_boeingharwell_boeing43 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
~harwell_boeingharwell_boeing82 ~harwell_boeing()
83 {
84 free(colptr);
85 free(rowind);
86 free(val);
87 }
88
nrowsharwell_boeing89 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
main(int argc,char * argv[])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(°ree[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 }
200