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

#include <iostream>
#include <iomanip>
#include <vector>
#include <map>

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);
	}
};

//Dense Matrix for comparison
class Matrix : public std::vector<double>
{
  size_t m;
  size_t n;

public:
  Matrix() : m(0), n(0) {}
  Matrix(size_t m_, size_t n_) : std::vector<double>(m_*n_), m(m_), n(n_) {}
  Matrix(size_t m_, size_t n_, double val) : std::vector<double>(m_*n_, val), m(m_), n(n_) {}

  void resize(size_t m_, size_t n_)
    {m = m_; n = n_; std::vector<double>::resize(m_*n_);}
  void reserve(size_t m_, size_t n_)
    {std::vector<double>::reserve(m_*n_);}
  void clear()
    {m = n = 0; std::vector<double>::clear();}

  size_t NRows() const {return m;}
  size_t NCols() const {return n;}

  double& operator()(size_t i, size_t j)
  {
    return operator[](i + j*m);
  }
  const double& operator()(size_t i, size_t j) const
  {
    return operator[](i + j*m);
  }
	void swap(Matrix &y)
	{
		std::vector<double>::swap(y);
		std::swap(n, y.n);
		std::swap(m, y.m);
	}
	void clearMemory()
	{
		Matrix 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;
		}
}

// a = b
void convert2dense(Matrix& a, SparseMatrix& b)
{
	a.resize(b.NRows(), b.NCols());
	fill(a.begin(), a.end(), 0.);
	for (int i = 0; i < b.NRows(); ++i)
		for (map<int, double>::const_iterator j = b[i].begin(); j != b[i].end(); ++j)
			a(i, (*j).first) = (*j).second;
}

void printMatrix(Matrix &a, char *comment)
{
	cout << comment << endl;
	for (int i = 0; i < a.NRows(); ++i)
	{
		for (int j = 0; j < a.NCols(); ++j)
			cout << setw(10) << a(i, j);
		cout << endl;
	}
	cout << endl;
}

int main()
{
//Initialize the matrix	
	SparseMatrix sm(4);
	sm.m = 4;
	sm[0][0] = 1;
	sm[0][2] = 1;
	sm[1][0] = 0.5;
	sm[1][1] = 0.8;
	sm[2][1] = 0.6;
	sm[2][2] = 0.7;
	sm[3][2] = 0.3;
	sm[3][3] = 0.2;
//Convert to dense and print the matrix
  Matrix dm;
	convert2dense(dm, sm);
	printMatrix(dm, "Original Matrix A");
// multiplication of dense matrices to test results
	Matrix res(dm.NRows(), dm.NRows());
// res = dm*dm^t		
	for (int i = 0; i < dm.NRows(); i++)
		for (int j = 0; j < dm.NRows(); j++)
		{
			res(i, j) = 0.;
			for (int k = 0; k < dm.NCols(); k++)
				res(i, j) += dm(i, k)*dm(j, k);
		}
	printMatrix(res, "A*A^t by multiplication of dense matrices");
// multiplication of sparse matrices	
	SparseMatrix sres;
	multiply(sm, sres);
//convert to dense to print it	
	convert2dense(res, sres);
	printMatrix(res, "A*A^t by multiplication of sparse matrices, only half was computed");
}