function main_cholesky

%DRIVER_cholesky_easy_to_read
%DRIVER_cholesky_efficient
DRIVER_LDL

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 - trans(R)*A*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 - trans(R)*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];
A = randn(7);
A = A*A';
m = size(A,2);

%%% Perform Cholesky factorization.
R = A;
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

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

function DRIVER_LDL

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

m = size(A,2);
L = eye(m);
D = A;

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

  %%% Print results
  fprintf(1,'D = \n')
  disp(D)
  fprintf(1,'L = \n')
  disp(L)
  fprintf(1,'|A - L*D*trans(L)| = %10.3e\n',max(max(abs(A - L*D*L'))))
  keyboard
  
end

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

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


