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