/*
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 "matrices.h"
extern "C" {
#include "taucs.h"
}

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

// 2) converting matrix to the TAUCS representation		
// 2a) creating TAUCS matrix		
	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);				
	}
// 2b) converting SparseMatrix to the TAUCS CCS		
// 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.write("Converting matrix");
	if (argc < 3)
		return 0;
// 3) Reading RHS
	ifstream in2(argv[2]);
	if (!in2)
	{
		cout << "cannot open file " << argv[2] << endl;
		return 4;
	}
	Matrix rhs;
	rhs.read(in2);
	if (A->n != rhs.NRows())
	{
		cout << "x is not compatible with the matrix" << endl;
		return 4;
	}
// 4) solve the linear system 
//	A x = rhs
	Matrix x(rhs.NRows(), rhs.NCols());

// 4a) call taucs_linsolve 
  void* F = NULL;
  char* options[] = {"taucs.factor.LLT=true", NULL};
  void* opt_arg[] = { NULL };

	taucs_logfile("stdout");
  int i = taucs_linsolve(A, &F, rhs.NCols(), &*x.begin(), &*rhs.begin(), options, opt_arg);
// 4b) check if successful 
  if (i != TAUCS_SUCCESS)
  {
    cout << "Solution error." << endl;
    if (i==TAUCS_ERROR)
      cout << "Generic error." << endl;

    if (i==TAUCS_ERROR_NOMEM)
      cout << "NOMEM error." << endl;

    if (i==TAUCS_ERROR_BADARGS)
      cout << "BADARGS error." << endl;

    if (i==TAUCS_ERROR_MAXDEPTH)
      cout << "MAXDEPTH error." << endl;

    if (i==TAUCS_ERROR_INDEFINITE)
      cout << "NOT POSITIVE DEFINITE error." << endl;
  }
  else
  {
// 4c) compute the norm and write the solution out 
    cout << "Solution successful." << endl;
		vector<double> prod(A->n);
		for (int i = 0; i < rhs.NCols(); ++i)
		{
			taucs_dccs_times_vec(A, &*x.column(i), &*prod.begin());
			double sum = 0.;
			for (int j = 0; j < A->n; ++j)
				sum += pow(rhs(j, i) - prod[j], 2);
			cout << "norm Ax - b for RHS " << i + 1 << " is " << sqrt(sum) << endl;
		}		
		cout << "Solution is in " << argv[2] << ".solution" << endl;
		ofstream outf((string(argv[2]) + ".solution").c_str());
		x.write(outf);
  }

// 5) clean up
  taucs_linsolve(NULL, &F, 0, NULL, NULL, NULL, NULL);
	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;	
}