#include "chapter6.cpp"

typedef double runtype;
typedef std::vector<runtype> Array;
typedef std::vector<Array> Matrix;

typedef std::vector<int> IntArray;



void main()
{
ofstream outfile1("ngauss.txt");
	int m = 10;
	Matrix A(m+1, Array(m+1, 0.0));
	Array b(m+1, 0.0);
	Array x(m+1, 0.0);
	int i, j, n;
	int temp_i, temp_j;

	for (n = 4; n <= 10; n++)
	{
		for (i = 1; i <= n; i++)
		{
			for (j = 1; j <= n; j++)
			{
				temp_i = i + 1;
				temp_j = j - 1;
				A[i][j] = pow(temp_i, temp_j);
			}
			b[i] = (pow(temp_i, n) - 1) / i;
		}
		naive_gauss(n, A, b, x);
		for (i=1; i<=n; i++)
			outfile1 << "n =  "<<setw(2)<< n << " x["<<setw(2)<< i << "] = " << setw(15) << x[i] << endl;
	}
	outfile1 << endl;
	A[1][1] = 3.0;  A[1][2] = -13.0;  A[1][3] = 9.0;  A[1][4] = 3.0;
	A[2][1] = -6.0; A[2][2] = 4.0;    A[2][3] = 1.0;  A[2][4] = -18.0;
	A[3][1] = 6.0;  A[3][2] = -2.0;   A[3][3] = 2.0;  A[3][4] = 4.0;
	A[4][1] = 12.0; A[4][2] = -8.0;   A[4][3] = 6.0;  A[4][4] = 10.0;
	b[1] = -19.0;
	b[2] = -34.0;
	b[3] = 16.0;
	b[4] = 26.0;
	m = 4;

	naive_gauss(4, A, b, x);
	for (i=1; i<=m; i++)
		outfile1 << " x["<<setw(2)<< i << "] = " << setw(15) << x[i] << endl;




ofstream outfile2("gauss.txt");

	IntArray l(m+1, 0);
   
	outfile2 << "Matrix A (4 x 4) = " <<endl;
	for (i = 1; i <= m; i++)
	{
		for (j = 1; j <= m; j++)
			outfile2 <<setw(5)<< A[i][j];
		outfile2 << endl;
	}
	gauss(m, A, l);
	outfile2 << "Matrix A after forward elimination = " << endl;
	for (i = 1; i <= m; i++)
	{
		for (j = 1; j <= m; j++)
			outfile2 <<setw(10)<< A[i][j];
		outfile2 << endl;
	}
	solve(m, A, l, b, x);
	outfile2 << "Vector X = " << endl;
	for (i = 1; i <= m; i++)
		outfile2 << x[i] << endl;


ofstream outfile3("tridiag.txt");

	m = 10; /* assuming square matrix */
	Array a(m+1, 0.0);
	Array c(m+1, 0.0);
	Array d(m+1, 0.0);

 /* setup matrix to test result */
	for (i = 1; i <= m; i++)
	{
		a[i] = 0.5;
		b[i] = 3.0;
		c[i] = 0.5;
		d[i] = 2.0;
	}
	b[1] = 2.5;
	b[m] = 2.5;
	tri (m, a, d, c, b, x);
	outfile3 << "Result from module tri (1)" << endl;
	for (i = 1; i <= m; i++)
	{
		outfile3 << x[i] << endl;
		d[i] = 2.0;
		c[i] = 0.5;
		b[i] = 3.0;
	}

	outfile3 << endl << "Result from module tri (2)" << endl;
	b[1] = 2.5;
	b[m] = 2.5;
	tri (m, c, d, c, b, x);
	for (i = 1; i <= m; i++)
		outfile3 << x[i] << endl;



ofstream outfile4("penta.txt");
   
	m = 10; /* assume square matrix */
	Array e(m+1, 0.0);
	Array f(m+1, 0.0);


  /* setup matrix to test result */
	for (i = 1; i <= m; i++)
	{
		a[i] = 0.25;
		b[i] = 2.0;
		c[i] = 0.25;
		d[i] = 1.0;
		e[i] = 0.25;
		f[i] = 0.25;
	}

	b[1] = 1.5;  
	b[2] = 1.75;      
	b[m - 1] = 1.75;
	b[m] = 1.5;
	penta(m, e, a, d, c, f, b, x);
	outfile4 << "Result from module Penta" << endl;
	for (i = 1; i <= m; i++)
		outfile4 << x[i] << endl;


ofstream outfile5("jacobi.txt");

	n = 3;
	Matrix a3(n+1, Array(n+1, 0.0));
	Array b3(n+1, 0.0);
	Array x3(n+1, 0.0);

	a3[1][1] = 2; a3[1][2] = -1; a3[1][3] = 0;
	a3[2][1] = -1;a3[2][2] = 3;  a3[2][3] = -1;
	a3[3][1] = 0; a3[3][2] = -1; a3[3][3] = 2;

	b3[1] = 1; b3[2] = 8; b3[3] = -5;
  
	x3[1] = 1; x3[2] = 2; x3[3] = 3;

	outfile5 << "Let A be:" << endl;
	for (i = 1; i <= n; i++)
	{
		for(j = 1; j <= n; j++)
			outfile5 << setw(10) << a3[i][j];
		outfile5 << endl;
	}
	outfile5 << "Let b be:" << endl;
	for (i = 1; i <= n; i++)
		outfile5 << b3[i] << endl;

	outfile5 << "Jacobi iteration:" << endl;
	jacobi(outfile5, n, a3, b3, x3);
}