function x = cholsolve(R,b)
% CHOLSOLVE  Solve a system after cholesky decomposition.
%
%   x = cholsolve(R,b) solves the system Ax = b, given
%   upper triangular matrix R obtained from the Cholesky
%   decomposition of positive definite matrix A.

[m, n] = size(R);
if m ~= n
	disp('The system is not square.');
	return;
end

% Solve the lower triangular system R'y = b
y = forsolve(R', b);
% Solve the upper triangular system Rx = y
% by back substitution
x = backsolve(R, y);
