Cholesky algorithm

Question

Good afternoon everyone, I'm in need of a factoring algorithm cholesky and algorithms to solve upper and lower triangular systems, but I'm not finding any work in that octave. Recalling that need the algorithm and can not just use the already available functions in the programs. I even tried to adapt this algorithm but it is in trouble in the windows octave

function [G] = cholesky(A,dim)
epsilon = 1.0e-10;
sum = 0.0;
                for k = 1:(i-1)
                sum = sum + G(k,i)*G(k,i);
                end
        delta = A(i,i) - sum;
        if( delta < 0.0 ), errorspd(delta,i); end
          G(i,i) = sqrt(delta);
         if( abs(G(i,i)) < epsilon );  end
         for j = (i+1):ndim
                        sum = 0.0;
                        for k = 1:(i-1)
                        sum = sum + G(k,i)*G(k,j);
                        end
            G(i,j) = ( A(i,j) - sum ) / G(i,i);
             end
        end

Thanks in advance

Answer 1

Here is my Matlab code for Cholesky, I hope it works also on Octave. It is taken step by step by the wikipedia Cholesky–Banachiewicz algorithm:

     function[L]=MyChol(A)

     [n,m]=size(A);
     L=eye(n);
     for k=1:n-1
           L(k,k)=sqrt(A(k,k));  %Computing the diagonal
           L(k+1:n,k)=(A(k+1:n,k))/L(k,k); %Computing the lower part
           A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-L(k+1:n,k)*L(k+1:n,k)'; %Putting the new A for the next step on the loop
     end
     L(n,n)=sqrt(A(n,n));
     end

The $L$ obtained is such that $A=LL^T$. For solving upper triangular systems of the form $Ux=b$, I use the next code based on the backward substitution

    function [x] = triup(U,b)

    n = length(b);
    x = zeros(n,1);

    % backward substitution
    x(n) = b(n)/U(n,n);
    for k = n-1:-1:1
        x(k) = ( b(k) -U(k,(k+1):n)*x((k+1):n))/U(k,k);
    end
    % end of backward substitution

    end

Meanwhile for lower ones of the form $Lx=b$, I use an algorithm base on the forward substitution

   function [x] = trilo( L,b )

   n=length(L);
   x=zeros(n,1);

   % forward substitution
   for k=1:n-1
       x(k)=b(k)/L(k,k); 
       b(k+1:n)=b(k+1:n)-x(k)*L(k+1:n,k);
   end
   x(n)=b(n)/L(n,n);
   % end of forward substitution

   end

As an example, if you want to solve a system $Ax=b$ where $A$ is symmetric positive definite matrix, then use

   L=MyChol(A);
   y=trilo(L,b);
   x=triup(L',y);

to obtain the solution $x$.