function [c, v] = condestdec(L,U,t)
%CONDESTDEC 1-norm condition number estimate using matrix decomposition
%   C = CONDESTDEC(L,U) computes a lower bound C for the 1-norm condition
%   number of a square matrix A = L*U where 
%   L is lower triangular, U is upper triangular (or row permutations of triangular
%   matrices).
%
%   C = CONDESTDEC(L,U,T) changes T, a positive integer parameter equal to
%   the number of columns in an underlying iteration matrix.  Increasing the
%   number of columns usually gives a better condition estimate but increases
%   the cost.  The default is T = 2, which almost always gives an estimate
%   correct to within a factor 2.
%
%   [C,V] = CONDESTDEC(L,U) also computes a vector V which is an approximate null
%   vector if C is large.  V satisfies NORM(A*V,1) = NORM(A,1)*NORM(V,1)/C.
%
%   Note: CONDESTDEC invokes RAND.  If repeatable results are required then
%   invoke RAND('STATE',J), for some J, before calling this function.
%
%   Uses block 1-norm power method of Higham and Tisseur.
%
%   See also NORMEST1, COND, NORM.

%  modified from condest by tvp.
%   Nicholas J. Higham, 9-8-99
%   Copyright 1984-2001 The MathWorks, Inc. 
%   $Revision: 5.16 $  $Date: 2001/04/15 12:01:38 $


if size(L,1) ~= size(L,2)
   error('Matrix must be square.')
end
if size(U,1) ~= size(U,2)
   error('Matrix must be square.')
end
if isempty(L), c = 0; v = []; return, end
if nargin < 3, t = []; end
warns = warning;
warning('off')


k = find(abs(diag(U))==0);
if ~isempty(k)
   c = Inf;
   n = max(size(L));
   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
else
   [Ainv_norm, temp, v] = normest1(@condestf,t,[],L,U);
   A_norm = norm(L*U,1);
   c = Ainv_norm*A_norm;
end
v = v/norm(v,1);
warning(warns)

function f = condestf(flag, X, L, U)
%CONDESTF   Function used by CONDEST.

if isequal(flag,'dim')
   f = max(size(L));
elseif isequal(flag,'real')
   f = isreal(L) & isreal(U);
elseif isequal(flag,'notransp')
   f = U\(L\X);
elseif isequal(flag,'transp')
   f = L'\(U'\X);
end
