// Evgenii B. Rudnyi, http://Evgenii.Rudnyi.Ru

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

using namespace::std;

//The class keeps the matrix by rows
class SparseMatrix : public vector< map<int, double> >
{
	void init()
	{
		m = 0;
	}
public:	
	size_t m;
	
	SparseMatrix() {init();}
  explicit SparseMatrix(size_t n_) : vector< map<int, double> >(n_) {init();}
	void clear()
	{
		vector< map<int, double> >::clear();
		init();
	}

 	size_t NRows() const {return size();}
  size_t NCols() const {return m;}
	
	void swap(SparseMatrix &y)
	{
		vector< map<int, double> >::swap(y);
		std::swap(m, y.m);
	}
	void clearMemory()
	{
		SparseMatrix empty;
	  swap(empty);
	}
};

// C = A A'
// A is assumed to be unsymmetric
void multiply(SparseMatrix &A, SparseMatrix &C)
{
// C_ij = Sum_k A_ik*B_kj
// B = A'  b_kj = a_jk
// C_ij = Sum_k A_ik*A_jk

// Let us compute only the upper half of C
  C.clear();
	int dim = A.NRows();
	C.resize(dim);
	C.m = dim;
	for (int i = 0; i < dim; ++i)
		for (int j = i; j < dim; ++j)
		{
			double sum = 0.;
			if (i == j)
			{
				for (map<int, double>::const_iterator k = A[i].begin(); k != A[i].end(); ++k)
					sum += (*k).second*(*k).second;
			}
			else
			{
				for (map<int, double>::const_iterator k = A[i].begin(); k != A[i].end(); ++k)
				{
					map<int, double>::const_iterator kk = A[j].find((*k).first);
					if (kk != A[j].end())
						sum += (*k).second*(*kk).second;
				}
			}
			if (sum != 0)
				C[i][j] = sum;
		}
}

// Symmetric matrix C above is actually equivalent to CCS required by TAUCS
// Attention! This is working only for symmetric matrices
taucs_ccs_matrix* convert2taucs(SparseMatrix &C) 
{
	taucs_ccs_matrix *A;
	int dim = C.size(); 
  int nnz = 0;
	for (int i = 0; i < dim; ++i)
	{
		int num = C[i].size();
			nnz += num;
	}
	A = taucs_dccs_create(dim, dim, nnz);
	A->flags |= TAUCS_SYMMETRIC | TAUCS_LOWER;
	int pos = 0;
	for (int i = 0; i < dim; ++i)
	{
		A->colptr[i] = pos;
		for (map<int, double>::const_iterator j = C[i].begin(); j != C[i].end(); ++j)
		{
			A->rowind[pos] = (*j).first;
			A->values.d[pos] = (*j).second;	
			++pos;
		}
	}
	A->colptr[dim] = nnz;
	return A;
}