Mercurial > matrix-functions
diff toolbox/gfpp.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 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolbox/gfpp.m Thu May 07 18:36:24 2015 +0200 @@ -0,0 +1,40 @@ +function A = gfpp(T, c) +%GFPP Matrix giving maximal growth factor for Gaussian elim. with pivoting. +% GFPP(T) is a matrix of order N for which Gaussian elimination +% with partial pivoting yields a growth factor 2^(N-1). +% T is an arbitrary nonsingular upper triangular matrix of order N-1. +% GFPP(T, C) sets all the multipliers to C (0 <= C <= 1) +% and gives growth factor (1+C)^(N-1). +% GFPP(N, C) (a special case) is the same as GFPP(EYE(N-1), C) and +% generates the well-known example of Wilkinson. + +% Reference: +% N.J. Higham and D.J. Higham, Large growth factors in +% Gaussian elimination with pivoting, SIAM J. Matrix Analysis and +% Appl., 10 (1989), pp. 155-164. + +if norm(T-triu(T),1) | any(~diag(T)) + error('First argument must be a nonsingular upper triangular matrix.') +end + +if nargin == 1, c = 1; end + +if c < 0 | c > 1 + error('Second parameter must be a scalar between 0 and 1 inclusive.') +end + +[m, m] = size(T); +if m == 1 % Handle the special case T = scalar + n = T; + m = n-1; + T = eye(n-1); +else + n = m+1; +end + +d = 1+c; +L = eye(n) - c*tril(ones(n), -1); +U = [T (d.^[0:n-2])'; zeros(1,m) d^(n-1)]; +A = L*U; +theta = max(abs(A(:))); +A(:,n) = (theta/norm(A(:,n),inf)) * A(:,n);