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

  clear;
  format short;

% Step 1: 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 2: 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
      L(i,k) = A(i,k);
      for s = 1:(k-1)
        L(i,k) = L(i,k)-  % Complete this line. Don't forget the semicolon.
      end
      L(i,k) = L(i,k)/L(k,k);     
    end
  end
  L                     % output L

% Step 3: Forward substitute 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 4: Back substitute to solve L'*x = y. 

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

  x                     % output x