
%function [a,x,rd,infopd]=chfac(a,d)
 function [infopd]=chfac(a)
% line changed and d modified below
%
%
%     given an n by n symmetric matrix a, an n by n diagonal
%     matrix d, this subroutine attempts to
%     compute an upper triangular matrix r such that
%
%                          t
%                 a + d = r r.
%
%
%     if the cholesky decomposition of a + d does not exist,
%     then there is a largest leading submatrix of a + d which
%     is not positive definite, and this subroutine determines
%     the order of this submatrix and a direction of negative
%     curvature for a + e. more precisely, infopd is the order
%     of this submatrix, and x is such that
%
%                 (x,(a + e)*x) .le. zero
%
%     where (x,y) denotes the inner product between vectors
%     x and y.
%
%
%     where
%
%
%     a is an n by n array. on input the full lower triangle
%       must contain the full lower triangle of the matrix a.
%       if the cholesky decomposition exists, then on output
%       the strict upper triangle contains the strict upper
%       triangle of the matrix r. otherwise, on output the
%       first infopd columns in the strict upper triangle
%       contains the strict upper triangle of the cholesky
%       factor for the order infopd leading submatrix of a + e.
%
%     d is an input array of length n which must contain
%       the diagonal elements of the matrix d.
%
%
%     x is an output array of length n.  if the
%       choleskey decomposition does not exist then
%       x is a vector
%       such that (x,(a + e)*x) is not positive.
%
%     rd is an output array of length n. if the cholesky
%       decomposition exists, then rdiag contains the positive
%       diagonal elements of the matrix r. otherwise the first
%       infopd elements of rdiag contain the positive diagonal
%       elements of the cholesky factor for the order infopd
%       leading submatrix of a + e, and rdiag(infopd+1) contains
%       (x,(a + e)*x).
%
%     infopd is an integer output variable which specifies the
%       order of the largest leading submatrix of a + d
%       which is positive definite. if no leading submatrix
%       is positive definite, then infopd = 0.
%
%     **********

%n=length(d); 
 n=size(a,1);
 d=zeros(n,1);
 a=tril(a); da=diag(diag(a)); a=a-da; a=a+a'+da; 
 infopd=0; phs=1; j=2; rd=zeros(n,1); x=rd; wa=rd;
 temp=a(1,1)+d(1);
 if temp<=0, phs=0;
 else,  rd(1)=sqrt(temp); infopd=1; end
 while j<=n & phs==1
     sum2=0;
     for i=1:j-1
         if i==1, sum1=0; else, sum1=a(1:i-1,j)'*a(1:i-1,i); end
         a(i,j)=(a(i,j)-sum1)/rd(i);
         sum2=sum2+a(i,j)^2;
     end
     temp=a(j,j)+d(j)-sum2;
     if temp<=0, phs=0;
     else, rd(j)=sqrt(temp); infopd=j; end
     j=j+1;
 end      
%  change bottom since we only need to know if pos def.

%'inside chfac'
%keyboard
 if infopd<n
     infopd=-1;
%    rd(infopd+1)=temp;
%    x(infopd+1)=1;
%    if infopd~=0
%       b=a(1:infopd,infopd+1);
%       B=triu(a(1:infopd,1:infopd))-diag(diag(a(1:infopd,1:infopd)))+diag(rd(1:infopd));
%       x(1:infopd)=-B\b;
%    end
 else
      infopd=1;
 end
