]>
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. | |
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" | |
92f5a8d4 TL |
33 | |
34 | void terminate(const char* msg) | |
35 | { | |
36 | std::cerr << msg << std::endl; | |
37 | abort(); | |
38 | } | |
7c673cae FG |
39 | |
40 | //copy and modify from mtl harwell boeing stream | |
41 | struct harwell_boeing | |
42 | { | |
43 | harwell_boeing(char* filename) { | |
44 | int Nrhs; | |
45 | char* Type; | |
46 | Type = new char[4]; | |
47 | isComplex = false; | |
92f5a8d4 TL |
48 | // Never called: |
49 | //readHB_info(filename, &M, &N, &nonzeros, &Type, &Nrhs); | |
7c673cae | 50 | colptr = (int *)malloc((N+1)*sizeof(int)); |
92f5a8d4 | 51 | if ( colptr == NULL ) terminate("Insufficient memory for colptr.\n"); |
7c673cae | 52 | rowind = (int *)malloc(nonzeros*sizeof(int)); |
92f5a8d4 | 53 | if ( rowind == NULL ) terminate("Insufficient memory for rowind.\n"); |
7c673cae FG |
54 | |
55 | if ( Type[0] == 'C' ) { | |
56 | isComplex = true; | |
57 | val = (double *)malloc(nonzeros*sizeof(double)*2); | |
92f5a8d4 | 58 | if ( val == NULL ) terminate("Insufficient memory for val.\n"); |
7c673cae FG |
59 | |
60 | } else { | |
61 | if ( Type[0] != 'P' ) { | |
62 | val = (double *)malloc(nonzeros*sizeof(double)); | |
92f5a8d4 | 63 | if ( val == NULL ) terminate("Insufficient memory for val.\n"); |
7c673cae FG |
64 | } |
65 | } | |
92f5a8d4 TL |
66 | // Never called: |
67 | //readHB_mat_double(filename, colptr, rowind, val); | |
7c673cae FG |
68 | |
69 | cnt = 0; | |
70 | col = 0; | |
71 | delete [] Type; | |
72 | } | |
73 | ||
74 | ~harwell_boeing() { | |
75 | free(colptr); | |
76 | free(rowind); | |
77 | free(val); | |
78 | } | |
79 | ||
80 | inline int nrows() const { return M; } | |
81 | ||
82 | int cnt; | |
83 | int col; | |
84 | int* colptr; | |
85 | bool isComplex; | |
86 | int M; | |
87 | int N; | |
88 | int nonzeros; | |
89 | int* rowind; | |
90 | double* val; | |
91 | }; | |
92 | ||
93 | int main(int argc, char* argv[]) | |
94 | { | |
95 | using namespace std; | |
96 | using namespace boost; | |
97 | ||
98 | if (argc < 2) { | |
99 | cout << argv[0] << " HB file" << endl; | |
100 | return -1; | |
101 | } | |
102 | ||
103 | int delta = 0; | |
104 | ||
105 | if ( argc >= 4 ) | |
106 | delta = atoi(argv[3]); | |
107 | ||
108 | harwell_boeing hbs(argv[1]); | |
109 | ||
110 | //must be BGL directed graph now | |
111 | typedef adjacency_list<vecS, vecS, directedS> Graph; | |
112 | ||
113 | int n = hbs.nrows(); | |
114 | ||
115 | cout << "n is " << n << endl; | |
116 | ||
117 | Graph G(n); | |
118 | ||
119 | int num_edge = 0; | |
120 | ||
121 | for (int i = 0; i < n; ++i) | |
122 | for (int j = hbs.colptr[i]; j < hbs.colptr[i+1]; ++j) | |
123 | if ( (hbs.rowind[j - 1] - 1 ) > i ) { | |
124 | add_edge(hbs.rowind[j - 1] - 1, i, G); | |
125 | add_edge(i, hbs.rowind[j - 1] - 1, G); | |
126 | num_edge++; | |
127 | } | |
128 | ||
129 | cout << "number of off-diagnal elements: " << num_edge << endl; | |
130 | ||
131 | typedef std::vector<int> Vector; | |
132 | ||
133 | Vector inverse_perm(n, 0); | |
134 | Vector perm(n, 0); | |
135 | ||
136 | Vector supernode_sizes(n, 1); // init has to be 1 | |
137 | ||
138 | boost::property_map<Graph, vertex_index_t>::type | |
139 | id = get(vertex_index, G); | |
140 | ||
141 | Vector degree(n, 0); | |
142 | ||
143 | minimum_degree_ordering | |
144 | (G, | |
145 | make_iterator_property_map(°ree[0], id, degree[0]), | |
146 | &inverse_perm[0], | |
147 | &perm[0], | |
148 | make_iterator_property_map(&supernode_sizes[0], id, supernode_sizes[0]), | |
149 | delta, id); | |
150 | ||
151 | if ( argc >= 3 ) { | |
152 | ifstream input(argv[2]); | |
153 | if ( input.fail() ) { | |
154 | cout << argv[3] << " is failed to open!. " << endl; | |
155 | return -1; | |
156 | } | |
157 | int comp; | |
158 | bool is_correct = true; | |
159 | int i; | |
160 | for ( i=0; i<n; i++ ) { | |
161 | input >> comp; | |
162 | if ( comp != inverse_perm[i]+1 ) { | |
163 | cout << "at i= " << i << ": " << comp | |
164 | << " ***is NOT EQUAL to*** " << inverse_perm[i]+1 << endl; | |
165 | is_correct = false; | |
166 | } | |
167 | } | |
168 | for ( i=0; i<n; i++ ) { | |
169 | input >> comp; | |
170 | if ( comp != perm[i]+1 ) { | |
171 | cout << "at i= " << i << ": " << comp | |
172 | << " ***is NOT EQUAL to*** " << perm[i]+1 << endl; | |
173 | is_correct = false; | |
174 | } | |
175 | } | |
176 | if ( is_correct ) | |
177 | cout << "Permutation and inverse permutation are correct. "<< endl; | |
178 | else | |
179 | cout << "WARNING -- Permutation or inverse permutation is not the " | |
180 | << "same ones generated by Liu's " << endl; | |
181 | ||
182 | } | |
183 | return 0; | |
184 | } |