编一个程序，用C++实现矩阵的cholesky分解

给，下面是Cholesky分解法的C++经典算法：
//-------------------------------------------------------------------
// Cholesky分解法
//-------------------------------------------------------------------
template <class T>
int cholesky(Matrix<T>& mat, double epsilon=EPSILON) {
size_t i, j, k;
for (i=0; i<mat.Rows(); ++i) {
// 计算第 i 轮主元
for (k=0; k<i; ++k) mat[i][i] -= mat[k][i]*mat[k][i];
mat[i][i] = sqrt(mat[i][i]);
// 计算第 i 轮主元结束
if (fabs(mat[i][i])<epsilon) break;
//计算第 i 列
// for (j=i+1; j<mat.Rows(); ++j) {
// for (k=0; k<i; ++k) mat[j][i] -= mat[j][k]*mat[i][k];
// }
// 计算第 i 列结束
// 计算第 i 行
for (j=i+1; j<mat.Cols(); ++j) {