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

void LoadVector(std::string file, std::vector<double> & v)
{
	int N, row;
	double val;
	std::ifstream input(file.c_str());
	v.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;
		v.resize(N);
		break;
	}
	assert( !v.empty() );
	row = 0;
	while( !(input >> val).eof() )
		v[row++] = val;
	input.close();
}

class CSRMatrix
{
	std::vector<int> ia, ja;
	std::vector<double> a;
public:
	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();
	}
	// y = alpha*A*x + beta*y
	void Multiply(double alpha, const std::vector<double> & x, double beta, std::vector<double> & y) const
	{
		assert(x.size() == y.size() && x.size() == Size());
		for(int i = 0; i < Size(); ++i)
		{
			y[i] *= beta;
			for(int j = ia[i]; j < ia[i+1]; ++j)
				y[i] += alpha*a[j]*x[ja[j]];
		}
	}
	int Size() const {return (int)ia.size()-1;}
};

int main(int argc, char ** argv)
{
	if( argc < 4 )
		std::cout << argv[0] << " matrix.mtx rhs.mtx sol.mtx [tol=1.0e-7]" << std::endl;
	else
	{
		double tol = argc > 4? atof(argv[4]) : 1.0e-7;
		CSRMatrix A;
		std::vector<double> x,b;
		A.Load(std::string(argv[1]));
		LoadVector(std::string(argv[2]),b);
		if( A.Size() == b.size() )
		{
			LoadVector(std::string(argv[3]),x);
			if( A.Size() == x.size() )
			{
				double normb = 0;
				for(int k = 0; k < (int)b.size(); ++k)
					normb += b[k]*b[k];
				A.Multiply(-1.0,x,1.0,b);
				double norm = 0;
				for(int k = 0; k < (int)b.size(); ++k)
					norm += b[k]*b[k];
				norm = sqrt(norm/normb);
				std::cout << "norm " << norm << std::endl;
				if( norm <= tol ) return 0;
			}
		}
	}
	return 1;
}