function [c, v] = condestdec(L,U)
%CONDEST 1-norm condition number estimate.
%   [C,V] = CONDESTDEC(L,U) computes a lower bound C for the 1-norm 
%   condition number of A with LU decomposition L*U, and a vector V such that 
%       NORM(A*V,1)=NORM(A,1)*NORM(V,1)/C.
%   V is an approximate null vector if C is large.
%
%   Uses Higham's modification of Hager's method.
%
%   See also COND, NORM, NORMEST.

%   Reference: N.J. Higham, FORTRAN codes for estimating the one-norm of
%   a real or complex matrix, with applications to condition estimation,
%   ACM Trans. Math. Soft., 14 (1988), pp. 381-396.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.9 $  $Date: 1997/11/21 23:38:21 $


warns = warning;
warning('off')
[n,n] = size(L);
k = find(abs(diag(U))==0);
if ~isempty(k)
   c = Inf;
   v = zeros(n,1);
   k = min(k);
   v(k) = 1;
   if k > 1
      v(1:k-1) = U(1:k-1,1:k-1)\(-U(1:k-1,k));
   end
   v = v/norm(v,1);
else
   x = ones(n,1)/n;
   v = U\(L\x);      % Av=x
   gamma = norm(v, 1);
   xi = sign(v); f = find(xi==0); xi(f) = 1;   % set sign(0) := 1
   lastsign = xi;
   x = L'\(U'\xi);   % A'x=xi
   x = real(x);      % Needed only for the complex case.
   
   maxits = 5*(~isempty(U));
   for k = 2:maxits
   
       % Get smallest index for which max(abs(z)) is attained.
       [big, j] = max(abs(x)); f = find(abs(x)==big); j = f(1);
       x = zeros(n,1); x(j) = 1;
   
       v = U\(L\x);      % Av=x
       xi = sign(v); f = find(xi==0); xi(f) = 1;   % set sign(0) := 1
   
       vold = gamma;
       gamma = norm(v, 1);
   
       if xi == lastsign | gamma <= vold
          break
       end
       lastsign = xi;
   
       x = L'\(U'\xi);   % A'x=xi
       x = real(x);      % Needed only for the complex case.
   
       if norm(x, inf) == x(j), break, end
       if k == maxits
          disp(['Warning: Not converged after ' int2str(maxits) ' iterations.'])
       end
   end
   
   x = (0:n-1)';
   x = (1-2*rem(x,2)) .* (1 + x/max(1,n-1));
   x = U\(L\x);
   temp = 2*norm(x,1)/(3*n);
   if temp > gamma
      v = x; gamma = temp;
   end
   
   anorm = norm(L*U,1);
   c = gamma*anorm;
   v = v/norm(v,1);
end
warning(warns)
