49 template <
class SPDMatrix,
class SymmMatrix>
59 if (M != L.dim(1) || N != L.dim(2))
65 typename SPDMatrix::element_type dot=0;
73 dot = dot + L(j,i)*L(j,i);
75 L(j,j) = A(j,j) - dot;
77 for (i=j+1; i<=N; i++)
81 dot = dot + L(i,k)*L(j,k);
82 L(i,j) = A(j,i) - dot;
85 if (L(j,j) <= 0.0)
return 1;
87 L(j,j) = sqrt( L(j,j) );
89 for (i=j+1; i<=N; i++)
90 L(i,j) = L(i,j) / L(j,j);
TNT_SUBSCRIPT_TYPE Subscript
int Cholesky_upper_factorization(SPDMatrix &A, SymmMatrix &L)