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