Simple matrix class with LU decomposition

Simple matrix class with LU decomposition.cpp

#pragma once

#include <tuple>

using std::tuple;
using std::make_tuple;

template<typename T, int m, int n>
class DenseMatrix
{
	T matrix[m * n];
public:
	T & operator() (int i, int j) { return matrix[i*n + j]; }
	const T& operator() (int i, int j) const { return matrix[i*n + j]; }

	template<typename TOtherMatrix>
	DenseMatrix<T, m, n> operator + (const TOtherMatrix& b);

	template<typename TOtherMatrix>
	DenseMatrix<T, m, n> operator - (const TOtherMatrix& b);
};

template<typename T, int m, int n>
template<typename TOtherMatrix>
DenseMatrix<T, m, n> DenseMatrix<T, m, n>::operator + (const TOtherMatrix& b)
{
	DenseMatrix<T, m, n> res;
	for (int i = 0; i < m; ++i)
	{
		for (int j = 0; j < n; ++j)
		{
			res(i, j) = (*this)(i, j) + b(i, j);
		}
	}
	return res;
}

template<typename T, int m, int n>
template<typename TOtherMatrix>
DenseMatrix<T, m, n> DenseMatrix<T, m, n>::operator - (const TOtherMatrix& b)
{
	DenseMatrix<T, m, n> res;
	for (int i = 0; i < m; ++i)
	{
		for (int j = 0; j < n; ++j)
		{
			res(i, j) = (*this)(i, j) - b(i, j);
		}
	}
	return res;
}

template<typename T, int m, int l, int n>
DenseMatrix<T, m, n> operator * (const DenseMatrix<T, m, l>& A, const DenseMatrix<T, l, n>& B)
{
	DenseMatrix<T, m, n> AB;
	for (int i = 0; i < m; ++i)
	{
		for (int j = 0; j < n; ++j)
		{
			T S = T();
			for (int k = 0; k < l; ++k)
			{
			    S += A(i,k)*B(k,j);
			}
			AB(i, j) = S;
		}
	}
	return AB;
}

template<typename T, int m, int n>
DenseMatrix<T, m, n> zero()
{
	DenseMatrix<T, m, n> Z;
	for (int i = 0; i < m; ++i)
	{
		for (int j = 0; j < n; ++j)
		{
			Z(i, j) = T();
		}
	}
	return Z;
}

template<typename T, int n>
tuple<DenseMatrix<T, n, n>, DenseMatrix<T, n, n> > lu_decompose(DenseMatrix<T, n, n>& A)
{
	DenseMatrix<T, n, n> L = zero<T,n,n>(), U = zero<T, n, n>();
	for (int i = 0; i < n; ++i)
	{
		U(0, i) = A(0, i);
		L(i, 0) = A(i, 0) / A(0, 0);
	}

	for (int i = 1; i < n; ++i)
	{
		for (int j = i; j < n; ++j)
		{
			T SU = T();
			T SL = T();
			for (int k = 0; k < i; ++k)
			{
				SU += L(i, k) * U(k, j);
				SL += L(j, k) * U(k, i);
			}
			U(i, j) = A(i, j) - SU;
			L(i, j) = (A(i, j) - SL) / U(i, i);
		}
	}
	return make_tuple(L, U);
}