function [c, v] = condest(A,t)
%CONDEST 1-norm condition number estimate.
%   C = CONDEST(A) computes a lower bound C for the 1-norm condition
%   number of a square matrix A.
%
%   C = CONDEST(A,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] = CONDEST(A) 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: CONDEST uses random numbers generated by RAND.  If repeatable
%   results are required,  use RNG to control MATLAB's random number
%   generator state.
%
%   CONDEST is based on the 1-norm condition estimator of Hager [1] and a
%   block oriented generalization of Hager's estimator given by Higham and
%   Tisseur [2].  The heart of the algorithm involves an iterative search
%   to estimate ||A^{-1}||_1 without computing A^{-1}. This is posed as the
%   convex, but nondifferentiable, optimization problem:
%
%         max ||A^{-1}x||_1 subject to ||x||_1 = 1.
%
%   See also NORMEST1, COND, NORM, RAND.

%   Reference:
%   [1] William W. Hager, Condition estimates,
%       SIAM J. Sci. Stat. Comput. 5, 1984, 311-316, 1984.
%
%   [2] Nicholas J. Higham and Fran\c{c}oise Tisseur,
%       A Block Algorithm for Matrix 1-Norm Estimation
%       with an Application to 1-Norm Pseudospectra,
%       SIAM J. Matrix Anal. App. 21, 1185-1201, 2000.
%
%   Nicholas J. Higham
%   Copyright 1984-2023 The MathWorks, Inc.

if size(A,1) ~= size(A,2)
    error(message('MATLAB:condest:NonSquareMatrix'))
end
if isempty(A)
    c = full(real(zeros('like',A)));
    v = zeros(0,1,'like',c);
    return
end
if nargin < 2
    t = [];
end
if issparse(A)
    [L,U,~,~] = lu(A,'vector');
else
    [L,U,~] = lu(A,'vector');
end
k = find(abs(diag(U))==0);

if isobject(A)
    % Can't call "nowarn" version of mldivide
    warns = warning('query','all');
    temp = onCleanup(@()warning(warns));
    warning('off','all');
end
if ~isempty(k)
    c = Inf;
    n = size(A,1);
    v = zeros(n,1);
    k = min(k);
    v(k) = 1;
    if k > 1
        v(1:k-1) = -localMldivide(U(1:k-1,1:k-1), U(1:k-1,k));
    end
else
    [Ainv_norm, ~, v] = normest1(@(flag,x)condestf(flag,x,L,U),t);
    A_norm = norm(A,1);
    c = Ainv_norm*A_norm;
end
v = v/norm(v,1);
end

function f = condestf(flag, X, L, U)
%CONDESTF   Function used by CONDEST.
if isequal(flag,'dim')
    f = size(L, 1);
elseif isequal(flag,'real')
    f = isreal(L) && isreal(U);
elseif isequal(flag,'notransp')
    f = localMldivide(U, localMldivide(L, X));
elseif isequal(flag,'transp')
    transp = true;
    f = localMldivide(L, localMldivide(U, X, transp), transp);
end
end

function x = localMldivide(A, b, transp)
if nargin < 3
    transp = false;
end
if isobject(A) || isobject(b)
    if transp
        x = A'\b;
    else
        x = A\b;
    end
else
    x = matlab.internal.math.nowarn.mldivide(A, b, transp);
end
end
