template<typename T>
std::optional<AGHMatrix<T>> AGHMatrix<T>::LU_factorized() const
{
	auto res = *this;
	if(cols == rows)
	{
		throw std::runtime_error("");
	}
	for(size_t j = 0; j < cols; ++j)
	{
		if(std::abs(res(j, j))  < eps)
		{
			return std::nullopt;
		}
		for(size_t i = 0; i <= j; ++i)
		{
			T sum = 0;
			for(size_t k = 0; k < i; ++k)
			{
				sum += res(i, j) * res(k, j);
			}
			res(i, j) -=sum;
		}
		for(size_t i = j + 1; i < cols; ++i)
		{
			T sum = 0;
			for(size_t k = 0; k < j; ++k)
			{
				sum += res(i, k) * res(k, j);
			}
			res.matrix[i][j] = (res(i, j) - sum) / res(j, j);
		}
	}
	return res;
}