function L=secfac(n,x,x0,g,g0,L,epsm,iter,rnf,iagflg)
s=x-x0;
y=g-g0;

if iter==1; l=sqrt(s'*y/(s'*s))*L; 
else; 
l=L;
end;
L=chol( l*( eye(n)-(l'*s*s'*l)/(s'*l*l'*s)+(l\y)*(l\y)'/(y'*s))*l')';

return
den1=s'*y;
snrm=norm(s);
ynrm=norm(y);
%        [den1,sqrt(epsm)*snrm*ynrm]
if den1>=sqrt(epsm)*snrm*ynrm;
  u=L'*s;
  den2=u'*u;
  alpha=sqrt(den1/den2);

  if iter==1;
    u=alpha*u;
    L=alpha*L;
    den2=den1;
    alpha=1;
  end;

  w=L*u;
  if iagflg==0;
    reltol=sqrt(rnf);
  else;
    reltol=rnf;
  end;
  if any(abs(y-w)<reltol*max(g0,g));
    return;
  end;
  w=y-alpha*w;
  alpha=alpha/den1;
  u=alpha*u;
  k=1;
  tem=find(u~=0);
  if ~isempty(tem); k=tem(length(tem)); end;
  if k>1;
    for i=1:k-1;
      j=k-i;
      if u(j)==0;
        L(j:n,[j,j+1])=L(j:n,[j+1,j]);
        u(j)=u(j+1);
      else;
        tem=sqrt(u(j)*u(j)+u(j+1)*u(j+1));
        t1=u(j)/tem;
        t2=-u(j+1)/tem;
        L(j:n,[j,j+1])=L(j:n,j)*[t1,t2]+L(j:n,j+1)*[-t2,t1];
        u(j)=tem;
      end;
    end;
  end;
  L(:,1)=L(:,1)+u(1)*w(:);
  if k>1;
    for i=1:k-1;
      if L(i,i)==0;
        L(i:n,[i,i+1])=L(i:n,[i+1,i]);
      else;
        tem=sqrt(L(i,i)*L(i,i)+L(i,i+1)*L(i,i+1));
        t1=L(i,i)/tem;
        t2=-L(i,i+1)/tem;
        L(i:n,[i,i+1])=L(i:n,i)*[t1,t2]+L(i:n,i+1)*[-t2,t1];
      end;
    end;
  end;
end;   

L
ltem

L-ltem
