//-*-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"
#include "iohb.h"

//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;
    readHB_info(filename, &M, &N, &nonzeros, &Type, &Nrhs);
    colptr = (int *)malloc((N+1)*sizeof(int));
    if ( colptr == NULL ) IOHBTerminate("Insufficient memory for colptr.\n");
    rowind = (int *)malloc(nonzeros*sizeof(int));
    if ( rowind == NULL ) IOHBTerminate("Insufficient memory for rowind.\n");

    if ( Type[0] == 'C' ) {
      isComplex = true;
      val = (double *)malloc(nonzeros*sizeof(double)*2);
      if ( val == NULL ) IOHBTerminate("Insufficient memory for val.\n");

    } else {
      if ( Type[0] != 'P' ) {   
        val = (double *)malloc(nonzeros*sizeof(double));
        if ( val == NULL ) IOHBTerminate("Insufficient memory for val.\n");
      }
    }

    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]);
  
  typedef double Type;
 
  harwell_boeing hbs(argv[1]);

  //must be BGL directed graph now
  typedef adjacency_list<vecS, vecS, directedS>  Graph;
  typedef graph_traits<Graph>::vertex_descriptor Vertex;

  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(&degree[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;
}