/*
Copyright (C) 2006-2009 Evgenii Rudnyi, http://Evgenii.Rudnyi.Ru/

The program reads the matrix and RHS in the Matrix Market format then run TAUCS.

This software is a copyrighted work licensed under the terms, described
in the file "FREE_LICENSE". 
*/

#include <fstream>
#include <iomanip>
#include "auxiliary.h"
extern "C" {
#include "taucs.h"
}

int main(int argc, char **argv)
{	
	if (argc < 2)
	{
		cout << "run_taucs matrix rhs" << endl;
  	return 1;	
	}
	try
	{
  SparseMatrix mat;
	ifstream in(argv[1]);
	if (!in)
	{
		cout << "cannot open file " << argv[1] << endl;
		return 2;
	}
	Timing t0;
 	mat.read(in);
  t0.print("Reading matrix");

//converting matrix to the TAUCS representation		
	Timing r1;
	taucs_ccs_matrix *A;
	A = taucs_dccs_create(mat.NRows(), mat.NCols(), mat.nnz);
	A->flags = TAUCS_DOUBLE;
	if (mat.IsSymmetric)
		A->flags |= TAUCS_SYMMETRIC | TAUCS_LOWER;
	else
	{
// SparseMatrix is by rows - it is necessary to change it by columns		
		SparseMatrix tmp(mat.NCols());
		for (int i = 0; i < mat.NRows(); ++i)
			for (map<int, double>::const_iterator j = mat[i].begin(); j != mat[i].end(); ++j)
				tmp[(*j).first][i] = (*j).second;
		tmp.nnz = mat.nnz;
		tmp.m = mat.NRows();
		mat.swap(tmp);				
	}
// now mat.NCols() means NRows and vice virsa	
	int pos = 0;
	for (int i = 0; i < mat.NRows(); ++i)
	{
		A->colptr[i] = pos;
		for (map<int, double>::const_iterator j = mat[i].begin(); j != mat[i].end(); ++j)
		{
			A->rowind[pos] = (*j).first;
			A->values.d[pos] = (*j).second;	
			++pos;
		}
	}
	A->colptr[mat.NRows()] = mat.nnz;
	mat.clearMemory();
	r1.print("Converting matrix");
	if (argc < 3)
		return 0;
// Preparing RHS
	ifstream in2(argv[2]);
	if (!in2)
	{
		cout << "cannot open file " << argv[2] << endl;
		return 4;
	}
	Matrix x;
	x.read(in2);
	if (A->n != x.NRows())
	{
		cout << "x is not compatible with the matrix" << endl;
		return 4;
	}
	Matrix res(A->m, 1);
//	res = A x
	taucs_dccs_times_vec(A, &*x.begin(), &*res.begin());
  cout << "Solution is in " << argv[2] << ".product" << endl;
	ofstream outf((string(argv[2]) + ".product").c_str());
	res.write(outf);
	taucs_dccs_free(A);
	}
	catch (BaseError &t)
	{
		cout << "Error is encountered:" << endl;
		cout << t << endl;
		return 5;
	}
	catch (...)
	{
		cout << "Not known error is encountered." << endl;
		return 6;
	}
  return 0;	
}