Mercurial > matrix-functions
diff toolbox/condex.m @ 2:c124219d7bfa draft
Re-add the 1995 toolbox after noticing the statement in the ~higham/mctoolbox/ webpage.
author | Antonio Pino Robles <data.script93@gmail.com> |
---|---|
date | Thu, 07 May 2015 18:36:24 +0200 |
parents | 8f23314345f4 |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolbox/condex.m Thu May 07 18:36:24 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