function A = condex(n, k, theta, classname)
%CONDEX "Counter-examples" to matrix condition number estimators.
%   GALLERY('CONDEX',N,K,THETA) is a "counter-example" matrix to a
%   condition estimator. It has order N and scalar parameter
%   THETA (default 100). The matrix, its natural size, and the
%   estimator to which it applies are specified by K as follows:
%      K = 1:   4-by-4,     LINPACK.
%      K = 2:   3-by-3,     LINPACK.
%      K = 3:   arbitrary,  LINPACK (independent of THETA)
%      K = 4:   N >= 4,     LAPACK (RCOND) (default).
%                           It is the inverse of this matrix
%                           that is a counter-example.
%   If N is not equal to the natural size of the matrix, then
%   the matrix is padded out with an identity matrix to order N.

%   Note: In practice, the K = 4 matrix is not usually a counter-example
%   because of the rounding errors in forming it.
%
%   References:
%   [1] A. K. Cline and R. K. Rew, A set of counter-examples to three
%       condition number estimators, SIAM J. Sci. Stat. Comput.,
%       4 (1983), pp. 602-611.
%   [2] N. J. Higham, FORTRAN codes for estimating the one-norm of a
%       real or complex matrix, with applications to condition
%       estimation (Algorithm 674), ACM Trans. Math. Soft., 14 (1988),
%       pp. 381-396.
%   [3] N. J. Higham, Accuracy and Stability of Numerical Algorithms,
%       Second edition, Society for Industrial and Applied Mathematics,
%       Philadelphia, 2002, Chap. 15.
%
%   Nicholas J. Higham
%   Copyright 1984-2005 The MathWorks, Inc.

if isempty(theta), theta = cast(100,classname); end
if isempty(k), k = 4; end

if k == 1    % Cline and Rew (1983), Example B.

   A = [1  -1  -2*theta     0
        0   1     theta  -theta
        0   1   1+theta  -(theta+1)
        0   0   0         theta];

elseif k == 2   % Cline and Rew (1983), Example C.

   A = [1   1-2/theta^2  -2
        0   1/theta      -1/theta
        0   0             1];

elseif k == 3    % Cline and Rew (1983), Example D.

    A = gallery('triw',n,-1,classname)';
    A(n,n) = -1;

elseif k == 4    % Higham (1988), p. 390; Higham (1996), p. 296.

    x = ones(n,3,classname);  %  First col is e
    x(2:n,2) = zeros(n-1,1);  %  Second col is e(1)

    % Third col is special vector b.
    x(:, 3) = (-1).^(0:n-1)' .* ( 1 + (0:n-1)'/(n-1) );

    Q = orth(x);  %  Q*Q' is now the orthogonal projector onto span(e(1),e,b)).
    P = eye(n,classname) - Q*Q';
    A = eye(n,classname) + theta*P;

end

% Pad out with identity as necessary.
m = length(A);
if m < n
   for i=n:-1:m+1, A(i,i) = 1; end
end
