function d=svdFrancis(A)

TOL=100*eps;
d=diag(A);
n=length(d);
for i=1:n-1
    e(i)=A(i,i+1);
end
maxit=500*n^2;

lambda(n)=abs(d(n));
for j=n-1:-1:1
  lambda(j)=abs(d(j))*lambda(j+1)/(lambda(j+1)+abs(e(j)));
end
mu(1)=abs(d(1));
for j=1:n-1
  mu(j+1)=abs(d(j+1))*mu(j)/(mu(j)+abs(e(j)));
end
sigmaLower=min(min(lambda),min(mu));
thresh=max(TOL*sigmaLower,maxit*realmin);


iUpper=n-1;
iLower=1;
for iterations=1:maxit

  for i=iUpper:-1:1
    iUpper=i;
    if abs(e(i))>thresh
      break;
    end
  end
  j=iUpper;
  for i=iLower:iUpper
    if abs(e(i))>thresh
      j=i;
      break;
    end
  end
  iLower=j;

  if (iUpper==iLower & abs(e(iUpper))<=thresh) | ...
        (iUpper<iLower)
    d=sort(abs(d));
    d(end:-1:1)=d;
    return
  end

  [d(iLower:iUpper+1),e(iLower:iUpper)]= ...
        IQR(d(iLower:iUpper+1),e(iLower:iUpper));
end
error('mysvd: too many iterations!')
