mudd.m

GPS software toolbox for GPS receiver development

%                           mudd.m
%  Scope:   This MATLAB macro computes the U-D factorization of a real
%           symmetric, positive (semi)definite matrix by using the
%           modified Cholesky decomposition.
%  Usage:   [ud,ier] = mudd(n,a,meps)
%  Description of parameters:
%           n    - input, real scalar, number of rows and columns in  A
%           a    - input, real array of length  n*(n+1)/2 , storing the
%                  upper triangular part of the real symmetric matrix  A,
%                  columnwise
%           meps - input, real scalar, machine dependent accuracy
%           ud   - output, real array of length  n*(n+1)/2, storing the
%                  upper triangular part of the pair  U-D ; the elements of
%                  the diagonal matrix  D  are stored on diagonal locations
%                  of the unit upper triangular matrix  U
%           ier  - output, integer scalar, index of error
%                  = 0    no computation error
%                  = 1    computation error occurred, a diagonal element
%                         is smaller than  meps 
%                  = 2    computation error occurred, the input matrix  a
%                         is not numerically positive definite            
%  Last update:  06/29/00
%  Copyright (C) 1996-00 by LL Consulting. All Rights Reserved.

function  [ud,ier] = mudd(n,a,meps)

ier = 0;
p = n * (n + 1) / 2;
ud = zeros(p,1);
ud(p) = a(p);
if n == 1 
   return
end 
if ud(p) <= meps
   ier = 1;
   return
end 
nm1 = n - 1;
d = 1. / ud(p);
for k = 1:nm1
   p = p - 1;
   ud(p) = a(p) * d;
end
for i = nm1:-1:1
   pp = p - 1;
   for j = i:-1:1
      p = p - 1;
      pi = pp;
      pk = pp;
      s = a(p);
      if p == pp
         for k = i:nm1
            pi = pi + k;
            pk = pk + k + 1;
            s = s - ud(pi) * ud(pi) * ud(pk);
         end
         ud(p) = s;
         if abs(s) <= meps
            ier = 2;
            return
         end 
         d = 1. / s;
      else             
         pj = p;
         for  k = i:nm1
            pj = pj + k;
            pi = pi + k;
            pk = pk + k + 1;
            s = s - ud(pj) * ud(pi) * ud(pk);
         end
         ud(p) = s * d;
      end
   end
end