mychol.m

LU decomposition routines in matlab

function [A,iflag] = mychol( A )
%
% Usage:
%        [A,iflag] = mychol( A )
%
% Given a symmetric N by N matrix A, CHOL attempts
% to compute the Cholesky decomposition of A.
% (Inner product form.)
% 
%
% input:
%        A:     The N by N matrix A.
%               Only the lower triangular part of A
%               is used.
%
% output:
%        A:     The upper triangular part of A
%               contains the Cholesky factor R.
%               The strict lower triangular part of A is unchanged.
%
%        iflag: error flag
%               iflag = 0  Cholesky decomposition could be computed
%                          and is stored in the upper triangular
%                          part of A.
%               iflag = 1  A is not squared
%               iflag > 1  Negative entry observed in step iflag-1.
%                          Matrix A is not positive definite.
%
%
% Matthias Heinkenschloss
% Jan 30, 2007

iflag = 0;

% get size of A and check dimensions
[m,n]    = size(A);
if ( m ~= n )
   iflag = 1;
   return
end


% Start Cholesky decomposition

% first step
if( A(1,1) <= 0 )
    iflag = 2;
    return
end

% Step 1
A(1,1)   = sqrt(A(1,1));
A(1,2:n) =  A(1,2:n) / A(1,1);

% Steps 2 to n
for i = 2:n

    t = A(i,i) - A(1:i-1,i)'*A(1:i-1,i);
    if( t <= 0 )
        iflag = i;
        return
    end
    A(i,i) = sqrt(t);
  
    for j = i+1:n
        A(i,j) = ( A(i,j) - A(1:i-1,i)'*A(1:i-1,j) ) / A(i,i);
    end
end
    
%end of mychol