# To unbundle, sh this file
echo augment.m 1>&2
cat >augment.m <<'End of augment.m'
function C = augment(A);
%AUGMENT  AUGMENT(A) is the square matrix [EYE(m) A; A' ZEROS(n)] of dimension
%         m+n, where A is m-by-n.  It is the symmetric and indefinite
%         coefficient matrix of the augmented system associated with a least
%         squares problem minimize NORM(A*x-b).
%         Special case: if A is a scalar, n say, then AUGMENT(A) is the same
%                       as AUGMENT(RAND(p,q)) where n = p+q and p = ROUND(n/2),
%                       that is, a random augmented matrix of dimension n is
%                       produced.
%         The eigenvalues of AUGMENT(A) are given in terms of the singular
%         values s(i) of A (where m>n) by
%                  1/2 +/- SQRT( s(i)^2 + 1/4 ),  i=1:n  (2n eigenvalues),
%                  1,  (m-n eigenvalues).
%         If m<n then the first expression provides 2m eigenvalues and the
%         remaining n-m eigenvalues are zero.

%   Reference:  G.H. Golub and C.F. Van Loan, Matrix Computations, Second
%   Edition, Johns Hopkins University Press, Baltimore, Maryland, 1989,
%   sec. 5.6.4.

[m, n] = size(A);

if max(m,n) == 1
   n = A;
   p = round(n/2);
   q = n - p;
   A = rand(p,q);
   m = p; n = q;
end

C = [eye(m) A; A' zeros(n)];
End of augment.m
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 routine is called by RANDSVD.
%         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).

%        References:
%        D.E. Knuth, The Art of Computer Programming, Volume 1,
%        Fundamental Algorithms, Second Edition, Addison-Wesley, Reading,
%        Massachusetts, 1973, p. 36.
%        E.E. Tyrtyshnikov, Cauchy-Toeplitz matrices and some applications,
%        Linear Algebra and Appl., 149 (1991), pp. 1-18.

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, SIAM J. Matrix Anal. Appl., 11 (1990),
%   pp. 23-41.

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.
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    CHOP(X, t) is the matrix obtained by rounding the elements of X
%        to t significant binary places.
%        Default is t = 24, corresponding to IEEE single precision.

if nargin < 2, t = 24; 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 unsymmetric, while for K=1 it
%          is symmetric.  CLEMENT(N, 1) is diagonally similar to CLEMENT(N).

%          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.

%          References:
%          P.A. Clement, A class of triple-diagonal matrices for test
%          purposes, SIAM Review, 1 (1959), pp. 50-52.
%          O. Taussky and J. Todd, Another look at a matrix of Mark Kac, 
%          Linear Algebra and Appl., 150 (1991), pp. 341-360.

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).

%  References:
%  J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press,
%       1965, p. 12.
%  C. Kenney and A.J. Laub, Controllability and stability radii for companion
%       form systems, Math. Control Signals Systems, 1 (1988), pp. 239-256.
%       (Gives explicit formulas for the singular values of COMPAN(P).)

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-2/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 Hankel 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 row 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.
%      The amount of diagonal dominance is given by (ignoring rounding errors):
%            COMP(C)*ONES(N,1) = THETA*(N+1)^2 * [1 0 0 ... 0 1]'.

%      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:
%           FIEDLER(N) has a dominant positive eigenvalue and all the other
%                      eigenvalues are negative (Szego 1936).
%           Explicit formulas for INV(A) and DET(A) are given in (Todd 1977)
%           and attributed to Fiedler.  These indicate that INV(A) is
%           tridiagonal except for nonzero (1,n) and (n,1) elements.
%           [I think these formulas are valid only if the elements of
%           C are in increasing or decreasing order---NJH.]

