--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/toolbox/condex.m	Wed May 06 14:56:53 2015 +0200
@@ -0,0 +1,63 @@
+function A = condex(n, k, theta)
+%CONDEX   `Counterexamples' to matrix condition number estimators.
+%         CONDEX(N, K, THETA) is a `counterexample' matrix to a condition
+%         estimator.  It has order N and scalar parameter THETA (default 100).
+%         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.
+%         The matrix, its natural size, and the estimator to which it applies
+%         are specified by K (default K = 4) as follows:
+%             K = 1:   4-by-4,     LINPACK (RCOND)
+%             K = 2:   3-by-3,     LINPACK (RCOND)
+%             K = 3:   arbitrary,  LINPACK (RCOND) (independent of THETA)
+%             K = 4:   N >= 4,     SONEST (Higham 1988)
+%         (Note that in practice the K = 4 matrix is not usually a
+%          counterexample because of the rounding errors in forming it.)
+
+%         References:
+%         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.
+%         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.
+
+if nargin < 3, theta = 100; end
+if nargin < 2, 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 = triw(n,-1)';
+    A(n,n) = -1;
+
+elseif k == 4    % Higham (1988), p. 390.
+
+    x = ones(n,3);            %  First col is e
+    x(2:n,2) = zeros(n-1,1);  %  Second col is e(1)
+
+    % Third col is special vector b in SONEST
+    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) - Q*Q';
+    A = eye(n) + theta*P;
+
+end
+
+% Pad out with identity as necessary.
+[m, m] = size(A);
+if m < n
+   for i=n:-1:m+1, A(i,i) = 1; end
+end