Mercurial > matrix-functions

comparison toolbox/triw.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 t = triw(n, alpha, k)
2 %TRIW Upper triangular matrix discussed by Wilkinson and others.
3 % TRIW(N, ALPHA, K) is the upper triangular matrix with ones on
4 % the diagonal and ALPHAs on the first K >= 0 superdiagonals.
5 % N may be a 2-vector, in which case the matrix is N(1)-by-N(2) and
6 % upper trapezoidal.
7 % Defaults: ALPHA = -1,
8 % K = N - 1 (full upper triangle).
9 % TRIW(N) is a matrix discussed by Kahan, Golub and Wilkinson.
10 %
11 % Ostrowski (1954) shows that
12 % COND(TRIW(N,2)) = COT(PI/(4*N))^2,
13 % and for large ABS(ALPHA),
14 % COND(TRIW(N,ALPHA)) is approximately ABS(ALPHA)^N*SIN(PI/(4*N-2)).
15 %
16 % Adding -2^(2-N) to the (N,1) element makes TRIW(N) singular,
17 % as does adding -2^(1-N) to all elements in the first column.
18
19 % References:
20 % G.H. Golub and J.H. Wilkinson, Ill-conditioned eigensystems and the
21 % computation of the Jordan canonical form, SIAM Review,
22 % 18(4), 1976, pp. 578-619.
23 % W. Kahan, Numerical linear algebra, Canadian Math. Bulletin,
24 % 9 (1966), pp. 757-801.
25 % A.M. Ostrowski, On the spectrum of a one-parametric family of
26 % matrices, J. Reine Angew. Math., 193 (3/4), 1954, pp. 143-160.
27 % J.H. Wilkinson, Singular-value decomposition---basic aspects,
28 % in D.A.H. Jacobs, ed., Numerical Software---Needs and Availability,
29 % Academic Press, London, 1978, pp. 109-135.
30
31 m = n(1); % Parameter n specifies dimension: m-by-n.
32 n = n(max(size(n)));
33
34 if nargin < 3, k = n-1; end
35 if nargin < 2, alpha = -1; end
36
37 if max(size(alpha)) ~= 1
38 error('Second argument must be a scalar.')
39 end
40
41 t = tril( eye(m,n) + alpha*triu(ones(m,n), 1), k);