Mercurial > matrix-functions

comparison toolbox/kahan.m @ 2:c124219d7bfa draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 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