#include <iostream>
#include <cstdlib>
#include <cmath>
#include <fstream>

using namespace std;

int main()
{
	int n;
	ifstream in("matrix.dat");
	in>>n;
	
	double *b  = new double[n];
	double *x  = new double[n];
	double **a = new double*[n];
	for(int i=0;i<n;i++)
		a[i] = new double[n];

	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
			in>>a[i][j];
		in>>b[i];
	}
	in.close();
	int niter;
	double omega;
	cout<<"Enter Number of Maximum Iterations"<<endl;
	cin>>niter;
	cout<<"Enter Relaxation Factor"<<endl;
	cin>>omega;
		
	for(int iter=0;iter<niter;iter++)
	{
		for(int i=0;i<n;i++)
		{
			double sum = 0.e0;
			for(int j=0;j<n;j++)
			{
				if(i!=j)
					sum = sum+a[i][j]*x[j];
			}
			x[i]=(b[i]-sum)/a[i][i]*omega+x[i]*(1.e0-omega);
		}
		double residual=0.e0;
		for(int i=0;i<n;i++)
		{
			double sum = 0.e0;
			for(int j=0;j<n;j++)
				sum = sum + a[i][j]*x[j];
			residual = residual+(sum-b[i])*(sum-b[i]);
		}
	cout<<" Iteration "<<iter<<" Residual "<<sqrt(residual)<<endl;
	}
	
	ofstream out("solution.dat");
	for(int i=0;i<n;i++)
	{
		out<<i<<"\t"<<x[i]<<endl;
	}
	out.close();
//	for(int i=0;i<n;i++)
//	{
//		for(int j=0;j<n;j++)
//			cout<<a[i][j]<<"\t";
//		cout<<b[i]<<endl;
//	}
}