Mercurial > matrix-functions

comparison toolbox/kahan.m @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 function U = kahan(n, theta, pert)
2 %KAHAN Kahan matrix - upper trapezoidal.
3 % KAHAN(N, THETA) is an upper trapezoidal matrix
4 % that has some interesting properties regarding estimation of
5 % condition and rank.
6 % The matrix is N-by-N unless N is a 2-vector, in which case it
7 % is N(1)-by-N(2).
8 % The parameter THETA defaults to 1.2.
9 % The useful range of THETA is 0 < THETA < PI.
10 %
11 % To ensure that the QR factorization with column pivoting does not
12 % interchange columns in the presence of rounding errors, the diagonal
13 % is perturbed by PERT*EPS*diag( [N:-1:1] ).
14 % The default is PERT = 25, which ensures no interchanges for KAHAN(N)
15 % up to at least N = 90 in IEEE arithmetic.
16 % KAHAN(N, THETA, PERT) uses the given value of PERT.
17
18 % The inverse of KAHAN(N, THETA) is known explicitly: see
19 % Higham (1987, p. 588), for example.
20 % The diagonal perturbation was suggested by Christian Bischof.
21 %
22 % References:
23 % W. Kahan, Numerical linear algebra, Canadian Math. Bulletin,
24 % 9 (1966), pp. 757-801.
25 % N.J. Higham, A survey of condition number estimation for
26 % triangular matrices, SIAM Review, 29 (1987), pp. 575-596.
27
28 if nargin < 3, pert = 25; end
29 if nargin < 2, theta = 1.2; end
30
31 r = n(1); % Parameter n specifies dimension: r-by-n.
32 n = n(max(size(n)));
33
34 s = sin(theta);
35 c = cos(theta);
36
37 U = eye(n) - c*triu(ones(n), 1);
38 U = diag(s.^[0:n-1])*U + pert*eps*diag( [n:-1:1] );
39 if r > n
40 U(r,n) = 0; % Extend to an r-by-n matrix.
41 else
42 U = U(1:r,:); % Reduce to an r-by-n matrix.
43 end