# To unbundle, sh this file echo bandred.m 1>&2 cat >bandred.m <<'End of bandred.m' function A = bandred(A, kl, ku) %BANDRED B = BANDRED(A, KL, KU) is a matrix orthogonally equivalent to A % with lower bandwidth KL and upper bandwidth KU % (i.e. B(i,j) = 0 if i>j+KL or j>i+KU). % The reduction is performed using Householder transformations. % If KU is omitted it defaults to KL. % This is a `standard' reduction. Cf. reduction to bidiagonal form % prior to computing the SVD. This code is a little wasteful in that % it computes certain elements which are immediately set to zero! if nargin == 2, ku = kl; end if kl == 0 & ku == 0 error('You''ve asked for a diagonal matrix. In that case use the SVD!') end % Check for special case where order of left/right transformations matters. % Easiest approach is to work on the transpose, flipping back at the end. flip = 0; if ku == 0 A = A'; temp = kl; kl = ku; ku = temp; flip = 1; end [m, n] = size(A); for j=1 : min( min(m,n), max(m-kl-1,n-ku-1) ) if j+kl+1 <= m [v, beta] = house(A(j+kl:m,j)); temp = A(j+kl:m,j:n); A(j+kl:m,j:n) = temp - beta*v*(v'*temp); A(j+kl+1:m,j) = zeros(m-j-kl,1); end if j+ku+1 <= n [v, beta] = house(A(j,j+ku:n)'); temp = A(j:m,j+ku:n); A(j:m,j+ku:n) = temp - beta*(temp*v)*v'; A(j,j+ku+1:n) = zeros(1,n-j-ku); end end if flip A = A'; end End of bandred.m echo cauchy.m 1>&2 cat >cauchy.m <<'End of cauchy.m' function C = cauchy(x, y) %CAUCHY C = CAUCHY(X,Y), where X,Y are N-vectors, is the N-by-N matrix % with C(i,j) = 1/(X(i)+Y(j)). By default, Y = X. % Special case: if X is a scalar CAUCHY(X) is the same as CAUCHY(1:X). % Explicit formulas are known for DET(C) (which is nonzero if X and Y % both have distinct elements) and the elements of INV(C). % Reference: D.E. Knuth, The Art of Computer Programming, Volume 1, % Fundamental Algorithms, Second Edition, Addison-Wesley, Reading, % Massachusetts, 1973, p. 36. n = max(size(x)); % Handle scalar x. if n == 1 n = x; x = 1:n; end if nargin == 1, y = x; end x = x(:); y = y(:); % Ensure x and y are column vectors. if any(size(x) ~= size(y)) error('Parameter vectors must be of same dimension.') end C = x*ones(1,n) + ones(n,1)*y.'; C = ones(n) ./ C; End of cauchy.m echo chebspec.m 1>&2 cat >chebspec.m <<'End of chebspec.m' function C = chebspec(n, k) %CHEBSPEC C = CHEBSPEC(N, K) is a Chebyshev spectral differentiation % matrix of order N. K = 0 (the default) or 1. % For K=0 (`no boundary conditions'), C is nilpotent, with % C^N=0 and it has the null vector ONES(N,1). C is similar to a % Jordan block of size N with eigenvalue zero. % For K=1, C is nonsingular and well-conditioned, and its eigenvalues % have negative real parts. % For both K, the computed eigenvector matrix X from EIG is % ill-conditioned (MESH(X) is interesting). % References: % C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, Spectral % Methods in Fluid Dynamics, Springer-Verlag, Berlin, 1987; p. 69. % L.N. Trefethen, M.R. Trummer, An instability phenomenon in spectral % methods, SIAM J. Numer. Anal., 24 (1987), pp. 1008-1023. % D. Funaro, Computing the inverse of the Chebyshev collocation % derivative, SIAM J. Sci. Stat. Comput., 9 (1988), pp. 1050-1057. if nargin == 1, k = 0; end % k = 1 case obtained from k = 0 case with one bigger n. if k == 1, n = n + 1; end n = n-1; C = zeros(n+1); one = ones(n+1,1); x = cos( (0:n)' * (pi/n) ); d = ones(n+1,1); d(1) = 2; d(n+1) = 2; C = (d * (one./d)') ./ (x*one'-one*x' +eye(C)); % eye(C) avoids div by zero % Now fix diagonal and signs. C(1,1) = (2*n^2+1)/6; for i=2:n+1 if rem(i,2) == 0 C(:,i) = -C(:,i); C(i,:) = -C(i,:); end if i < n+1 C(i,i) = -x(i)/(2*(1-x(i)^2)); else C(n+1,n+1) = -C(1,1); end end if k == 1 C = C(2:n+1,2:n+1); end End of chebspec.m echo chebvand.m 1>&2 cat >chebvand.m <<'End of chebvand.m' function C = chebvand(m,p) %CHEBVAND Vandermonde-like matrix for the Chebyshev polynomials. % C = CHEBVAND(P), where P is a vector, produces the (primal) % Chebyshev Vandermonde matrix based on the points P, % i.e. C(i,j) = T_{i-1}(P(j)), where T_{i-1} is the Chebyshev % polynomial of degree i-1. % CHEBVAND(M,P) is a rectangular version of CHEBVAND(P) with M rows. % Special case: If P is a scalar then P equally spaced points on % [0,1] are used. % Reference: N.J. Higham, Stability analysis of algorithms for solving % confluent Vandermonde-like systems, to appear in SIAM J. Matrix Anal. Appl. if nargin == 1, p = m; end n = max(size(p)); % Handle scalar p. if n == 1 n = p; p = seqa(0,1,n); end p = p(:).'; % Ensure p is a row vector. C = ones(m,n); if m == 1, return, end C(2,:) = p; % Use Chebyshev polynomial recurrence. for i=3:m C(i,:) = 2.*p.*C(i-1,:) - C(i-2,:); end End of chebvand.m echo chop.m 1>&2 cat >chop.m <<'End of chop.m' function c = chop(x, t) %CHOP Rounds matrix x to t significant binary places. % Default is t=23, corresponding to IEEE single precision. if nargin < 2, t = 23; end [m, n] = size(x); % Use the representation: % x(i,j) = 2^e(i,j) * .d(1)d(2)...d(s) * sign(x(i,j)) % On the next line `+(x==0)' avoids passing a zero argument to LOG, which % would cause a warning message to be generated. y = abs(x) + (x==0); e = floor(log(y)./log(2) + 1); p = (2.*ones(m,n)).^(t-e); c = round(p.*x)./p; End of chop.m echo chow.m 1>&2 cat >chow.m <<'End of chow.m' function A = chow(n, alpha, delta) %CHOW A = CHOW(N, ALPHA, DELTA) is a Toeplitz lower Hessenberg matrix % A = H(ALPHA) + DELTA*EYE, where H(i,j) = ALPHA^(i-j+1). % H(ALPHA) has p=FLOOR((N+1)/2) zero eigenvalues, the rest being % 4*ALPHA*COS( k*PI/(N+2) )^2, k=1:N-p. % Defaults: ALPHA = 1, DELTA = 0. % References: % T.S. Chow, A class of Hessenberg matrices with known % eigenvalues and inverses, SIAM Review, 11 (1969), pp. 391-395. % G. Fairweather, On the eigenvalues and eigenvectors of a class of % Hessenberg matrices, SIAM Review, 13 (1971), pp. 220-221. if nargin < 3, delta = 0; end if nargin < 2, alpha = 1; end A = toeplitz( alpha.^(1:n)', [alpha 1 zeros(1,n-2)] ) + delta*eye(n); End of chow.m echo circul.m 1>&2 cat >circul.m <<'End of circul.m' function C = circul(v) %CIRCUL C = CIRCUL(V) is the circulant matrix whose first row is V. % (A circulant matrix has the property that each row is obtained % from the previous one by cyclically permuting the entries one step % forward; it is a special Toeplitz matrix in which the diagonals % `wrap round'.) % Special case: if V is a scalar then C = CIRCUL(1:V). % The eigensystem of C (N-by-N) is known explicitly. If t is an Nth % root of unity, then the inner product of V with W = [1 t t^2 ... t^N] % is an eigenvalue of C, and W(N:-1:1) is an eigenvector of C. % Reference: P.J. Davis, Circulant Matrices, John Wiley, 1977. n = max(size(v)); if n == 1 n = v; v = 1:n; end v = v(:).'; % Make sure v is a row vector. C = toeplitz( [ v(1) v(n:-1:2) ], v ); End of circul.m echo clement.m 1>&2 cat >clement.m <<'End of clement.m' function A = clement(n, k) %CLEMENT CLEMENT(N, K) is a tridiagonal matrix with zero diagonal entries % and known eigenvalues. It is singular if N is odd. About 64 % percent of the entries of the inverse are zero. The eigenvalues % are plus and minus the numbers N-1, N-3, N-5, ..., (1 or 0). % For K=0 (the default) the matrix is nonsymmetric, while for K=1 it % is symmetric. % Similar properties hold for TRIDIAG(X,Y,Z) where Y = ZEROS(N,1). % The eigenvalues still come in plus/minus pairs but they are not % known explicitly. % Reference: P.A. Clement, A class of triple-diagonal matrices for test % purposes, SIAM Review, 1 (1959), pp. 50-52. if nargin == 1, k = 0; end n = n-1; x = n:-1:1; z = 1:n; if k == 0 A = diag(x, -1) + diag(z, 1); else y = sqrt(x.*z); A = diag(y, -1) + diag(y, 1); end End of clement.m echo comp.m 1>&2 cat >comp.m <<'End of comp.m' function C = comp(A, k) %COMP Comparison matrices. % COMP(A) is DIAG(B) - TRIL(B,-1) - TRIU(B,1), where B = ABS(A). % COMP(A, 1) is A with each diagonal element replaced by its % absolute value, and each off-diagonal element replaced by minus % the absolute value of the largest element in absolute value in % its row. However, if A is triangular COMP(A, 1) is too. % COMP(A, 0) is the same as COMP(A). % COMP(A) is often denoted by M(A) in the literature. % Reference (e.g.): N.J. Higham, A survey of condition number estimation for % triangular matrices, SIAM Review, 29 (1987), pp. 575-596. if nargin == 1, k = 0; end [m, n] = size(A); p = min(m, n); if k == 0 % This code uses less temporary storage than the `high level' definition above. C = -abs(A); for j=1:p C(j,j) = abs(A(j,j)); end elseif k == 1 C = A'; for j=1:p C(k,k) = 0; end mx = max(abs(C)); C = -mx'*ones(1,n); for j=1:p C(j,j) = abs(A(j,j)); end if all( A == tril(A) ), C = tril(C); end if all( A == triu(A) ), C = triu(C); end end End of comp.m echo compan.m 1>&2 cat >compan.m <<'End of compan.m' function A = compan(p) %COMPAN COMPAN(P), where P is an (n+1)-vector, is the n-by-n companion matrix % whose first row is -P(2:n+1)/P(1). % Special case: if P is a scalar then COMPAN(P) is the P-by-P matrix % COMPAN(1:P+1). % Reference: J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford % University Press, 1965, p. 12. n = max(size(p)); % Handle scalar p. if n == 1 n = p+1; p = 1:n; end p = p(:)'; % Ensure p is a row vector. % Construct matrix of order n-1. if n == 2 A = 1; else A = diag(ones(1,n-2),-1); A(1,:) = -p(2:n)/p(1); end End of compan.m echo condex.m 1>&2 cat >condex.m <<'End of condex.m' function A = condex(n, k, theta) %CONDEX CONDEX(N, K, THETA) is a `counter-example' matrix to a condition % estimator. It has order N and scalar parameter THETA (default 100). % If N is not equal to the `natural' size of the matrix then % the matrix is padded out with an identity matrix to order N. % The matrix, its natural size, and the estimator to which it applies % are specified by K (default K = 4) as follows: % K=1: 4-by-4, LINPACK (RCOND) % K=2: 3-by-3, LINPACK (RCOND) % K=3: arbitrary, LINPACK (RCOND) % K=4: N >= 4, SONEST (Higham 1988) % (Note that in practice the K=4 matrix is not usually a % counter-example because of the rounding errors in forming it.) % References: % A.K. Cline and R.K. Rew, A set of counter-examples to three condition number % estimators, SIAM J. Sci. Stat. Comput., 4 (1983), pp. 602-611. % N.J. Higham, FORTRAN codes for estimating the one-norm of a real or complex % matrix, with applications to condition estimation, ACM Trans. Math. % Soft., 14 (1988), pp. 381-396. if nargin < 3, theta = 100; end if nargin < 2, k = 4; end if k == 1 % Cline and Rew (1983), Example B. A = [1 -1 -2*theta 0 0 1 theta -theta 0 1 1+theta -(theta+1) 0 0 0 theta]; elseif k == 2 % Cline and Rew (1983), Example C. A = [1 1-1/theta^2 -2 0 1/theta -1/theta 0 0 1]; elseif k == 3 % Cline and Rew (1983), Example D. A = triw(n,-1)'; A(n,n) = -1; elseif k == 4 % Higham (1988), p. 390. x = ones(n,3); % First col is e x(2:n,2) = zeros(n-1,1); % Second col is e(1) % Third col is special vector b in SONEST x(:, 3) = (-1).^[0:n-1]' .* ( 1 + [0:n-1]'/(n-1) ); Q = orth(x); % Q*Q' is now the orthogonal projector onto span(e(1),e,b)). P = eye(n) - Q*Q'; A = eye(n) + theta*P; end % Pad out with identity as necessary. [m, m] = size(A); if m < n for i=n:-1:m+1, A(i,i) = 1; end end End of condex.m echo cycol.m 1>&2 cat >cycol.m <<'End of cycol.m' function A = cycol(n, k) %CYCOL A = CYCOL([M N], K) is an M-by-N matrix of the form A = B(1:M,1:N) % where B = [C C C...] and C=RAND(M,K). Thus A's columns repeat % cyclically, and A has rank at most K. K need not divide N. % K defaults to ROUND(N/4). % CYCOL(N,K), where N is a scalar, is the same as CYCOL([N N], K). % This type of matrix can lead to underflow problems for Gaussian % elimination: see NA Digest Volume 89, Issue 3 (January 22, 1989). m = n(1); % Parameter n specifies dimension: m-by-n. n = n(max(size(n))); if nargin < 2, k = max(round(n/4),1); end A = rand(m, k); for i=2:ceil(n/k) A = [A A(:,1:k)]; end A = A(:, 1:n); End of cycol.m echo dingdong.m 1>&2 cat >dingdong.m <<'End of dingdong.m' function A = dingdong(n) %DINGDONG A = DINGDONG(N) is the symmetric n-by-n matrix with % A(i,j) = 0.5/(n-i-j+1.5). % The eigenvalues of A cluster around PI/2 and -PI/2. % Invented by F.N. Ris. % Reference: J.C. Nash, Compact Numerical Methods for Computers: Linear % Algebra and Function Minimisation, Adam Hilger, Bristol, and John Wiley, % New York, 1979, p. 210. p= -2*(1:n) + (n+1.5); A = cauchy(p); End of dingdong.m echo dorr.m 1>&2 cat >dorr.m <<'End of dorr.m' function [c, d, e] = dorr(n, theta) %DORR [C, D, E] = DORR(N, THETA) returns the vectors defining a diagonally % dominant, tridiagonal M-matrix which is ill-conditioned for small % values of the parameter THETA>=0. If only one output parameter is % supplied then C = TRIDIAG(C,D,E), i.e. the full matrix is returned. % The columns of INV(C) vary greatly in norm. THETA defaults to 0.01. % Reference: F.W. Dorr, An example of ill-conditioning in the numerical % solution of singular perturbation problems, Math. Comp., 25 (1971), % pp. 271-283. if nargin < 2, theta = 0.01; end c = zeros(n,1); e = c; d = c; % All length n for convenience. Make c,e of length n-1 later. h = 1/(n+1); m = floor( (n+1)/2 ); term = theta/h^2; i = (1:m)'; c(i) = -term*ones(i); e(i) = c(i) - (0.5-i*h)/h; d(i) = -(c(i) + e(i)); i = (m+1:n)'; e(i) = -term*ones(n-m,1); c(i) = e(i) + (0.5-i*h)/h; d(i) = -(c(i) + e(i)); c = c(2:n); e = e(1:n-1); if nargout <= 1 c = tridiag(c, d, e); end End of dorr.m echo dramadah.m 1>&2 cat >dramadah.m <<'End of dramadah.m' function A = dramadah(n, k) %DRAMADAH An anti-Hadamard matrix A is a matrix with elements 0 or 1 for % which MU(A) := NORM(INV(A),'FRO') is maximal. % A = DRAMADAH(N, K) is an N-by-N (0,1) matrix for which MU(A) is % relatively large, although not necessarily maximal. % Available types (the default is K=1): % K = 1: A is Toeplitz, with ABS(DET(A)) = 1, and MU(A) > c(1.75)^N, % where c is a constant. % K = 2: A is upper triangular and Toeplitz. % The inverses of both types have integer entries. % Reference: R.L. Graham and N.J.A. Sloane, Anti-Hadamard matrices, % Linear Algebra and Appl., 62 (1984), pp. 113-137. if nargin < 2, k = 1; end if k == 1 % Toeplitz c = ones(n,1); for i=2:4:n m = min(1,n-i); c(i:i+m) = zeros(m+1,1); end r = zeros(n,1); r(1:4) = [1 1 0 1]; if n < 4, r = r(1:n); end A = toeplitz(c,r); else % Upper triangular and Toeplitz c = zeros(n,1); c(1) = 1; r = ones(n,1); for i = 3:2:n r(i) = 0; end A = toeplitz(c,r); end End of dramadah.m echo fiedler.m 1>&2 cat >fiedler.m <<'End of fiedler.m' function A = fiedler(c) %FIEDLER A = FIEDLER(C), where C is an n-vector, is the n-by-n symmetric % matrix with elements ABS(C(i)-C(j)). % Special case: if C is a scalar, then A = FIEDLER(1:C) % (i.e. A(i,j) = ABS(i-j)). % Properties: % INV(A) is tridiagonal except for nonzero (1,n) and (n,1) elements. % DET(A) = (-1)^(n-1)*2^(n-2)*PROD( C(2:n)-C(1:n-1) ) % FIEDLER(N) has a dominant positive eigenvalue and all the other % eigenvalues are negative (Szego 1936). % Reference: G. Szego, Solution to problem 3705, Amer. Math. Monthly, % 43 (1936), pp. 246-259. n = max(size(c)); % Handle scalar c. if n == 1 n = c; c = 1:n; end c = c(:).'; % Ensure c is a row vector. A = ones(n,1)*c; A = abs(A-A.'); % NB. array transpose End of fiedler.m echo forsythe.m 1>&2 cat >forsythe.m <<'End of forsythe.m' function A = forsythe(n, alpha, lambda) %FORSYTHE FORSYTHE(N, ALPHA, LAMBDA) is the N-by-N matrix equal to % JORDAN(N, LAMBDA) except it has an ALPHA in the (N,1) position. % It has the characteristic equation % DET(A-t*EYE) = (LAMBDA-t)^N - (-1)^N ALPHA. % ALPHA defaults to SQRT(EPS) and LAMBDA to 0. if nargin < 2, alpha = sqrt(eps); end if nargin < 3, lambda = 0; end A = jordan(n, lambda); A(n,1) = alpha; End of forsythe.m echo frank.m 1>&2 cat >frank.m <<'End of frank.m' function F = frank(n, k) %FRANK F = FRANK(N, K) is the Frank matrix of order n. It is upper % Hessenberg with determinant 1. K = 0 is the default; if K = 1 the % elements are reflected about the anti-diagonal (1,n)--(n,1). % The eigenvalues of F may be obtained in terms of the zeros of the % Hermite polynomials. They are positive and occur in reciprocal % pairs. Thus if N is odd, 1 is an eigenvalue. % F has floor(n/2) ill-conditioned eigenvalues---the smaller ones. % For large N, DET(FRANK(N)') comes out far from 1---see Frank (1958) % for an explanation. % References: % W.L. Frank, Computing eigenvalues of complex matrices by determinant % evaluation and by methods of Danilewski and Wielandt, J. Soc. Indust. % Appl. Math., 6 (1958), pp. 378-392 (see pp. 385, 388). % J.H. Wilkinson, Error analysis of floating-point computation, % Numer. Math., 2 (1960), pp. 319-340 (section 8). % J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University % Press, 1965 (pp. 92-93). % G.H. Golub and J.H. Wilkinson, Ill-conditioned eigensystems and the % computation of the Jordan canonical form, SIAM Review, 18 (1976), % pp. 578-619 (section 13). % The next two references give details of the eigensystem: % P.J. Eberlein, A note on the matrices denoted by \$B_n\$, SIAM J. Appl. % Math., 20 (1971), pp. 87-92. % J.M. Varah, A generalization of the Frank matrix, SIAM J. Sci. Stat. % Comput., 7 (1986), pp. 835-839. if nargin == 1, k = 0; end F = min( ones(n,1)*(1:n), (1:n)'*ones(1,n) ); % Take upper Hessenberg part. F = triu(F,-1); if k == 0 p=n:-1:1; F = F(p,p)'; end End of frank.m echo fv.m 1>&2 cat >fv.m <<'End of fv.m' function [f, e] = fv(b, nk, thmax, noplot) %FV Field of values (or numerical range). % FV(A, NK, THMAX); evaluates and plots the field of values of the NK % largest leading principal submatrices of A, using THMAX equally spaced % angles in the complex plane. The defaults are NK = 1 and THMAX = 20. % The eigenvalues of A are displayed as 'x'. % Alternative usage: [F, E] = FV(A, NK, THMAX, 1) suppresses the plot % and returns the field of values plot data in F, with A's eigenvalues % in E. % Theory: % Field of values FV(A) = set of all Rayleigh quotients. FV(A) is a % convex set containing the eigenvalues of A. When A is normal FV(A) is % the convex hull of the eigenvalues of A (but not vice versa). % z = x'Ax/(x'x), z' = x'A'x/(x'x) % => REAL(z) = x'Hx/(x'x), H = (A+A')/2 % so MIN(EIG(H)) <= REAL(z) <= MAX(EIG(H)) % with equality for x = corresponding eigenvectors of H. For these x, % RQ(A,x) is on the boundary of FV(A). % Reference: A.S. Householder, The Theory of Matrices in Numerical % Analysis, Blaisdell, New York, 1964, section 3.3. % Original version by A. Ruhe. Modified by N.J. Higham. % The plot is `incorrect' for skew-symmetric matrices! if nargin < 2, nk = 1; end if nargin < 3, thmax = 20; end iu = sqrt(-1); [n, p] = size(b); f = []; z = zeros(2*thmax+1,1); for m = 1:nk ns = n+1-m; a = b(1:ns, 1:ns); for i = 0:thmax th = i/thmax*pi; ath = exp(iu*th)*a; % Rotate A through angle th. h = 0.5*(ath + ath'); % Hermitian part of rotated A. [x, d] = eig(h); [lmbh, k] = sort(real(diag(d))); z(1+i) = rq(a,x(:,k(1))); % RQ's of A corr. to eigenvalues of H z(1+i+thmax) = rq(a,x(:,k(ns))); % with smallest/largest real part. end f = [f; z]; end if nargin < 4 axis('square') if norm(imag(f),1) == 0 % In case of real f (e.g. Hermitian A) choose [-1,1] for y-range. axis([ min(f) max(f) -1 1]); end plot(real(f), imag(f),'w') % Plot the field of values hold on e = eig(b); plot(real(e), imag(e), 'xg') % Plot the eigenvalues too. axis('normal'), hold off if norm(imag(f),1) == 0 % Return to auto-scaling axis; end end End of fv.m echo gallery.m 1>&2 cat >gallery.m <<'End of gallery.m' function [A, e] = gallery(n) %GALLERY Famous, and not so famous, test matrices. % A = GALLERY(N) is an N-by-N matrix with some special property. % Only the following values of N are currently available: % N = 3 is badly conditioned. % N = 4 is the Wilson matrix. Symmetric pos def, integer inverse. % N = 5 is an interesting ei'value problem: defective, nilpotent. % N = 8 is the Rosser matrix, a classic symmetric eigenvalue problem. % [A, e] = GALLERY(8) returns the exact eigenvalues in e. % N = 21 is Wilkinson's tridiagonal W21+, another eigenvalue problem. % Original version supplied with MATLAB. Modified by N.J. Higham. % References: % J.R. Westlake, A Handbook of Numerical Matrix Inversion and Solution % of Linear Equations, John Wiley, New York, 1968. % J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University % Press, 1965. if n == 3 A = [ -149 -50 -154 537 180 546 -27 -9 -25 ]; elseif n == 4 A = [10 7 8 7 7 5 6 5 8 6 10 9 7 5 9 10]; elseif n == 5 % disp('Try to find the EXACT eigenvalues and eigenvectors.') % Matrix devised by Cleve Moler. Its Jordan form has just one block, with % eigenvalue zero. Proof: A^k is nonzero for k<5, zero for k=5. % TRACE(A)=0. No simple form for null vector. A = [ -9 11 -21 63 -252 70 -69 141 -421 1684 -575 575 -1149 3451 -13801 3891 -3891 7782 -23345 93365 1024 -1024 2048 -6144 24572 ]; elseif n == 8 A = [ 611. 196. -192. 407. -8. -52. -49. 29. 196. 899. 113. -192. -71. -43. -8. -44. -192. 113. 899. 196. 61. 49. 8. 52. 407. -192. 196. 611. 8. 44. 59. -23. -8. -71. 61. 8. 411. -599. 208. 208. -52. -43. 49. 44. -599. 411. 208. 208. -49. -8. 8. 59. 208. 208. 99. -911. 29. -44. 52. -23. 208. 208. -911. 99. ]; % Exact eigenvalues from Westlake (1968), p.150 (ei'vectors given too): a = sqrt(10405); b = sqrt(26); e = [-10*a, 0, 510-100*b, 1000, 1000, 510+100*b, ... 1020, 10*a]'; elseif n == 21 % W21+, Wilkinson (1965), p.308. E = diag(ones(n-1,1),1); m = (n-1)/2; A = diag(abs(-m:m)) + E + E'; else error('Sorry, that value of N is not available.') end End of gallery.m echo gear.m 1>&2 cat >gear.m <<'End of gear.m' function A = gear(n, i, j) %GEAR A = GEAR(N,I,J) is the N-by-N matrix with ones on the sub- and % super-diagonals, SIGN(I) in the (1,ABS(I)) position, SIGN(J) % in the (N,N+1-ABS(J)) position, and zeros everywhere else. % Defaults: I=N, j=-N. % All eigenvalues are of the form 2*COS(a) and the eigenvectors % are of the form [SIN(w+a), SIN(w+2a), ..., SIN(w+Na)]. % The values of a and w are given in the reference below. % A can have double and triple eigenvalues and can be defective. % GEAR(N) is singular. % Reference: C.W. Gear, A simple set of test matrices for eigenvalue % programs, Math. Comp., 23 (1969), pp. 119-125. if nargin == 1, i = n; j = -n; end if ~(i~=0 & abs(i)<=n & j~=0 & abs(j)<=n) error('Invalid I and J parameters') end A = diag(ones(n-1,1),-1) + diag(ones(n-1,1),1); A(1, abs(i)) = sign(i); A(n, n+1-abs(j)) = sign(j); End of gear.m echo gfpp.m 1>&2 cat >gfpp.m <<'End of gfpp.m' function A = gfpp(T, c) %GFPP 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 T ~= triu(T) | 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); End of gfpp.m echo hadamard.m 1>&2 cat >hadamard.m <<'End of hadamard.m' function H = hadamard(n) %HADAMARD HADAMARD(N) is a Hadamard matrix of order N. % An N-by-N Hadamard matrix exists only if N=2 or REM(N,4)=0. % This function handles only the cases where N, N/12 or N/20 % is a power of 2. % This is an expanded version of a routine supplied with MATLAB. % Reference: S.W. Golomb and L.D. Baumert, The search for Hadamard % matrices, Amer. Math. Monthly, 70 (1963) pp. 12-17. if n ~= 2 & rem(n,4) ~= 0 error('For an NxN Hadamard matrix to exist, N must be 2 or a multiple of 4.') end k2 = log(n)/log(2); k12 = log(n/12)/log(2); k20 = log(n/20)/log(2); if k2 == fix(k2) H = [1 1 1 -1]; k = k2-1; elseif k12 == fix(k12) H = toeplitz( [-1 -1 1 -1 -1 -1 1 1 1 -1 1 ], ... [-1 1 -1 1 1 1 -1 -1 -1 1 -1] ); H = [ones(1,12); ones(11,1) H]; k = k12; elseif k20 == fix(k20) H = hankel( [ -1 -1 1 1 -1 -1 -1 -1 1 -1 1 -1 1 1 1 1 -1 -1 1], ... [1 -1 -1 1 1 -1 -1 -1 -1 1 -1 1 -1 1 1 1 1 -1 -1] ); H = [ones(1,20); ones(19,1) H]; k = k20; else error('Sorry, can''t handle that value of N.') end % Kronecker product construction. for i=1:k H = [H H H -H]; end End of hadamard.m echo hanowa.m 1>&2 cat >hanowa.m <<'End of hanowa.m' function A = hanowa(n, d) %HANOWA HANOWA(N, d) is the N-by-N block 2x2 matrix (thus N must be even) % [d*EYE(N) -DIAG(1:N) % DIAG(1:N) d*EYE(N)] % It has complex eigenvalues lambda(k) = d +/- k*i (1 <= k <= n/2). % Parameter d defaults to -1. % Reference: E. Hairer, S.P. Norsett and G. Wanner, Solving Ordinary % Differential Equations I: Nonstiff Problems, Springer-Verlag, % Berlin, 1987. (pp. 86-87) if nargin == 1, d = -1; end m = n/2; if round(m) ~= m error('N must be even.') end A = [ d*eye(m) -diag(1:m) diag(1:m) d*eye(m)]; End of hanowa.m echo hilb.m 1>&2 cat >hilb.m <<'End of hilb.m' function H = hilb(n) %HILB HILB(N) is the N-by-N Hilbert matrix with elements 1/(i+j-1), % which is a famous example of a badly conditioned matrix. % See INVHILB (standard MATLAB routine) for the exact inverse, which % has integer entries. % HILB(N) is symmetric positive definite and totally positive. % This routine is shorter and faster than the one supplied with MATLAB. % References: D.E. Knuth, The Art of Computer Programming, % Volume 1, Fundamental Algorithms, Second Edition, Addison-Wesley, % Reading, Massachusetts, 1973, p. 37. % M.-D. Choi, Tricks or treats with the Hilbert matrix, Amer. Math. % Monthly, 90 (1983), pp. 301-312. if n == 1 H = 1; else H = cauchy( (1:n) - .5); end End of hilb.m echo house.m 1>&2 cat >house.m <<'End of house.m' function [v, beta] = house(x) %HOUSE Householder matrix. % If [v, beta] = HOUSE(x) then H = EYE - beta*v*v' is a Householder % matrix such that Hx = -sign(x(1))*norm(x)*e_1. % NB: x is assumed to be nonzero. x can be real or complex. % sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0). % Theory: (textbook references Golub & Van Loan 1983, 38-43; % Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50). % Hx = y: (I - beta*v*v')x = -s*e_1. % Must have |s| = norm(x), v = x+s*e_1, and % x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)). % So take s = sign(x(1))*norm(x) (which avoids cancellation). % v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2 % = 2*norm(x)*(norm(x) + |x(1)|). s = norm(x) * (sign(x(1)) + (x(1)==0)); % Modification for sign(0)=1. v = x; v(1) = v(1) + s; beta = 1/(s'*v(1)); % NB the conjugated s. % beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real. % But beta as above can be non-real (due to rounding) only when x is complex. End of house.m echo invol.m 1>&2 cat >invol.m <<'End of invol.m' function A = invol(n) %INVOL A = INVOL(N) is an N-by-N involutory (A*A=EYE(N)) and % ill-conditioned matrix. % It is a diagonally scaled version of HILB(N). % NB. B=(EYE(N)-A)/2 and B=(EYE(N)+A)/2 are idempotent (B*B=B). % Reference: A.S. Householder and J.A. Carpenter, The singular values % of involutory and of idempotent matrices, Numer. Math. 5 (1963), % pp. 234-237. A = hilb(n); d = -n; A(:, 1) = d*A(:, 1); for i = 1:n-1 d = -(n+i)*(n-i)*d/(i*i); A(i+1, :) = d*A(i+1, :); end End of invol.m echo ipjfact.m 1>&2 cat >ipjfact.m <<'End of ipjfact.m' function A = ipjfact(n, k) %IPJFACT A = IPJFACT(N, K) is the matrix with % a(i,j) = (i+j)! (K=0, default) % a(i,j) = 1/(i+j)! (K=1) % Both are Hankel matrices. % Suggested by P.R. Graves-Morris. if nargin == 1, k = 0; end c = cumprod(2:n+1); d = cumprod(n+1:2*n) * c(n-1); A = hankel(c, d); if k == 1 A = ones(A)./A; end End of ipjfact.m echo jordan.m 1>&2 cat >jordan.m <<'End of jordan.m' function J = jordan(n, lambda) %JORDAN JORDAN(N, LAMBDA) is the N-by-N Jordan block with eigenvalue LAMBDA. % LAMBDA = 1 is the default. if nargin == 1, lambda = 1; end J = lambda*eye(n) + diag(ones(n-1,1),1); End of jordan.m echo kahan.m 1>&2 cat >kahan.m <<'End of kahan.m' function U = kahan(n, theta) %KAHAN KAHAN(N, THETA) is an upper trapezoidal matrix % involving a parameter THETA, which has some interesting % properties regarding estimation of condition and rank. % The matrix is N-by-N unless N is a 2-vector, in which case it % is N(1)-by-N(2). THETA defaults to 0.25. % The inverse of KAHAN(N, THETA) is known explicitly: see % Higham (1987, p. 588), for example. % References: % W. Kahan, Numerical linear algebra, Canadian Math. Bulletin, % 9 (1966), pp. 757-801. % N.J. Higham, A survey of condition number estimation for % triangular matrices, SIAM Review, 29 (1987), pp. 575-596. if nargin < 2, theta = 0.25; end r = n(1); % Parameter n specifies dimension: r-by-n. n = n(max(size(n))); s = sin(theta); c = cos(theta); U = eye(n) - c*triu(ones(n), 1); U = diag(s.^[1:n])*U; if r > n U(r,n) = 0; % Extend to an r-by-n matrix. else U = U(1:r,:); % Reduce to an r-by-n matrix. end End of kahan.m echo kms.m 1>&2 cat >kms.m <<'End of kms.m' function A = kms(n, rho) %KMS A = KMS(N, RHO) is the N-by-N Kac-Murdock-Szego Toeplitz matrix with % A(i,j) = RHO^(ABS((i-j))) (for real RHO). % If RHO is complex, then the same formula holds except that elements % below the diagonal are conjugated. % RHO defaults to 0.5. % Properties: % A has an LDL' factorization with % L = INV(TRIW(N,-RHO,1)'), % D(i,i) = (1-ABS(RHO)^2)*EYE(N) except D(1,1) = 1. % A is positive definite if and only if 0 < ABS(RHO) < 1. % INV(A) is tridiagonal. % Reference: W.F. Trench, Numerical solution of the eigenvalue problem % for Hermitian Toeplitz matrices, SIAM J. Matrix Analysis and Appl., % 10 (1989), pp. 135-146 (and see the references therein). if nargin < 2, rho = 0.5; end A = (1:n)'*ones(1,n); A = abs(A - A'); A = rho .^ A; if imag(rho) A = conj(tril(A,-1)) + triu(A); end End of kms.m echo krylov.m 1>&2 cat >krylov.m <<'End of krylov.m' function B = krylov(A, x, j) %KRYLOV KRYLOV(A, x, j) is the Krylov matrix % [x, Ax, A^2x, ..., A^(j-1)x], where A is an n-by-n matrix and % x is an n-vector. % Defaults: x = ONES(n,1), j = n. % KRYLOV(n) is the same as KRYLOV(RAND(n)). % Reference: G.H. Golub and C.F. Van Loan, Matrix Computations, Johns Hopkins % University Press, Baltimore, Maryland, 1983, p. 224. [n, n] = size(A); if n == 1 % Handle special case A = scalar. n = A; A = rand(n); end if nargin < 3, j = n; end if nargin < 2, x = ones(n,1); end B = ones(n,j); B(:,1) = x(:); for i=2:j B(:,i) = A*B(:,i-1); end End of krylov.m echo lauchli.m 1>&2 cat >lauchli.m <<'End of lauchli.m' function A = lauchli(n, mu) %LAUCHLI LAUCHLI(N, MU) is the (N+1)-by-N matrix [ONES(1,N); MU*EYE(N))]. % It is a well-known example in least squares and other problems % that indicates the dangers of forming A'*A. % MU defaults to SQRT(EPS). % Reference: P. Lauchli, Jordan-Elimination und Ausgleichung nach % kleinsten Quadraten, Numer. Math, 3 (1961), pp. 226-240. if nargin < 2, mu = sqrt(eps); end A = [ones(1,n); mu*eye(n)]; End of lauchli.m echo lehmer.m 1>&2 cat >lehmer.m <<'End of lehmer.m' function A = lehmer(n) %LEHMER A = LEHMER(N) is the symmetric positive definite N-by-N matrix with % A(i,j) = i/j for j>=i. % A is totally nonnegative. INV(A) is tridiagonal, and explicit % formulas are known for its entries. % N <= COND(A) <= 4*N*N. % References: % M. Newman and J. Todd, The evaluation of matrix inversion % programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476. % Solutions to problem E710 (proposed by D.H. Lehmer): The inverse % of a matrix, Amer. Math. Monthly, 53 (1946), pp. 534-535. A = ones(n,1)*(1:n); A = A./A'; A = tril(A) + tril(A,-1)'; End of lehmer.m echo lotkin.m 1>&2 cat >lotkin.m <<'End of lotkin.m' function A = lotkin(n) %LOTKIN A = LOTKIN(N) is the Hilbert matrix with its first row altered to % all ones. A is unsymmetric, ill-conditioned, and has many negative % eigenvalues of small magnitude. The inverse has integer entries % and is known explicitly. % Reference: M. Lotkin, A set of test matrices, MTAC, 9 (1955), % pp. 153-161. A = hilb(n); A(1,:) = ones(1,n); End of lotkin.m echo matrix.m 1>&2 cat >matrix.m <<'End of matrix.m' function A = matrix(k, n) %MATRIX MATRIX(K, N) is the N-by-N instance of the matrix number K % in the collection [1] (including some of the matrices provided % with MATLAB), with all other parameters set to their default. % N.B. Only those matrices which take an arbitrary dimension N % are included (thus GALLERY is omitted, for example). % MATRIX(K) is a string containing the name of the K'th matrix. % MATRIX(0) is the number of matrices, i.e. the upper limit for K. % Thus to set A to each N-by-N test matrix in turn use a loop like % for k=1:matrix(0) % A = matrix(k, N); % Aname = matrix(k); % The name of the matrix % end % MATRIX(-1) returns the version number and date of the collection. % [1] N.J. Higham, A collection of test matrices in MATLAB, Technical Report, % Department of Computer Science, Cornell University, Ithaca, NY 14853, % 1989; also Numerical Analysis Report, Department of Mathematics, % University of Manchester, Manchester, M13 9PL, England, 1989. % Matrices omitted are: gallery, hadamard, hanowa, lauchli, wathen, wilk. % Matrices provided with MATLAB that are included here: invhilb, magic. % Set up string array a few lines at a time to avoid `input buffer line % overflow'. matrices = [ 'cauchy '; 'chebspec'; 'chebvand'; 'chow '; 'circul '; ... 'clement '; 'compan '; 'condex '; 'cycol '; 'dingdong'; ... 'dorr '; 'dramadah'; 'fiedler '; 'forsythe'; 'frank '; ... 'gear '; 'gfpp '; 'hilb '; 'invhilb '; 'invol ';]; matrices = [matrices 'ipjfact '; 'jordan '; 'kahan '; 'kms '; 'krylov '; ... 'lehmer '; 'lotkin '; 'magic '; 'minij '; 'moler '; ... 'ohess '; 'orthog '; 'pascal '; 'pei '; 'rando '; ... 'randsvd '; 'riemann '; 'tridiag '; 'triw '; 'vand ';]; if nargin == 1 if k == 0 [A, temp] = size(matrices); elseif k > 0 A = matrices(k,:); else A = 'Version 1.0, July 4 1989'; % Version number and date of collection. end else A = eval( [matrices(k,:) '(n)'] ); end End of matrix.m echo minij.m 1>&2 cat >minij.m <<'End of minij.m' function A = minij(n) %MINIJ A = MINIJ(N) is the N-by-N symmetric positive definite matrix with % A(i,j) = min(i,j). % Properties, variations: % INV(A) is tridiagonal: it is minus the second difference matrix % except its (N,N) element is 1. % 2*A-ONES(A) (Givens' matrix) has tridiagonal inverse and % eigenvalues .5*sec^2([2r-1)PI/4N], r=1:N. % (N+1)ONES(A)-A has elements max(i,j), and also has a tridiagonal % inverse. A = min( ones(n,1)*(1:n), (1:n)'*ones(1,n) ); End of minij.m echo moler.m 1>&2 cat >moler.m <<'End of moler.m' function A = moler(n, alpha) %MOLER A = MOLER(N, ALPHA) is the symmetric positive definite N-by-N matrix % U'*U where U = TRIW(N,-1, ALPHA). % For ALPHA = -1 (the default) A(i,j) = min(i,j)-2, A(i,i) = i. % A has one small eigenvalue. % Nash (1979) attributes the ALPHA = -1 matrix to Moler. % Reference: J.C. Nash, Compact Numerical Methods for Computers: Linear % Algebra and Function Minimisation, Adam Hilger, Bristol, and John Wiley, % New York, 1979. pp. 76, 210. if nargin == 1, alpha = -1; end A = triw(n, alpha)'*triw(n, alpha); End of moler.m echo ohess.m 1>&2 cat >ohess.m <<'End of ohess.m' function H = ohess(x) %OHESS H = OHESS(N) is an N-by-N real, random, orthogonal upper Hessenberg % matrix. % Alternatively, H = OHESS(X), where X is an arbitrary real N-vector % (N>1) constructs H non-randomly using the elements of X as parameters. % In both cases H is constructed via a product of N-1 Givens rotations. % Note: See Gragg (1986) for how to represent an N-by-N (complex) % unitary Hessenberg matrix with positive subdiagonal elements in terms % of 2N-1 real parameters (the Schur parametrization). % This M-file handles the real case only and is intended simply as a % convenient way to generate random or non-random orthogonal Hessenberg % matrices. % Reference: W.B. Gragg, The QR algorithm for unitary Hessenberg matrices, % J. Comp. Appl. Math., 16 (1986), pp. 1-8. if any(imag(x)), error('Parameter must be real.'), end n = max(size(x)); if n == 1 % Handle scalar x. n = x; rand('uniform') x = rand(n-1,1)*2*pi; H = eye(n); rand('normal'); H(n,n) = sign(rand); else H = eye(n); H(n,n) = sign(x(n)) + (x(n)==0); % Second term ensures H(n,n) nonzero. end for i=n:-1:2 % Apply Givens rotation through angle x(i-1). theta = x(i-1); c = cos(theta); s = sin(theta); H( [i-1 i], : ) = [ c*H(i-1,:)+s*H(i,:) -s*H(i-1,:)+c*H(i,:) ]; end End of ohess.m echo orthog.m 1>&2 cat >orthog.m <<'End of orthog.m' function q = orthog(n, k) %ORTHOG Orthogonal and nearly orthogonal matrices. % Q = ORTHOG(N, K) selects the K'th type of matrix of order N. % K > 0 for exactly orthogonal matrices, K < 0 for diagonal scalings of % orthogonal matrices. % Available types: (K = 1 is the default) % K = 1: q(i,j) = SQRT(2/(n+1)) * SIN( i*j*PI/(n+1) ) % Symmetric eigenvector matrix for second difference matrix. % K = 2: q(i,j) = 2/SQRT(2*n+1)) * SIN( 2*i*j*PI/(2*n+1) ) % Symmetric. % K = 3: q(r,s) = EXP(2*PI*i*(r-1)*(s-1)/n) / SQRT(n) (i=SQRT(-1)) % Unitary, the Fourier matrix. % K = -1: q(i,j) = COS( (i-1)*(j-1)*PI/(n-1) ) % Chebyshev Vandermonde-like matrix, based on extrema of T(n-1). % K = -2: q(i,j) = COS( (i-1)*(j-1/2)*PI/n) ) % Chebyshev Vandermonde-like matrix, based on zeros of T(n). % 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 nargin == 1, k = 1; end if k == 1 % E'vectors second difference matrix m = (1:n)'*(1:n) * (pi/(n+1)); q = sin(m) * sqrt(2/(n+1)); elseif k == 2 m = (1:n)'*(1:n) * (2*pi/(2*n+1)); q = sin(m) * (2/sqrt(2*n+1)); elseif k == 3 % Vandermonde based on roots of unity m = 0:n-1; q = exp(m'*m*2*pi*sqrt(-1)/n) / sqrt(n); elseif k == -1 % extrema of T(n-1) m = (0:n-1)'*(0:n-1) * (pi/(n-1)); q = cos(m); elseif k == -2 % zeros of T(n) m = (0:n-1)'*(.5:n-.5) * (pi/n); q = cos(m); end End of orthog.m echo pascal.m 1>&2 cat >pascal.m <<'End of pascal.m' function P = pascal(n, k) %PASCAL PASCAL(N) is the Pascal matrix of order N: a symmetric positive % definite matrix with integer entries, made up from Pascal's % triangle. Its inverse has integer entries. % [Conjecture by NJH: the Pascal matrix is totally positive.] % PASCAL(N,1) is the lower triangular Cholesky factor (up to signs % of columns) of the Pascal matrix. It is involutary (is its own % inverse). % PASCAL(N,2) is a transposed and permuted version of PASCAL(N,1) % which is a cube root of the identity. if nargin == 1, k = 0; end P = diag( (-1).^[0:n-1] ); P(:, 1) = ones(n,1); % Generate the Pascal Cholesky factor (up to signs). for j=2:n-1 for i=j+1:n P(i,j) = P(i-1,j) - P(i-1,j-1); end end if k == 0 P = P*P'; elseif k == 2 P = P'; P = P(n:-1:1,:); for i=1:n-1 P(i,:) = -P(i,:); P(:,i) = -P(:,i); end if n/2 == round(n/2), P = -P; end end End of pascal.m echo pei.m 1>&2 cat >pei.m <<'End of pei.m' function P = pei(n, alpha) %PEI PEI(N, ALPHA), where ALPHA is a scalar, is the Pei matrix with % parameter ALPHA, i.e. the symmetric matrix ALPHA*EYE(N) + ONES(N). % If ALPHA is omitted then ALPHA = 1 is used. % The matrix is symmetric, and singular for ALPHA = 0, -N. % Reference: M.L. Pei, A test matrix for inversion procedures, % Comm. ACM, 5 (1962), p. 508. if nargin == 1, alpha = 1; end P = alpha*eye(n) + ones(n); End of pei.m echo qmult.m 1>&2 cat >qmult.m <<'End of qmult.m' function B = qmult(A) %QMULT QMULT(A) is Q*A where Q is a random real orthogonal matrix from the % Haar distribution, of dimension the number of rows in A. % Special case: if A is a scalar then QMULT(A) is the same as % QMULT(EYE(A)). N.B. This routine sets RAND('NORMAL'). % Reference: G.W. Stewart, The efficient generation of random % orthogonal matrices with an application to condition estimators, % SIAM J. Numer. Anal., 17 (1980), 403-409. [n, m] = size(A); % Handle scalar A. if max(n,m) == 1 n = A; A = eye(n); end d = zeros(n); rand('normal') for k = n-1:-1:1 % Generate random Householder transformation. x = rand(n-k+1,1); s = norm(x); sgn = sign(x(1)) + (x(1)==0); % Modification for sign(1)=1. s = sgn*s; d(k) = -sgn; x(1) = x(1) + s; beta = s*x(1); % Apply the transformation to A. y = x'*A(k:n,:); A(k:n,:) = A(k:n,:) - x*(y/beta); end % Tidy up signs. for i=1:n-1 A(i,:) = d(i)*A(i,:); end A(n,:) = A(n,:)*sign(rand); B = A; End of qmult.m echo rando.m 1>&2 cat >rando.m <<'End of rando.m' function A = rando(n, k) %RANDO A = RANDO(N, K) is a random N-by-N matrix with elements from one of % the following discrete distributions (default K=1): % K=1: A(i,j) = 0 or 1 with equal probability, % K=2: A(i,j) = -1 or 1 with equal probability, % K=3: A(i,j) = -1, 0 or 1 with equal probability. % N may be a 2-vector, in which case the matrix is N(1)-by-N(2). if nargin < 2, k = 1; end m = n(1); % Parameter n specifies dimension: m-by-n. n = n(max(size(n))); rand('uniform') if k == 1 % {0, 1} A = floor( rand(m,n) + .5 ); elseif k == 2 % {-1, 1} A = 2*floor( rand(m,n) + .5 ) - 1; elseif k == 3 % {-1, 0, 1} A = round( 3*rand(m,n) - 1.5 ); end End of rando.m echo randsvd.m 1>&2 cat >randsvd.m <<'End of randsvd.m' function A = randsvd(n, kappa, mode, kl, ku) %RANDSVD Random matrices with pre-assigned singular values. % RANDSVD(N, KAPPA, MODE, KL, KU) is a (banded) random matrix of order N % with COND(A) = KAPPA and singular values from the distribution MODE. % N may be a 2-vector, in which case the matrix is N(1)-by-N(2). % Available types: % MODE = 1: one large singular value, % MODE = 2: one small singular value, % MODE = 3: geometrically distributed singular values, % MODE = 4: arithmetically distributed singular values, % MODE = 5: random singular values with unif. dist. logarithm. % If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS). % If MODE < 0 then the effect is as for ABS(MODE) except that in the % original matrix of singular values the order of the diagonal entries is % reversed: small to large instead of large to small. % KL and KU are the lower and upper bandwidths respectively; if they are % omitted a full matrix is produced. If only KL is present, KU defaults % to KL. % This routine is similar to the more comprehensive Fortran routine xLATMS % in the following reference: % J.W. Demmel and A. McKenney, A test matrix generation suite, % LAPACK Working Note #9, Courant Institute of Mathematical Sciences, % New York, 1989. if nargin < 2, kappa = sqrt(1/eps); end if nargin < 3, mode = 3; end if nargin < 4, kl = n-1; end % Full matrix. if nargin < 5, ku = kl; end % Same upper and lower bandwidths. if kappa < 1 error('Condition number must be at least 1!') end p = min(n); m = n(1); % Parameter n specifies dimension: m-by-n. n = n(max(size(n))); if p == 1 % Handle case where A is a vector. rand('normal') A = rand(m, n); A = A/norm(A); return end j = abs(mode); % Set up vector sigma of singular values. if j == 3 factor = kappa^(-1/(p-1)); sigma = factor.^[0:p-1]; elseif j == 4 sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa); elseif j == 5 rand('uniform') sigma = exp( -rand(p,1)*log(kappa) ); elseif j == 2 sigma = ones(p,1); sigma(p) = 1/kappa; elseif j == 1 sigma = ones(p,1)./kappa; sigma(1) = 1; end % Convert to diagonal matrix of singular values. if mode < 0 sigma = sigma(p:-1:1) end sigma = diag(sigma); if m ~= n sigma(m, n) = 0; % Expand to m-by-n diagonal matrix. end if kl == 0 & ku == 0 % Diagonal matrix requested - nothing more to do. A = sigma; return end % A = U*sigma*V, where U, V are random orthogonal matrices from the % Haar distribution. A = qmult(sigma'); A = qmult(A'); if kl < n-1 | ku < n-1 % Bandwidth reduction. A = bandred(A, kl, ku); end End of randsvd.m echo riemann.m 1>&2 cat >riemann.m <<'End of riemann.m' function A = riemann(n) %RIEMANN A = RIEMANN(N) is a N-by-N matrix associated with the % Riemann hypothesis, which is true if and only if % DET(A) = O( N! N^(-1/2+epsilon) ) for every epsilon>0 % ('!' denotes factorial). % A = B(2:N+1, 2:N+1), where B(i,j) = i-1 if i divides j and % -1 otherwise. % Properties include, with M = N+1: % Each eigenvalue E(i) satisfies ABS(E(i)) <= M - 1/M. % i <= E(i) <= i+1 with at most M-SQRT(M) exceptions. % All integers in the interval (M/3, M/2] are eigenvalues. % Reference: F. Roesler, Riemann's hypothesis as an % eigenvalue problem, Linear Algebra and Appl., 81 (1986), % pp. 153-198. n = n+1; i = (2:n)'*ones(1,n-1); j = i'; A = i .* (~rem(j,i)) - ones(n-1); End of riemann.m echo rq.m 1>&2 cat >rq.m <<'End of rq.m' function z = rq(a,x) %RQ RQ(A,x) is the Rayleigh quotient of A and x, x'*A*x/(x'*x). z = x'*a*x/(x'*x); End of rq.m echo see.m 1>&2 cat >see.m <<'End of see.m' function see(a, k) %SEE Pictures of a matrix and its inverse. % SEE(A) displays MESH(A), MESH(INV(A)), PLOT(A), and FV((A)) % in four subplot windows. % SEE(A,1) uses the full screen for each picture, with a PAUSE in % between each one. if nargin < 2, k = 0; end clg if k == 1, mesh(a), pause b = inv(a); mesh(b), pause plot(a), pause fv(a); else subplot(221) mesh(a) b = inv(a); mesh(b) plot(a) fv(a); % subplot end End of see.m echo seqa.m 1>&2 cat >seqa.m <<'End of seqa.m' function y = seqa(a, b, n) %SEQA Generate an additive sequence. % Y = SEQA(A, B, N) produces a row vector comprising N equally % spaced numbers starting at A and finishing at B. % If N is omitted then 10 points are generated. if nargin == 2, n = 10; end if n <= 1 y = a; return end y = [a+(0:n-2)*(b-a)/(n-1), b]; End of seqa.m echo seqcheb.m 1>&2 cat >seqcheb.m <<'End of seqcheb.m' function x = seqcheb(n, k) %SEQCHEB Sequence of points related to Chebyshev polynomials, T_N. % X = SEQCHEB(N, K) produces a row vector of length N. % K = 1: zeros of T_N, (the default) % K = 2: extrema of T_{N-1}. if nargin == 1, k = 1; end if k == 1 % Zeros of T_n i = 1:n; j = .5*ones(i); x = cos( (i-j) * (pi/n) ); elseif k == 2 % Extrema of T_(n-1) i = 0:n-1; x = cos( i * (pi/(n-1)) ); end End of seqcheb.m echo seqm.m 1>&2 cat >seqm.m <<'End of seqm.m' function y = seqm(a, b, n) %SEQM Generate a multiplicative sequence. % Y = SEQM(A, B, N) produces a row vector comprising N % logarithmically equally spaced numbers, starting at A~=0 % and finishing at B~=0. % If A*B < 0 and N > 2 then complex results are produced. % If N is omitted then 10 points are generated. if nargin == 2, n = 10; end if n <= 1 y = a; return end p = [0:n-2]/(n-1); r = (b/a).^p; y = [a*r, b]; End of seqm.m echo show.m 1>&2 cat >show.m <<'End of show.m' function show(x) %SHOW SHOW(X) displays X in 'FORMAT +' form. format + disp(x) format End of show.m echo skew.m 1>&2 cat >skew.m <<'End of skew.m' function S = skew(A) %SKEW SKEW(A) is the skew-symmetric (Hermitian) part of A, (A-A')/2. % It is the nearest skew-symmetric (Hermitian) matrix to A in both the % 2- and the Frobenius norms. S = (A - A')./2; End of skew.m echo sparse.m 1>&2 cat >sparse.m <<'End of sparse.m' function sparse(A, s, x, y) % SPARSE SPARSE(A, s) plots the sparsity pattern of a matrix A, showing % s = 'x', 'o', '+' or '.' for every nonzero element. % SPARSE(A,s,x,y) overlays the current graphics screen taking (x,y) % as the origin. if nargin <= 2 x = 1; y = 1; end if nargin == 1, s = 'x'; end [m, n] = size(A); A = A(m:-1:1,:)'; [m, n] = size(A); if nargin <= 2 clg hold off axis( [1, m, 1, n]) else hold on end for j=n:-1:1 plot( (j+y-1) .* [zeros(x-1,1); (A(:,j)~=0)], ['w' s] ); hold on end hold off End of sparse.m echo sparsify.m 1>&2 cat >sparsify.m <<'End of sparsify.m' function A = sparsify(A, p) %SPARSIFY S = SPARSIFY(A, P) is A with elements randomly set to zero % (S = S' if A = A', i.e. symmetry is preserved). % Each element has probability P of being zeroed. % Thus on average 100*P percent of the elements of A will be zeroed. % Default: P = 0.25. if nargin < 2, p = 0.25; end if p<0 | p>1, error('Second parameter must be between 0 and 1 inclusive'), end rand('uniform') if norm(A-A',1) A = A .* (rand(A) > p); % Unsymmetric case else A = triu(A) .* (rand(A) > p); % Preserve symmetry A = (A + A')/2; end End of sparsify.m echo sub.m 1>&2 cat >sub.m <<'End of sub.m' function S = sub(A, i, j) %SUB Principal submatrix of A. % SUB(A,i,j) is A(i:j,i:j). % SUB(A,i) is the leading principal submatrix of order i, A(1:i,1:i), % if i>0, and the trailing principal submatrix of order ABS(i) if i<0. if nargin == 2 if i >= 0 S = A(1:i, 1:i); else n = min(size(A)); S = A(n+i+1:n, n+i+1:n); end else S = A(i:j, i:j); end End of sub.m echo symm.m 1>&2 cat >symm.m <<'End of symm.m' function S = symm(A) %SYMM SYMM(A) is the symmetric (Hermitian) part of A, (A+A')/2. % It is the nearest symmetric (Hermitian) matrix to A in both the % 2- and the Frobenius norms. S = (A + A')./2; End of symm.m echo tridiag.m 1>&2 cat >tridiag.m <<'End of tridiag.m' function T = tridiag(n, x, y, z) %TRIDIAG TRIDIAG(X,Y,Z) is the tridiagonal matrix with subdiagonal X, % diagonal Y, and superdiagonal Z. % X and Z must be vectors of dimension one less than Y. % Alternatively TRIDAG(N,C,D,E), where C, D, and E are all % scalars, yields the Toeplitz tridiagonal matrix of order N % with subdiagonal elements C, diagonal elements D, and superdiagonal % elements E. This matrix has eigenvalues (Todd 1978) % D + 2*SQRT(C*E)*COS(k*PI/(N+1)), k=1:N. % TRIDIAG(N) is the same as TRIDIAG(N,-1,2,-1), which is % a symmetric positive definite M-matrix (the negative of the % second difference matrix). % Reference: J. Todd, Basic Numerical Mathematics, Vol.2: Numerical Algebra, % Academic Press, New York, 1978, p. 155. if nargin == 1, x = -1; y = 2; z = -1; end if nargin == 3, z = y; y = x; x = n; end x = x(:); y = y(:); z = z(:); % Force column vectors. if max( [ size(x) size(y) size(z) ]' ) == 1 x = x*ones(n-1,1); z = z*ones(n-1,1); y = y*ones(n,1); else [nx, m] = size(x); [ny, m] = size(y); [nz, m] = size(z); if (ny - nx - 1) | (ny - nz -1) error('Dimensions of vector arguments are incorrect.') end end T = diag(x, -1) + diag(y) + diag(z, 1); End of tridiag.m echo triw.m 1>&2 cat >triw.m <<'End of triw.m' function t = triw(n, alpha, k) %TRIW TRIW(N, ALPHA, K) is the upper triangular matrix with ones on % the diagonal and ALPHAs on the first K >= 0 superdiagonals. % N may be a 2-vector, in which case the matrix is N(1)-by-N(2) and % upper trapezoidal. % Defaults: ALPHA = -1, % K = N-1 (full upper triangle). % TRIW(N) is a matrix discussed by Wilkinson. m = n(1); % Parameter n specifies dimension: m-by-n. n = n(max(size(n))); if nargin < 3, k = n-1; end if nargin < 2, alpha = -1; end if max(size(alpha)) ~= 1 error('Second argument must be a scalar.') end t = tril( eye(m,n) + alpha*triu(ones(m,n), 1), k); End of triw.m echo vand.m 1>&2 cat >vand.m <<'End of vand.m' function V = vand(m, p) %VAND V = VAND(P), where P is a vector, produces the (primal) % Vandermonde matrix based on the points P, i.e. V(i,j) = P(j)^(i-1). % VAND(M,P) is a rectangular version of VAND(P) with M rows. % Special case: If P is a scalar then P equally spaced points on [0,1] % are used. % Reference: N.J. Higham, Stability analysis of algorithms for solving % confluent Vandermonde-like systems, to appear in SIAM J. Matrix Anal. Appl. if nargin == 1, p = m; end n = max(size(p)); % Handle scalar p. if n == 1 n = p; p = seqa(0,1,n); end if nargin == 1, m = n; end p = p(:).'; % Ensure p is a row vector. V = ones(m,n); for i=2:m V(i,:) = p.*V(i-1,:); end End of vand.m echo wathen.m 1>&2 cat >wathen.m <<'End of wathen.m' function A = wathen(nx, ny, k) %WATHEN A = WATHEN(NX,NY) is a random N-by-N finite element matrix % where N = 3*NX*NY + 2*NX + 2*NY + 1. % A is precisely the "consistent mass matrix" for a regular NX-by-NY % grid of 8-node (serendipity) elements in 2 space dimensions. % A is symmetric positive definite for any (positive) values of % the "density", RHO(NX,NY), which is chosen randomly in this routine. % In particular, if D=DIAG(DIAG(A)), then % 0.25 <= EIG(INV(D)*A) <= 4.5 % for any positive integers NX and NY and any densities RHO(NX,NY). % This diagonally scaled matrix is returned by WATHEN(NX,NY,1). % Reference: A.J.Wathen, Realistic eigenvalue bounds for the Galerkin % mass matrix, IMA J. Numer. Anal., 7 (1987), pp. 449-457. % BEWARE - this is a sparse matrix and it quickly gets large! if nargin < 2, error('Two dimensioning arguments must be specified.'), end if nargin < 3, k = 0; end e1 = [6,-6,2,-8;-6,32,-6,20;2,-6,6,-6;-8,20,-6,32]; e2 = [3,-8,2,-6;-8,16,-8,20;2,-8,3,-8;-6,20,-8,16]; e = [e1,e2;e2',e1]/45; n = 3*nx*ny+2*nx+2*ny+1; A = zeros(n); rand('uniform') RHO = 100*rand(nx,ny); for j=1:ny for i=1:nx nn(1) = 3*j*nx+2*i+2*j+1; nn(2) = nn(1)-1; nn(3) = nn(2)-1; nn(4) = (3*j-1)*nx+2*j+i-1; nn(5) = 3*(j-1)*nx+2*i+2*j-3; nn(6) = nn(5)+1; nn(7) = nn(6)+1; nn(8) = nn(4)+1; em = e*RHO(i,j); for krow=1:8 for kcol=1:8 A(nn(krow),nn(kcol)) = A(nn(krow),nn(kcol))+em(krow,kcol); end end end end if k == 1 A = diag(diag(A)) \ A; end End of wathen.m echo wilk.m 1>&2 cat >wilk.m <<'End of wilk.m' function [A, b] = wilk(n) %WILK Various specific matrices devised/discussed by Wilkinson. % [A, b] = WILK(N) is the matrix or system of order N. % N=3: upper triangular system Ux=b illustrating inaccurate solution. % N=4: lower triangular system Lx=b, ill-conditioned. % N=5: HILB(6)(1:5,2:6)*1.8144. Symmetric positive definite. % N=21: W21+, tridiagonal. Eigenvalue problem. % References: % J.H. Wilkinson, Error analysis of direct methods of matrix inversion, % J. Assoc. Comput. Mach., 8 (1961), pp. 281-330. % J.H. Wilkinson, Rounding Errors in Algebraic Processes, Notes on Applied % Science No. 32, Her Majesty's Stationery Office, London, 1963. % J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University % Press, 1965. if n == 3 % Wilkinson (1961) p.323. A = [ 1e-10 .9 -.4 0 .9 -.4 0 0 1e-10]; b = [ 0 0 1]'; elseif n == 4 % Wilkinson (1963) p.105. A = [0.9143e-4 0 0 0 0.8762 0.7156e-4 0 0 0.7943 0.8143 0.9504e-4 0 0.8017 0.6123 0.7165 0.7123e-4]; b = [0.6524 0.3127 0.4186 0.7853]'; elseif n == 5 % Wilkinson (1965), p.234. A = hilb(6); A = A(1:5, 2:6)*1.8144; elseif n == 21 % Taken from gallery.m. Wilkinson (1965), p.308. E = diag(ones(n-1,1),1); m = (n-1)/2; A = diag(abs(-m:m)) + E + E'; else error('Sorry, that value of N is not available.') end End of wilk.m