/*
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).
void terminate(const char* msg)
{
- std::cerr << msg << std::endl;
- abort();
+ std::cerr << msg << std::endl;
+ abort();
}
-//copy and modify from mtl harwell boeing stream
+// 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");
- }
+ 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;
}
- // 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);
+ ~harwell_boeing()
+ {
+ free(colptr);
+ free(rowind);
+ free(val);
+ }
- Vector degree(n, 0);
+ inline int nrows() const { return M; }
+
+ int cnt;
+ int col;
+ int* colptr;
+ bool isComplex;
+ int M;
+ int N;
+ int nonzeros;
+ int* rowind;
+ double* val;
+};
- 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);
+int main(int argc, char* argv[])
+{
+ using namespace std;
+ using namespace boost;
- if ( argc >= 3 ) {
- ifstream input(argv[2]);
- if ( input.fail() ) {
- cout << argv[3] << " is failed to open!. " << endl;
- return -1;
+ if (argc < 2)
+ {
+ cout << argv[0] << " HB file" << 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;
- }
+
+ 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;
}
- 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;
+ return 0;
}
projects / ceph.git / blobdiff