function main_cholesky

%DRIVER_cholesky_easy_to_read
DRIVER_cholesky_efficient

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_cholesky_easy_to_read

A = [ 4, 1,-2, 1;...
      1, 3, 2, 1;...
     -2, 2, 4, 1;...
      1, 1, 1, 4];

m      = size(A,2);
A_orig = A;
R      = eye(m);

%%% Perform Cholesky factorization.
for k = 1:(m-1)
     
  J      = (k+1):m;
  alpha  = sqrt(A(k,k));
  ell    = (1/alpha)*A(k,J)';
  R(k,k) = alpha;
  R(k,J) = ell;  
  A(J,J) = A(J,J) - ell*ell';
  A(k,k) = 1;
  A(J,k) = zeros(length(J),1);
  A(k,J) = zeros(1,length(J));

  %%% Print results
  fprintf(1,'A = \n')
  disp(A)
  fprintf(1,'R = \n')
  disp(R)
  fprintf(1,'|A_orig - R*A*trans(R)| = %10.3e\n',max(max(abs(A_orig - R'*A*R))))
  keyboard
  
end
R(m,m) = sqrt(A(m,m));

fprintf(1,'\n\nUpon completion, we have: R = \n')
disp(R)
fprintf(1,'|A_orig - R*trans(R)| = %10.3e\n',max(max(abs(A_orig - R'*R))))
keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_cholesky_efficient

A = [ 4, 1,-2, 1;...
      1, 3, 2, 1;...
     -2, 2, 4, 1;...
      1, 1, 1, 4];

m = size(A,2);
R = A;

%%% Perform Cholesky factorization.
for k = 1:(m-1)
     
  J      = (k+1):m;
  alpha  = sqrt(R(k,k));
  ell    = (1/alpha)*R(k,J)';
  R(k,k) = alpha;
  R(k,J) = ell';  
  R(J,J) = R(J,J) - ell*ell'; % We really only need the lower half here!
  
end
R(m,m) = sqrt(R(m,m));
R      = triu(R);

fprintf(1,'Upon completion, we have: R = \n')
disp(R)
fprintf(1,'|A - trans(R)*R| = %10.3e\n',max(max(abs(A - R'*R))))
keyboard
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

