123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184 |
- //-*-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"
- void terminate(const char* msg)
- {
- std::cerr << msg << std::endl;
- abort();
- }
- //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;
- // Never called:
- //readHB_info(filename, &M, &N, &nonzeros, &Type, &Nrhs);
- colptr = (int *)malloc((N+1)*sizeof(int));
- if ( colptr == NULL ) terminate("Insufficient memory for colptr.\n");
- rowind = (int *)malloc(nonzeros*sizeof(int));
- if ( rowind == NULL ) terminate("Insufficient memory for rowind.\n");
- if ( Type[0] == 'C' ) {
- isComplex = true;
- val = (double *)malloc(nonzeros*sizeof(double)*2);
- if ( val == NULL ) terminate("Insufficient memory for val.\n");
- } else {
- if ( Type[0] != 'P' ) {
- val = (double *)malloc(nonzeros*sizeof(double));
- if ( val == NULL ) terminate("Insufficient memory for val.\n");
- }
- }
- // Never called:
- //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]);
-
- harwell_boeing hbs(argv[1]);
- //must be BGL directed graph now
- typedef adjacency_list<vecS, vecS, directedS> Graph;
- 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;
- }
|