#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include <algorithm>
#include <limits.h>
#include <assert.h>




class CSRMatrix
{
	std::vector<int> ia, ja;
	std::vector<double> a;
public:
	CSRMatrix() {}
	CSRMatrix(const std::vector<int> & ia, const std::vector<int> & ja, const std::vector<double> & a)
	: ia(ia), ja(ja), a(a) {}
	CSRMatrix(const CSRMatrix & b) : ia(b.ia), ja(b.ja), a(b.a) {}
	CSRMatrix & operator=(CSRMatrix const & b) { ia = b.ia; ja = b.ja; a = b.a; return *this;}
	void Load(std::string file)
	{
		int N,M,L;
		int row, col, last_row;
		double val;
		std::ifstream input(file.c_str());
		ia.clear();
		ja.clear();
		a.clear();
		while( !input.eof() )
		{
			int c = input.peek();
			if( isspace(c) ) continue;
			if( c == '%' ) 
			{
				input.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
				continue;
			}
			input >> N >> M >> L;
			assert(N == M);
			ia.resize(N+1);
			ja.reserve(L);
			a.reserve(L);
			break;
		}
		assert( !ia.empty() );
		last_row = 0;
		ia[0] = 0;
		while( !(input >> row >> col >> val).eof() )
		{
			assert(row == last_row || row == last_row+1);
			if( row != last_row )
			{
				ia[row] = ia[last_row];
				last_row = row;
			}
			ia[row]++;
			ja.push_back(col-1);
			a.push_back(val);
		}
		input.close();
	}
	void Save(std::string file) const
	{
		std::ofstream output(file.c_str());
		output << "%%MatrixMarket matrix coordinate real general" << std::endl;
		output << Size() << " " << Size() << " " << ja.size() << std::endl;
		output << std::scientific;
	       	output.precision(16);
		for(int i = 0; i < Size(); ++i)
		{
			for(int j = ia[i]; j < ia[i+1]; ++j)
				output << i+1 << " " << ja[j]+1 << " " << a[j] << std::endl;
		}
		output.close();
	}
	int    Size() const {return (int)ia.size()-1;}
	int    RowSize(int i) const {return ia[i+1]-ia[i];}
	void   SortColumns()
	{
		for(int k = 0; k < Size(); ++k)
		{
			int i, j, n = ia[k+1] - ia[k], s = ia[k];  
			for (i = 0; i < n-1; i++)
			{
				for (j = 0; j < n-i-1; j++)
				{
					if (ja[s+j] > ja[s+j+1])
					{
						std::swap(ja[s+j],ja[s+j+1]);
						std::swap( a[s+j], a[s+j+1]);
					}
				}
			}
		}
	}
	friend class CSCTraversal;
	class CSCTraversal
	{
		std::vector<int> ia0;
		std::vector<int> first;
		std::vector<int> next;
		const CSRMatrix & A;
		int col, Cols;
		const int EOL = INT_MAX-1;
	public:
		CSCTraversal(const CSRMatrix & A, int Cols = 0) : A(A), Cols(Cols) 
		{
			ia0 = A.ia;
			if( Cols == 0 )
			{
				for(int i = 0; i < (int)A.ja.size(); ++i) 
					Cols = std::max(Cols,A.ja[i]+1);
			}
			first.resize(Cols,EOL);
			next.resize(A.Size(),INT_MAX);
			int i,k,j;
			for(k = A.Size(); k > 0; --k)
			{
				i = k-1;
				j = A.ja[ia0[i]];
				next[i] = first[j];
				first[j] = i;
			}
			col = 0;
		}
		int GetFirst() const {return first[col];}
		int GetNext(int i) const {return next[i];}
		int GetPosition(int i) const {return ia0[i];}
		int Last() const {return EOL;}
		void NextColumn()
		{
			assert(col < Cols);
			int k, l, j, cur, nxt;
			k = first[col++];
			while(k != EOL)
			{
				l = next[k];
				ia0[k]++;
				if( A.ia[k+1] - ia0[k] > 0 )
				{
					j = A.ja[ia0[k]];
					//ordered insert to linked-list
					if( first[j] > k )
					{
						next[k] = first[j];
						first[j] = k;
					}
					else
					{
						cur = nxt = first[j];
						while( nxt < k )
						{
							cur = nxt;
							nxt = next[cur];
						}
						next[cur] = k;
						next[k] = nxt;
					}
				}
				k = l;
			}
		}
	};
	CSRMatrix Transpose() const
	{
		int Cols = 0;
		for(int i = 0; i < (int)ja.size(); ++i) 
			Cols = std::max(Cols,ja[i]+1);
		std::vector< int > iat(Cols+1), jat;
		std::vector<double> at;
		CSCTraversal At(*this,Cols);
		//iterate over columns
		iat[0] = 0;
		for(int i = 0; i < Cols; ++i)
		{
			//fill row of transposed matrix
			int k = At.GetFirst();
			while(k != At.Last())
			{
				jat.push_back(k);
				at.push_back(a[At.GetPosition(k)]);
				k = At.GetNext(k);
			}
			iat[i+1] = jat.size();
			//advance to next column
			At.NextColumn();
		}
		return CSRMatrix(iat,jat,at);
	}
	bool Compare(CSRMatrix const & b, double eps) const
	{
		bool ret = true;
		double diftol = 0;
		if( Size() == b.Size() )
		{
			for(int i = 0; i < Size(); ++i)
			{
				if( RowSize(i) == b.RowSize(i) )
				{
					for(int j = ia[i]; j < ia[i+1]; ++j)
					{
						if( ja[j] == b.ja[j] )
						{
							double div = std::max(std::numeric_limits<double>::min(),fabs(a[j])+fabs(b.a[j]));
							double dif = fabs(a[j] - b.a[j]);
							if( dif > eps*div ) 
							{
								std::cout << "at ("<<i<<","<<ja[j]<<"), values are " << a[j] << " and " << b.a[j] << " relative difference " << dif/div << std::endl;
								ret = false;
							}
							diftol = std::max(diftol,dif/div);
						}
						else 
						{
							std::cout << "at row " << i << ", columns are " << ja[j] << " and " << b.ja[j] << std::endl;
							ret = false;
						}
					}
				}
				else 
				{
					std::cout << "at row " <<i << ", row sizes are " << RowSize(i) << " and " << b.RowSize(i) << std::endl;
					ret = false;
				}
			}
			if( ret ) 
				std::cout << "Matrices are same" << std::endl;
			else
				std::cout << "Relative difference tolerance " << diftol << std::endl;
		}
		else 
		{
			std::cout << "Matrix sizes " << Size() << " and " << b.Size() << std::endl;
			ret = false;
		}
		return ret;
	}
};


int main(int argc, char ** argv)
{
	if( argc < 2 )
		std::cout << argv[0] << " matrix.mtx [tol=1.0e-5]" << std::endl;
	else
	{
		double tol = 1.0e-5;
		CSRMatrix A;
		A.Load(std::string(argv[1]));
		A.SortColumns(); //required for comparison
		if( argc > 2 ) tol = atof(argv[2]);
		if( A.Compare(A.Transpose(),tol) ) return 0;
	}
	return 1;
}