%    References: G. Szego, Solution to problem 3705, Amer. Math. Monthly,
%    43 (1936), pp. 246-259.
%    J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
%    Birkhauser, Basel, and Academic Press, New York, 1977, p. 159.

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 polynomial
%                  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 = 16.
%       (For a `publication quality' picture, set THMAX higher, say 32.)
%       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.   Note that NORM(F,INF) approximates the numerical radius,
%       max {abs(z): z is in the field of values of A}.

%       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).

%       References:
%       R.A. Horn and C.R. Johnson, Topics in Matrix Analysis, Cambridge
%            University Press, 1991, section 1.5.
%       A.S. Householder, The Theory of Matrices in Numerical Analysis,
%            Blaisdell, New York, 1964, section 3.3.
%       C.R. Johnson, Numerical determination of the field of values of a
%            general complex matrix, SIAM J. Numer. Anal., 15 (1978),
%            pp. 595-602.

%       Original version by A. Ruhe.  Modified by N.J. Higham.

if nargin < 2, nk = 1; end
if nargin < 3, thmax = 16; end

iu = sqrt(-1);
[n, p] = size(B);
f = [];
z = zeros(2*thmax+1,1);
e = eig(B);

% Filter out cases where B is Hermitian or skew-Hermitian, for efficiency.
if norm(skew(B),1) == 0

   f = [min(e) max(e)];             

elseif norm(symm(B),1) == 0

   e = imag(e);
   f = [min(e) max(e)];
   e = iu*e; f = iu*f;
   
else

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

end

if nargin < 4

   % Set x and y axis ranges so both have the same length.
   xmin = min(real(f)); xmax = max(real(f));
   ymin = min(imag(f)); ymax = max(imag(f));
   if xmax-xmin >= ymax-ymin
      ymid = (ymin + ymax)/2;
      ymin =  ymid - (xmax-xmin)/2; ymax = ymid + (xmax-xmin)/2;
   else
      xmid = (xmin + xmax)/2;
      xmin = xmid - (ymax-ymin)/2; xmax = xmid + (ymax-ymin)/2;
   end
   axis('square')
   % Scale ranges by 1+2*alpha to give extra space around edges of plot.
   alpha = 0.1;
   axis( [xmin-alpha*(xmax-xmin) xmax+alpha*(xmax-xmin) ...
          ymin-alpha*(ymax-ymin) ymax+alpha*(ymax-ymin) ] );

   plot(real(f), imag(f),'w')      % Plot the field of values
   hold on
   plot(real(e), imag(e), 'xg')    % Plot the eigenvalues too.
   axis([1 2 3 4]); axis; axis('normal');  % Normal ratio and auto-scaling.
   hold off

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 gersh.m 1>&2
cat >gersh.m <<'End of gersh.m'
function  [G, e] = gersh(A, nolplot)
%GERSH    GERSH(A) draws the Gershgorin disks for the matrix A.
%         The eigenvalues are plotted as 'x'.
%         Alternative usage: [G, E] = GERSH(A, 1) suppresses the plot
%         and returns the data in G, with A's eigenvalues in E.
%         GERSH(LESP(N)) is interesting!

n = max(size(A));
m = 40;
G = zeros(m,n);

d = diag(A);
r = sum( abs( A.'-diag(d) ) )';
e = eig(A);

radvec = exp(i * seqa(0,2*pi,m)');

for j=1:n
    G(:,j) = d(j)*ones(radvec) + r(j)*radvec;
end

if nargin < 2

   % Set x and y axis ranges so both have the same length.
   xmin = min(real(G(:))); xmax = max(real(G(:)));
   ymin = min(imag(G(:))); ymax = max(imag(G(:)));
   if xmax-xmin >= ymax-ymin
      ymid = (ymin + ymax)/2;
      ymin =  ymid - (xmax-xmin)/2; ymax = ymid + (xmax-xmin)/2;
   else
      xmid = (xmin + xmax)/2;
      xmin = xmid - (ymax-ymin)/2; xmax = xmid + (ymax-ymin)/2;
   end
   axis('square')
   % Scale ranges by 1+2*alpha to give extra space around edges of plot.
   alpha = 0.1;
   axis( [xmin-alpha*(xmax-xmin) xmax+alpha*(xmax-xmin) ...
          ymin-alpha*(ymax-ymin) ymax+alpha*(ymax-ymin) ] );

   for j=1:n
       plot(real(G(:,j)), imag(G(:,j)),'-c5')      % Plot the disks.
       hold on
   end
   plot(real(e), imag(e), 'xg')    % Plot the eigenvalues too.
   axis([1 2 3 4]); axis; axis('normal');  % Normal ratio and auto-scaling.
   hold off

end
End of gersh.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, that is, a matrix H
%          with elements 1 or -1 such that H*H' = N*EYE(N).
%          An N-by-N Hadamard matrix with N>2 exists only if 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
   msg = 'This function handles only the cases where N, N/12 or N/20 ';
   msg = [msg 'is a power of 2.'];
   error(msg)
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=2M must be even)
%                      [d*EYE(M)   -DIAG(1:M)
%                       DIAG(1:M)   d*EYE(M)]
%        It has complex eigenvalues lambda(k) = d +/- k*i  (1 <= k <= M).
%        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.
%       COND(HILB(N)) grows like EXP(3.5*N).
%       See INVHILB (standard MATLAB routine) for the exact inverse, which
%       has integer entries.
%       HILB(N) is symmetric positive definite, totally positive, and a
%       Hankel matrix.

%       This routine is shorter and faster than the one supplied with MATLAB.

%       References:
%       M.-D. Choi, Tricks or treats with the Hilbert matrix, Amer. Math.
%           Monthly, 90 (1983), pp. 301-312.
%       M. Newman and J. Todd, The evaluation of matrix inversion
%           programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476.
%       D.E. Knuth, The Art of Computer Programming,
%           Volume 1, Fundamental Algorithms, Second Edition, Addison-Wesley,
%           Reading, Massachusetts, 1973, p. 37.

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: If x = 0 then v = 0, beta = 1 is returned.
%            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)|).
                       
[m, n] = size(x);
if n > 1, error('Argument must be a column vector.'), end

s = norm(x) * (sign(x(1)) + (x(1)==0));    % Modification for sign(0)=1.
v = x;
if s == 0, beta = 1; return, end           % Quit if x is the zero vector.
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 1.2.
%       The useful range of THETA is 0 < THETA < PI.

%       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 = 1.2; 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.^[0:n-1])*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.
%        J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
%           Birkhauser, Basel, and Academic Press, New York, 1977, p. 154.

A = ones(n,1)*(1:n);
A = A./A';
A = tril(A) + tril(A,-1)';
End of lehmer.m
echo lesp.m 1>&2
cat >lesp.m <<'End of lesp.m'
function T = lesp(n)
%LESP   LESP(N) is an N-by-N tridiagonal matrix whose eigenvalues are real 
%       and smoothly distributed in the interval [-2*N-3.5, -4.5].
%       The sensitivities of the eigenvalues increase exponentially as
%       the eigenvalues grow more negative.
%       The matrix is similar to the symmetric tridiagonal matrix with the same
%       diagonal entries and with off-diagonal entries -1, via a similarity
%       transformation with D = diag(1!,2!,...,N!).

%       References:
%       H.W.J. Lenferink and M.N. Spijker, On the use of stability regions in
%            the numerical analysis of initial value problems, Report
%            TW-89-07, Department of Mathematics and Computer Science,
%            University of Leiden.
%       L.N. Trefethen, Psedospectra of matrices, Technical Report 91/10,
%            Oxford University Computing Laboratory, Oxford, 1991;
%            to appear in the Proceedings of the 14th Dundee conference.

x = 2:n;
T = tridiag( ones(x)./x, -(2*[x n+1]+1), x);
End of lesp.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 makejcf.m 1>&2
cat >makejcf.m <<'End of makejcf.m'
function A = makejcf(n, e, m, X)
%MAKEJCF    MAKEJCF(N, E, M) is a matrix having the Jordan canonical form
%           whose ith Jordan block is of dimension M(i) with eigenvalue E(i),
%           and where N = SUM(M).  The default is M = ONES(E).
%           The matrix is constructed by applying a random similarity
%           transformation to the Jordan form.
%           Alternatively, the matrix used in the similarity transformation
%           can be specified as a fifth parameter.
%           In particular, MAKEJCF(N, E, M, EYE(N)) returns the Jordan form
%           itself.
%           NB: The JCF is very sensitive to rounding errors.

if nargin < 3, m = ones(e); end

if any( size(e(:)) ~= size(m(:)) )
   error('Parameters E and M must be of same dimension.')
end

if sum(m) ~= n, error('Block dimensions must add up to N.'), end

A = zeros(n);
j = 1;
for i=1:max(size(m))
    if m(i) > 1
        Jb = jordan(m(i),e(i));
    else
        Jb = e(i);  % JORDAN fails in n=1 case.
    end
    A(j:j+m(i)-1,j:j+m(i)-1) = Jb;
    j = j + m(i);
end

if nargin < 4
   X = rand(n);
end
A = X\A*X;
End of makejcf.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,2] (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.
%           MATRIX with no arguments lists the names of the M-files in the
%           collection.

%  [1] N.J. Higham, A collection of test matrices in MATLAB, Technical Report
%      TR 89-1025, Department of Computer Science, Cornell University, Ithaca,
%      NY 14853, July 1989; also Numerical Analysis Report No. 172, Department
%      of Mathematics, University of Manchester, Manchester, M13 9PL, England,
%      1989.
%  [2] N.J. Higham, Algorithm 694: A collection of test matrices in MATLAB,
%      ACM Trans. Math. Soft., 17 (1991), pp. 289-305.

% 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 = 'augment ';

matrices = [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  '; 'lesp    '; 'lotkin  '; 'magic   '; 'minij   '; ...
'moler   '; 'ohess   '; 'orthog  '; 'parter  '; 'pascal  '; ...
'pdtoep  '; 'pei     '; 'prolate '; 'rando   '; 'randsvd ';];

matrices = [matrices
'riemann '; 'tridiag '; 'triw    '; 'vand    ';];

if nargin == 0

fprintf(' Test matrices:                                                  \n')
fprintf('                                                                 \n')
fprintf(' augment  condex   gallery  jordan   minij    prolate  wathen    \n')
fprintf(' cauchy   cycol    gear     kahan    moler    rando    wilk      \n')
fprintf(' chebspec dingdong gfpp     kms      ohess    randsvd            \n')
fprintf(' chebvand dorr     hadamard krylov   orthog   riemann            \n')
fprintf(' chow     dramadah hanowa   lauchli  parter   symm               \n')
fprintf(' circul   fiedler  hilb     lehmer   pascal   tridiag            \n')
fprintf(' clement  forsythe invol    lesp     pdtoep   triw               \n')
fprintf(' compan   frank    ipjfact  lotkin   pei      vand               \n')
fprintf('                                                                 \n')
fprintf(' Utility routines:                                               \n')
fprintf('                                                                 \n')
fprintf(' bandred  ps       show                                          \n')
fprintf(' chop     qmult    skew                                          \n')
fprintf(' comp     rq       sparsify                                      \n')
fprintf(' fv       see      sub                                           \n')
fprintf(' gersh    seesp                                                  \n')
fprintf(' house    seqa                                                   \n')
fprintf(' makejcf  seqcheb                                                \n')
fprintf(' matrix   seqm                                                   \n')

elseif nargin == 1
   if k == 0
      [A, temp] = size(matrices);
   elseif k > 0
      A = matrices(k,:);
   else
      % Version number and date of collection.
%     A = 'Version 1.0, July 4 1989';
%     A = 'Version 1.1, November 15 1989';
%     A = 'Version 1.2, May 30 1990';
      A = 'Version 1.3, November 14 1991';
   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.

%    References:
%    J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
%       Birkhauser, Basel, and Academic Press, New York, 1977, p. 158.
%    D.E. Rutherford, Some continuant determinants arising in physics and
%       chemistry---II, Proc. Royal Soc. Edin., 63, A (1952), pp. 232-241.
%       (For the eigenvalues of Givens' matrix.)

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, 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.  Q^4 is the identity.
%       K = 4:  Helmert matrix: a permutation of a lower Hessenberg matrix,
%               whose first row is ONES(1:N)/SQRT(N).
%       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).

%       References:
%       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.
%       P. Morton, On the eigenvectors of Schur's matrix, J. Number Theory,
%            12 (1980), pp. 122-127. (Re. ORTHOG(N, 3))
%       H.O. Lancaster, The Helmert Matrices, Amer. Math. Monthly, 72 (1965),
%            pp. 4-12.

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 == 4                          %  Helmert matrix

   Q = tril(ones(n));
   Q(1,2:n) = ones(1,n-1);
   for i=2:n
       Q(i,i) = -(i-1);
   end
   Q = diag( sqrt( [n 1:n-1] .* [1:n] ) ) \ Q;

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 parter.m 1>&2
cat >parter.m <<'End of parter.m'
function A = parter(n)
%PARTER    PARTER(N) is the matrix with (i,j) element 1/(i-j+0.5).
%          It is a Cauchy matrix and a Toeplitz matrix.
%          Most of the singular values cluster near PI.

%          At the Second SIAM Conference on Linear Algebra, Raleigh, N.C.,
%          1985, Cleve Moler noted that most of the singular values of
%          PARTER(N) are very close to PI.  An explanation of the phenomenon
%          was given by Parter; see also the paper by Tyrtyshnikov.

%          References:
%          The MathWorks Newsletter, Volume 1, Issue 1, March 1986, page 2.
%          S.V. Parter, On the distribution of the singular values of Toeplitz
%               matrices, Linear Algebra and Appl., 80 (1986), pp. 115-130.
%          E.E. Tyrtyshnikov, Cauchy-Toeplitz matrices and some applications,
%               Linear Algebra and Appl., 149 (1991), pp. 1-18.

A = cauchy( (1:n)+0.5, -(1:n) );
End of parter.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.

%    J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
%    Birkhauser, Basel, and Academic Press, New York, 1977, p. 172.

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 pdtoep.m 1>&2
cat >pdtoep.m <<'End of pdtoep.m'
function T = pdtoep(n, m, w, theta)
%PDTOEP   PDTOEP(N, M, W, THETA) is an N-by-N symmetric positive (semi-)
%         definite (SPD) Toeplitz matrix, comprised of the sum of M rank 2
%         (or, for certain THETA, rank 1) SPD Toeplitz matrices.
%         Specifically,
%                 T = W(1)*T(THETA(1)) + ... + W(M)*T(THETA(M)),
%         where T(THETA(k)) has (i,j) element COS(2*PI*THETA(k)*(i-j)).
%         Defaults: M = N, W = RAND(M,1), THETA = RAND(M,1) with
%                   RAND('UNIFORM').

%         Reference:
%         G. Cybenko and C.F. Van Loan, Computing the minimum eigenvalue of
%         a symmetric positive definite Toeplitz matrix, SIAM J. Sci. Stat.
%         Comput., 7 (1986), pp. 123-131.

if nargin < 2, m = n; end
if nargin < 3, rand('uniform'); w = rand(m,1); end
if nargin < 4, rand('uniform'); theta = rand(m,1); end

if max(size(w)) ~= m | max(size(theta)) ~= m
   error('Arguments W and THETA must be vectors of length M.')
end

T = zeros(n);
E = 2*pi*( (1:n)'*ones(1,n) - ones(n,1)*(1:n) );

for i=1:m
    T = T + w(i) * cos( theta(i)*E );
end
End of pdtoep.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 symmetric matrix
%       ALPHA*EYE(N) + ONES(N).
%       If ALPHA is omitted then ALPHA = 1 is used.
%       The matrix is 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 prolate.m 1>&2
cat >prolate.m <<'End of prolate.m'
function A = prolate(n, w)
%PROLATE   A = PROLATE(N, W) is the N-by-N prolate matrix with parameter W.
%          It is a symmetric Toeplitz matrix.
%          If 0 < W < 0.5 then
%             - A is positive definite
%             - the eigenvalues of A are distinct, lie in (0, 1), and
%               tend to cluster around 0 and 1.
%          W defaults to 0.25.

%  Reference: The Prolate matrix, talk given by by J.M. Varah
%  at the SIAM meeting on Signals, Systems, and Control held in
%  San Francisco, November 5-8, 1990.

if nargin == 1, w = 0.25; end

a = zeros(n,1);
a(1) = 2*w;
a(2:n) = sin( 2*pi*w*(1:n-1) ) ./ ( pi*(1:n-1) );

A = toeplitz(a);
End of prolate.m
echo ps.m 1>&2
cat >ps.m <<'End of ps.m'
function x = ps(A, m, tol, rl, noprint)
%PS     PS(A, M, TOL, RL) plots an approximation to the pseudospectrum
%       of the square matrix A, using M random perturbations of 2-norm TOL.
%       M defaults to 10 and TOL to 1E-3.
%       The perturbations are complex, unless RL is 1.
%       The eigenvalues of A are plotted as 'x'.
%       Alternative usage: X = PS(A, M, TOL, RL, 1) suppresses the plot and
%       returns the plot data in X.

%       For a given TOL, the pseudospectrum of A is the set of
%       pseudo-eigenvalues of A, that is, the set
%       { e : e is an eigenvalue of A+E, for some E with NORM(E) <= TOL }.

%       References:
%       L.N. Trefethen, Psedospectra of matrices, Technical Report 91/10, 
%            Oxford University Computing Laboratory, Oxford, 1991;
%            to appear in the Proceedings of the 14th Dundee conference.
%       L.N. Trefethen, Approximation theory and numerical linear algebra, in
%            Algorithms for Approximation II, J.C Mason and M.G. Cox, eds.,
%            Chapman and Hall, 1990, pp. 336-360

if nargin < 4, rl = 0; end
if nargin < 3, tol = 1e-3; end
if nargin < 2, m = 10; end

rand('normal');
   
n = max(size(A));
x = zeros(m*n,1);
i = sqrt(-1);

for j = 1:m
   if rl
      dA = rand(n);
   else
      dA = rand(n) + i*rand(n);
   end
   dA = tol/norm(dA)*dA;
   e = eig(A + dA);
   x((j-1)*n+1:j*n) = e;
end

if nargin < 5

   % Set axis ranges so both have the same length.
   xmin = min(real(x)); xmax = max(real(x));
   ymin = min(imag(x)); ymax = max(imag(x));
   if xmax-xmin >= ymax-ymin
      ymid = (ymin + ymax)/2;
      ymin =  ymid - (xmax-xmin)/2; ymax = ymid + (xmax-xmin)/2;
   else
      xmid = (xmin + xmax)/2;
      xmin = xmid - (ymax-ymin)/2; xmax = xmid + (ymax-ymin)/2;
   end
   axis('square')
   % Scale ranges by 1+2*alpha to give extra space around edges of plot.
   alpha = 0.1;
   axis( [xmin-alpha*(xmax-xmin) xmax+alpha*(xmax-xmin) ...
          ymin-alpha*(ymax-ymin) ymax+alpha*(ymax-ymin) ] );

   plot(real(x),imag(x),'.c5');
   hold on
   e = eig(A);
   plot(real(e),imag(e),'xw');
   hold off
   axis([1 2 3 4]); axis; axis('normal');  % Normal ratio and auto-scaling.

end
End of ps.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').

%       This routine is called by RANDSVD.

%       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.
%      Special case: if KAPPA < 0 then a random full symmetric positive
%                    definite matrix is produced with COND(A) = -KAPPA and
%                    eigenvalues distributed according to MODE.
%                    KL and KU, if present, are ignored.

%      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 abs(kappa) < 1, error('Condition number must be at least 1!'), end
posdef = 0; if kappa < 0, posdef = 1; kappa = -kappa; end  % Special case.

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    % In this case cond(A) <= kappa.
   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 posdef                % Handle special case.
   Q = qmult(p);
   A = Q'*sigma*Q;
   A = (A + A')/2;       % Ensure matrix is symmetric.
   return
end

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 readme.doc 1>&2
cat >readme.doc <<'End of readme.doc'
                Collection of Test Matrices in MATLAB

Note that in PLOT statements colours have been chosen as ones suitable
for a PC colour display, because these M-files were developed on a PC and
the default PLOT colours on a PC are rather dark.   These may not be
suitable for a display with a light background (most workstations), and in
this case PLOT statements in the following M-files should be suitably
modified: FV, GERSH, PS, SEE, SEESP.

Summary of changes for Version 1.3 (November 14, 1991).

New matrices: LESP, PARTER, PDTOEP, PROLATE.

New utility routines GERSH, MAKEJCF, PS.

Error corrected in diagonal scaling in KAHAN.  Default value of THETA changed
to 1.2, which gives less extreme ill-conditioning.

(1,2) element of k=2 case of CONDEX corrected.

Minor changes to comments of AUGMENT, CAUCHY, CLEMENT, DORR, HADAMARD, HILB,
RIEMANN.

SPARSE renamed to SEESP because a future release of MATLAB will contain a
routine of this name; see [3].

Bug in RANDSVD corrected: missing semi-colon in code for MODE<0.

In CHOP, default is now t=24, which gives the correct unit roundoff
for IEEE single precision.

SEE now prints a message if the matrix is not square so the user knows why no
field of values is displayed.

SEESP now uses 'x' for n<100,  and '.' for bigger dimensions.  It also uses
axis('square').

FV and PS share new code for ensuring that the plot has some space around the
edges and that the axes are on the same scale (previously the field of values
would be stretched at the whim of PLOT's automatic scaling).

Default THMAX in FV changed to 16 from 20.  Note that a larger value of THMAX
may be necessary for a correct plot, e.g. for the matrix HANOWA(8)^2.

---------
Summary of the changes in going from the original Version 1.0 (July 4, 1989)
[1] to Version 1.2 (May 30, 1990) [2].   Those marked with a `*' were already
(partially) present at Version 1.1.

(1*)  Some improvements to the documentation in the comment lines of the
M-files, including correction of errors, improved wording, and extra
references.

(2*)  A few minor bugs or inefficiencies corrected.

(3*)  Added a new matrix, bringing the total to 45:

%AUGMENT  AUGMENT(A) is the square matrix [EYE(m) A; A' ZEROS(n)] of dimension
%         m+n, where A is m-by-n.  It is the symmetric and indefinite
%         coefficient matrix of the augmented system associated with a least
%         squares problem minimize NORM(A*x-b).
%         Special case: if A is a scalar, n say, then AUGMENT(A) is the same
%                       as AUGMENT(RAND(p,q)) where n = p+q and p = ROUND(n/2),
%                       that is, a random augmented matrix of dimension n is
%                       produced.

(4*)  Added an extra option to ORTHOG:

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 = 4:  Helmert matrix: a permutation of a lower Hessenberg matrix,
%               whose first row is ONES(1:N)/SQRT(N).

(5*)  Changed second parameter in SPARSE from a character selecting the plot
type to a tolerance such that elements smaller than the tolerance are
regarded as zero.

(6)  Added extra functionality to RANDSVD:

function A = randsvd(n, kappa, mode, kl, ku)
%RANDSVD  Random matrices with pre-assigned singular values.
...
%      Special case: if KAPPA < 0 then a random full symmetric positive
%                    definite matrix is produced with COND(A) = -KAPPA and
%                    eigenvalues distributed according to MODE.
%                    KL and KU, if present, are ignored.

(7)  Changed SEE.  It now works for rectangular matrices.  Instead of the
inverse, the pseudo-inverse is plotted in the second window.  The plot of
columns of A versus the row index in the third window is replaced in this
new version by a plot of the singular values of A on a logarithmic scale.

(8)  HOUSE now avoids division by zero when it is given a zero vector
(BANDRED is the only routine that calls HOUSE).

---------
References

[1] N.J. Higham, A collection of test matrices in MATLAB, Technical Report
    TR 89-1025, Department of Computer Science, Cornell University, Ithaca,
    NY 14853, July 1989; also Numerical Analysis Report No. 172, Department
    of Mathematics, University of Manchester, Manchester, M13 9PL, England,
    1989.

[2] N.J. Higham, Algorithm 694: A collection of test matrices in MATLAB,
    ACM Trans. Math. Soft., 17 (1991), pp. 289-305.

[3] J.R. Gilbert, C.B. Moler and R.S. Schreiber, Sparse matrices in MATLAB:
    Design and implementation, Technical Report CSL-91-4, Xerox Palo Alto
    Research Center, Palo Alto, 1991; to appear in SIAM J. Matrix Anal. Appl.

   Nick Higham
   Dept. of Mathematics
   University of Manchester
   Manchester
   M13 9PL
   UK
   na.nhigham@na-net.ornl.gov
   November 30, 1991
End of readme.doc
echo riemann.m 1>&2
cat >riemann.m <<'End of riemann.m'
function A = riemann(n)
%RIEMANN    A = RIEMANN(N) is an 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).

%        This routine is called by FV.

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 (pseudo-) inverse.
%       SEE(A) displays MESH(A), MESH(PINV(A)), SEMILOGY(SVD(A),'x'),
%       and (if A is square) FV(A) in four subplot windows.
%       SEE(A,-1) plots just the eigenvalues in the last window, which is
%       much quicker than plotting the field of values.
%       SEE(A,1) uses the full screen for each picture, with a PAUSE in
%       between each one.

if nargin < 2, k = 0; end
[m, n] = size(A);
square = (m == n);

clg
B = pinv(A);
s = svd(A);
zs = (s == zeros(s));
if any( zs )
   s( zs ) = [];  % Remove zero singular values for semilogy plot.
end

if k == 1,
   mesh(A),                       pause
   mesh(B),                       pause
   semilogy(s, 'xg')
   hold on, semilogy(s, '--'), hold off, pause
   if any(zs), title('Zero(s) omitted'), end
else
   subplot(221)
   mesh(A)
   mesh(B)
   subplot(223)
   semilogy(s, 'xg')
   hold on, semilogy(s, '--'), hold off
   if any(zs), subplot(223), title('Zero(s) omitted'), subplot(224), end
end

if square
   if k == -1
      axis('square');
      e = eig(A);
      plot(real(e), imag(e), 'xg')
   else
      fv(A);
   end
else
   if k ==0
      subplot(224)
   else
      clg
   end
   xlabel('(Matrix not square)')
end
subplot
End of see.m
echo seesp.m 1>&2
cat >seesp.m <<'End of seesp.m'
function seesp(A, tol, x, y)
% SEESP   SEESP(A) plots the sparsity pattern of a matrix A, showing an
%         `x' or '.' (depending on SIZE(A)) for every nonzero element.
%         SEESP(A, TOL) plots the elements bigger than TOL in absolute value.
%         SEESP(A,TOL,x,y) overlays the current graphics screen taking (x,y)
%         as the origin.

if nargin <= 2
   x = 1; y = 1;
end
if nargin == 1, tol = 0; 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

if m*n >= 1e4
   s = 'w.';
else
   s = 'wx';
end
axis('square')
for j=n:-1:1
    plot( (j+y-1) .* [zeros(x-1,1); (abs(A(:,j))>tol)], s );
    hold on
end
hold off
axis('normal')
End of seesp.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, that is, with `+', `-' and
%       blank representing positive, negative and zero elements respectively.

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 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 is square and 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

% Is A square and symmetric?
symm = 0;
if min(size(A)) == max(size(A))
   if norm(A-A',1) == 0, symm = 1; end
end

rand('uniform')
if ~symm
   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 TRIDIAG(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,
%  Birkhauser, Basel, and Academic Press, New York, 1977, 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, SIAM J. Matrix Anal. Appl., 11 (1990),
%   pp. 23-41.

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
