% cholesky.m: Cholesky factorization for solving the linear system 
%             of equations  A*x = b. 

  clear;
  format short;

% Step 0: Assign the matrix A and the vector b. 

  n = 4;
  A = [ 0.05, 0.07, 0.06, 0.05;  
        0.07, 0.10, 0.08, 0.07;
        0.06, 0.08, 0.10, 0.09;
        0.05, 0.07, 0.09, 0.10 ];
  b = [ 0.23, 0.32, 0.33, 0.31 ];

% Step 1: Cholesky factorization: A = L*L' where L is a lower 
%         triangular matrix and L' is the tranpose matrix of L.

  L = zeros(n,n);       % initialize the n x n matrix L
  U = zeros(n,n);       % initialize the n x n matrix U 
  for k = 1:n
    L(k,k) = A(k,k);
    for s = 1:(k-1)
      L(k,k) = L(k,k)-L(k,s)*L(k,s);
    end
    L(k,k) = sqrt(L(k,k));
    for i = (k+1):n
%     (fill out a few lines here to compute L(i,k))
    end
  end
  L                     % output L

% Step 2: Forward subsitiution to solve  L*y = b. 

  y = zeros(n,1);       % initialize the vector y
  y(1) = b(1)/L(1,1);
  for i = 2:n
    y(i) = b(i);
    for j = 1:(i-1)
      y(i) = y(i)-L(i,j)*y(j);
    end
    y(i) = y(i)/L(i,i);
  end
  y                     % output y

% Step 3: Back subsitiution to solve  U*x = y. 

  x = zeros(n,1);       % initialize the vector x
  x(n) = y(n)/L(n,n);
  for i = (n-1):-1:1
%   (fill out a few lines here to compute x(i))
  end 
  x                     % output x