Mercurial > matrix-functions

comparison toolbox/gfpp.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 A = gfpp(T, c)
2 %GFPP Matrix giving maximal growth factor for Gaussian elim. with pivoting.
3 % GFPP(T) is a matrix of order N for which Gaussian elimination
4 % with partial pivoting yields a growth factor 2^(N-1).
5 % T is an arbitrary nonsingular upper triangular matrix of order N-1.
6 % GFPP(T, C) sets all the multipliers to C (0 <= C <= 1)
7 % and gives growth factor (1+C)^(N-1).
8 % GFPP(N, C) (a special case) is the same as GFPP(EYE(N-1), C) and
9 % generates the well-known example of Wilkinson.
10
11 % Reference:
12 % N.J. Higham and D.J. Higham, Large growth factors in
13 % Gaussian elimination with pivoting, SIAM J. Matrix Analysis and
14 % Appl., 10 (1989), pp. 155-164.
15
16 if norm(T-triu(T),1) | any(~diag(T))
17 error('First argument must be a nonsingular upper triangular matrix.')
18 end
19
20 if nargin == 1, c = 1; end
21
22 if c < 0 | c > 1
23 error('Second parameter must be a scalar between 0 and 1 inclusive.')
24 end
25
26 [m, m] = size(T);
27 if m == 1 % Handle the special case T = scalar
28 n = T;
29 m = n-1;
30 T = eye(n-1);
31 else
32 n = m+1;
33 end
34
35 d = 1+c;
36 L = eye(n) - c*tril(ones(n), -1);
37 U = [T (d.^[0:n-2])'; zeros(1,m) d^(n-1)];
38 A = L*U;
39 theta = max(abs(A(:)));
40 A(:,n) = (theta/norm(A(:,n),inf)) * A(:,n);