Mercurial > matrix-functions
comparison 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 |
comparison
equal
deleted
inserted
replaced
| 1:e471a92d17be | 2:c124219d7bfa |
|---|---|
| 1 function A = condex(n, k, theta) | |
| 2 %CONDEX `Counterexamples' to matrix condition number estimators. | |
| 3 % CONDEX(N, K, THETA) is a `counterexample' matrix to a condition | |
| 4 % estimator. It has order N and scalar parameter THETA (default 100). | |
| 5 % If N is not equal to the `natural' size of the matrix then | |
| 6 % the matrix is padded out with an identity matrix to order N. | |
| 7 % The matrix, its natural size, and the estimator to which it applies | |
| 8 % are specified by K (default K = 4) as follows: | |
| 9 % K = 1: 4-by-4, LINPACK (RCOND) | |
| 10 % K = 2: 3-by-3, LINPACK (RCOND) | |
| 11 % K = 3: arbitrary, LINPACK (RCOND) (independent of THETA) | |
| 12 % K = 4: N >= 4, SONEST (Higham 1988) | |
| 13 % (Note that in practice the K = 4 matrix is not usually a | |
| 14 % counterexample because of the rounding errors in forming it.) | |
| 15 | |
| 16 % References: | |
| 17 % A.K. Cline and R.K. Rew, A set of counter-examples to three | |
| 18 % condition number estimators, SIAM J. Sci. Stat. Comput., | |
| 19 % 4 (1983), pp. 602-611. | |
| 20 % N.J. Higham, FORTRAN codes for estimating the one-norm of a real or | |
| 21 % complex matrix, with applications to condition estimation | |
| 22 % (Algorithm 674), ACM Trans. Math. Soft., 14 (1988), pp. 381-396. | |
| 23 | |
| 24 if nargin < 3, theta = 100; end | |
| 25 if nargin < 2, k = 4; end | |
| 26 | |
| 27 if k == 1 % Cline and Rew (1983), Example B. | |
| 28 | |
| 29 A = [1 -1 -2*theta 0 | |
| 30 0 1 theta -theta | |
| 31 0 1 1+theta -(theta+1) | |
| 32 0 0 0 theta]; | |
| 33 | |
| 34 elseif k == 2 % Cline and Rew (1983), Example C. | |
| 35 | |
| 36 A = [1 1-2/theta^2 -2 | |
| 37 0 1/theta -1/theta | |
| 38 0 0 1]; | |
| 39 | |
| 40 elseif k == 3 % Cline and Rew (1983), Example D. | |
| 41 | |
| 42 A = triw(n,-1)'; | |
| 43 A(n,n) = -1; | |
| 44 | |
| 45 elseif k == 4 % Higham (1988), p. 390. | |
| 46 | |
| 47 x = ones(n,3); % First col is e | |
| 48 x(2:n,2) = zeros(n-1,1); % Second col is e(1) | |
| 49 | |
| 50 % Third col is special vector b in SONEST | |
| 51 x(:, 3) = (-1).^[0:n-1]' .* ( 1 + [0:n-1]'/(n-1) ); | |
| 52 | |
| 53 Q = orth(x); % Q*Q' is now the orthogonal projector onto span(e(1),e,b)). | |
| 54 P = eye(n) - Q*Q'; | |
| 55 A = eye(n) + theta*P; | |
| 56 | |
| 57 end | |
| 58 | |
| 59 % Pad out with identity as necessary. | |
| 60 [m, m] = size(A); | |
| 61 if m < n | |
| 62 for i=n:-1:m+1, A(i,i) = 1; end | |
| 63 end |