> restart; 
> with(LinearAlgebra):
> cholesky:=proc(A) 
> local G,n,j,k,r,g; 
> # option trace;
> n:=RowDimension(A); 
> G:=Matrix(n,n);
> for j from 1 to n  do
>     g:=A[j,j];
>     for r from 1 to j-1 do g:=g - G[j,r]^2 end do;
>     G[j,j]:=sqrt(g);
>     for k from j+1 to n do 
>        g:=A[k,j];
>        for r from 1 to j-1 do g:=g-G[j,r]*G[k,r] end do;
>        G[k,j]:=g/G[j,j]
>        end do;
> end do; 
> evalm(G)
> end proc:
> B:=Matrix([[1,2,3],[0,4,5],[0,0,6]]);
> A:=MatrixMatrixMultiply(Transpose(B),B);
> C:=cholesky(A);

>