// The code from Alejandro modifed by EBR to use low-level TAUCS functions

#include <iostream>
#include <vector>
extern "C" {
#include <taucs.h>
}

using namespace::std;

int main(){
	vector<double> an(10);
	vector<int> jn(10);
	vector<int> ia(10);
	vector<double> f(10); // right-hand size vector object

	// create CCS matrix structure using vector class
	an[0] = 1.0;
	an[1] = 0.5;
	an[2] = 1.0;
	an[3] = 0.5;
	an[4] = 1.0;
	an[5] = 0.5;
	an[6] = 1.0;

	jn[0] = 0;
	jn[1] = 1;
	jn[2] = 1;
	jn[3] = 2;
	jn[4] = 2;
	jn[5] = 3;
	jn[6] = 3;

	ia[0] = 0;
	ia[1] = 2;
	ia[2] = 4;
	ia[3] = 6;
	ia[4] = 7;

	// create right-hand size vector object
	f[0] = 1.0;
	f[1] = 2.0;
	f[2] = 3.0;
	f[3] = 4.0;

	// resize vectors.
	an.resize(7);
	jn.resize(7);
	ia.resize(5);
	f.resize(4);
	int dim = 4;

	// create TAUCS matrix from vector objects an, jn and ia
	taucs_ccs_matrix  A; // a matrix to solve Ax=b in CCS format
	A.n = dim;
	A.m = dim;
	A.flags = (TAUCS_DOUBLE | TAUCS_SYMMETRIC | TAUCS_LOWER);
	A.colptr = &ia[0];
	A.rowind = &jn[0];
	A.values.d = &an[0];

	// we need solution and RHS twice: original and reordered
	// create TAUCS right-hand size
	taucs_double* b = &f[0]; // right hand side vector to solve Ax=b
	vector <double> bodv(dim);
	taucs_double* bod = &*bodv.begin(); 

	// allocate TAUCS solution vector
	vector <double> xv(dim);
	taucs_double* x = &*xv.begin(); 
	vector <double> xodv(dim);
	taucs_double* xod = &*xodv.begin(); 

	//Using TAUCS low-level routines
  int*         perm;
  int*         invperm;
	taucs_ccs_matrix*  Aod;
	void* F;

	taucs_logfile("stdout");
	// 1) Reordering
	taucs_ccs_order(&A, &perm, &invperm, "metis");
  Aod = taucs_ccs_permute_symmetrically(&A, perm, invperm);
  taucs_vec_permute(dim, TAUCS_DOUBLE, b, bod, perm);

	// 2) Factoring
  F = taucs_ccs_factor_llt_mf(Aod);	

	// 3) Back substitution and reodering the solution back
	taucs_supernodal_solve_llt(F, xod, bod);	
  taucs_vec_ipermute(dim, TAUCS_DOUBLE, xod, x, perm);

	// 4) Clean-up
	taucs_supernodal_factor_free(F);
 	taucs_ccs_free(Aod);			

	cout << "Solution" << endl;
  for (unsigned j = 0; j < xv.size(); j++)
		cout << x[j] << endl;
}