%function cholesky_factorization ( )

  n = 3;

  l = [ 2,  0,  0; ...
       -1,  5,  0; ...
        4, -3,  6 ];

  a = l * l'
%
%  Compute the L factor.
%
  l = zeros ( n, n );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  L | L'' | L * L''\n' );
  fprintf ( 1, '\n' );
  [l,l',l*l']
  pause

  for i = 1 : n

    l(i,i) = sqrt ( a(i,i) - l(i,1:i-1) * l(i,1:i-1)' );
    for k = i + 1 : n
      l(k,i) = ( a(k,i) - l(k,1:i-1) * l(i,1:i-1)' ) / l(i,i);
    end

    fprintf ( 1, '  L | U | L * U after step %d\n', i );
    [ l, l', l * l' ]
    pause

  end


