procedure Cholesky (n : int, S : array 1..*,1..* of real,
                             var L : array 1..*,1..* of real)

        %  This procedure computes the lower triangular Cholesky factor L
        %  of a real n x n matrix S.  S must be symmetric positive-definite.
        %  If it is not, the fatal errors of divide-by-0 or square root of
        %  a negative number could occur.
        %
        %  Upon return, S = L x transpose(L).
        %
        %  Elements above the main diagonal of S are not accessed.

        var sum : real

        % L is constructed column-wise
        for j : 1 .. n
            % elements above main diagonal are zero
            for i : 1 .. j-1
                L(i,j) := 0.0
            end for
            % compute element on main diagonal
            sum := 0
            for k : 1 .. j-1
                sum := sum + L(j,k)**2
            end for
            L(j,j) := sqrt(S(j,j) - sum)
            % compute elements below main diagonal
            for i : j+1 .. n
                sum := 0
                for k : 1 .. j-1
                    sum := sum + L(i,k)*L(j,k)
                end for
                L(i,j) := (S(i,j) - sum)/L(j,j)
            end for
        end for

    end Cholesky