]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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. | |
f67539c2 | 13 | |
7c673cae FG |
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" | |
92f5a8d4 TL |
33 | |
34 | void terminate(const char* msg) | |
35 | { | |
f67539c2 TL |
36 | std::cerr << msg << std::endl; |
37 | abort(); | |
92f5a8d4 | 38 | } |
7c673cae | 39 | |
f67539c2 | 40 | // copy and modify from mtl harwell boeing stream |
7c673cae FG |
41 | struct harwell_boeing |
42 | { | |
f67539c2 TL |
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; | |
7c673cae | 80 | } |
7c673cae | 81 | |
f67539c2 TL |
82 | ~harwell_boeing() |
83 | { | |
84 | free(colptr); | |
85 | free(rowind); | |
86 | free(val); | |
87 | } | |
7c673cae | 88 | |
f67539c2 TL |
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 | }; | |
7c673cae | 101 | |
f67539c2 TL |
102 | int main(int argc, char* argv[]) |
103 | { | |
104 | using namespace std; | |
105 | using namespace boost; | |
7c673cae | 106 | |
f67539c2 TL |
107 | if (argc < 2) |
108 | { | |
109 | cout << argv[0] << " HB file" << endl; | |
110 | return -1; | |
7c673cae | 111 | } |
f67539c2 TL |
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; | |
7c673cae | 197 | } |
f67539c2 | 198 | return 0; |
7c673cae | 199 | } |