#include <mexhelper.h>

static bool verbose = false;

template<class T> myVector<T> gaussElimination (matrix<T> A, myVector<T> b)
{
	if (verbose)
		std::cerr << "A:\n" << A << "\nb:\n" << b << "\n\n";
	int N = A.N();

	for (int j=1;j<N;j++)
		for (int i=j;i<N;i++)
		{
			T m = A.get(i, j-1) / A.get(j-1, j-1);
			A.set(i, A.get(i) - A.get(j-1) * m);
			b[i] = b[i] - m*b[j-1];
		}

	myVector<T> x(N);
	for (int i=0;i<N;i++)
		x[i] = 0;

	x[N] = b[N] / A.get(N, N);
	for (int j=N-1;j>=0;j--)
		x[j] = (b[j] - A.get(j, j+1, N) * x.get(j+1, N)) / A.get(j, j);

	return x;
}