#include <cstdlib>
#include <iostream>
#include <boost/numeric/bindings/ublas/matrix.hpp>
#include <boost/numeric/bindings/ublas/vector.hpp>
#include <boost/numeric/bindings/ublas/matrix_sparse.hpp>
#include <boost/numeric/bindings/umfpack/umfpack.hpp>
#include <boost/numeric/ublas/io.hpp>

using namespace std;
namespace ublas = boost::numeric::ublas;
namespace umf = boost::numeric::bindings::umfpack;

int main()
{
	int n = 5;
	ublas::compressed_matrix<double, ublas::column_major, 0, ublas::unbounded_array<int>, ublas::unbounded_array<double> > A (n,n,12); 
	ublas::vector<double> B (n), X (n);

	A(0,0) = 2.; A(0,1) = 3; 
	A(1,0) = 3.; A(1,2) = 4.; A(1,4) = 6;
	A(2,1) = -1.; A(2,2) = -3.; A(2,3) = 2.;
	A(3,2) = 1.;
	A(4,1) = 4.; A(4,2) = 2.; A(4,4) = 1.; 

	B(0) = 8.; B(1) = 45.; B(2) = -3.; B(3) = 3.; B(4) = 19.; 

	umf::symbolic_type<double> Symbolic;
    umf::numeric_type<double> Numeric;
    umf::control_type<double> Control;
    umf::info_type<double> Info;
    
    int do_recip, status;
    
    Control[UMFPACK_PRL] = 6;

	umf::symbolic (A, Symbolic); 
	umf::numeric (A, Symbolic, Numeric); 
	
	// Testing bindings
    umf::report_numeric(Numeric, Control);
    int lnz, unz, n_row, n_col, nz_udiag;
    umf::getLunz(&lnz, &unz, &n_row, &n_col, &nz_udiag, Numeric); // This works!
    int lnz1 = max(lnz, 1);
    int unz1 = max(unz, 1);
    ublas::vector<int> P(n_row), Q(n_col);
    ublas::vector<double> Dx(min(n_row, n_col)), Rs(n_row);
    ublas::compressed_matrix<double, ublas::column_major, 0, ublas::unbounded_array<int>, ublas::unbounded_array<double> > Lx(n_row, min(n_row, n_col), lnz1), Ux(min(n_row, n_col), n_col, unz1);
    umf::getNumeric(Lx, Ux, P, Q, Dx, &do_recip, Rs, Numeric);
    cout << "Lx = " << Lx << endl; // Shouldn't be all zeros! :(
    
    // From UMFPACK's demo
    int* Lp = (int*) malloc((n+1) * sizeof (int));
    int* Lj = (int*) malloc(lnz1 * sizeof (int));
    double* Lxx = (double*) malloc(lnz1 * sizeof (double));
    int* Up = (int*) malloc((n+1) * sizeof (int));
    int* Ui = (int*) malloc(unz1 * sizeof (int));
    double* Uxx = (double*) malloc(unz1 * sizeof (double));
    int* PP = (int*) malloc(n * sizeof (int));
    int* QQ = (int*) malloc(n * sizeof (int));
    double* Dxx = (double*) malloc(min(n_row, n_col) * sizeof (double));
    double* Rss  = (double *) malloc (n * sizeof (double));
    if (!Lp || !Lj || !Lxx || !Up || !Ui || !Uxx || !PP || !QQ || !Rss)
    {
		cerr << "out of memory" << endl;
    }
    status = umfpack_di_get_numeric(Lp, Lj, Lxx, Up, Ui, Uxx, PP, QQ, Dxx, &do_recip, Rss, Numeric.ptr);
    if (status < 0)
    {
		cerr << "umfpack_di_get_numeric failed" << endl;
    }
    cout << "Lxx = ";
    for ( int i = 0; i < lnz1; i++ )
    	cout << Lxx[i] << " "; // Lxx is correct! But Lx isn't!
    cout << endl;
    
    umf::solve (A, X, B, Numeric);

	cout << "X = " << X << endl;  // output: [5](1,2,3,4,5)
	
	free(Lp);
	free(Lj);
	free(Lxx);
	free(Up);
	free(Ui);
	free(Uxx);
	free(PP);
	free(QQ);
	free(Dxx);
	free(Rss);

	return 0;
}