########################################################################
##
## Copyright (C) 1989-2024 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or <https://octave.org/copyright/>.
##
## This file is part of Octave.
##
## Octave is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <https://www.gnu.org/licenses/>.
##
########################################################################

## -*- texinfo -*-
## @deftypefn  {} {} gallery (@var{name})
## @deftypefnx {} {} gallery (@var{name}, @var{args})
## Create interesting matrices for testing.
##
## @end deftypefn
##
## @deftypefn  {} {@var{c} =} gallery ("cauchy", @var{x})
## @deftypefnx {} {@var{c} =} gallery ("cauchy", @var{x}, @var{y})
## Create a Cauchy matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{c} =} gallery ("chebspec", @var{n})
## @deftypefnx {} {@var{c} =} gallery ("chebspec", @var{n}, @var{k})
## Create a Chebyshev spectral differentiation matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{c} =} gallery ("chebvand", @var{p})
## @deftypefnx {} {@var{c} =} gallery ("chebvand", @var{m}, @var{p})
## Create a @nospell{Vandermonde}-like matrix for the Chebyshev polynomials.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("chow", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("chow", @var{n}, @var{alpha})
## @deftypefnx {} {@var{a} =} gallery ("chow", @var{n}, @var{alpha}, @var{delta})
## Create a Chow matrix -- a singular Toeplitz lower Hessenberg matrix.
##
## @end deftypefn
##
## @deftypefn {} {@var{c} =} gallery ("circul", @var{v})
## Create a circulant matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("clement", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("clement", @var{n}, @var{k})
## Create a tridiagonal matrix with zero diagonal entries.
##
## @end deftypefn
##
## @deftypefn  {} {@var{c} =} gallery ("compar", @var{a})
## @deftypefnx {} {@var{c} =} gallery ("compar", @var{a}, @var{k})
## Create a comparison matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("condex", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("condex", @var{n}, @var{k})
## @deftypefnx {} {@var{a} =} gallery ("condex", @var{n}, @var{k}, @var{theta})
## Create a @nospell{"counterexample"} matrix to a condition estimator.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("cycol", [@var{m} @var{n}])
## @deftypefnx {} {@var{a} =} gallery ("cycol", @var{n})
## @deftypefnx {} {@var{a} =} gallery (@dots{}, @var{k})
## Create a matrix whose columns repeat cyclically.
##
## @end deftypefn
##
## @deftypefn  {} {[@var{c}, @var{d}, @var{e}] =} gallery ("dorr", @var{n})
## @deftypefnx {} {[@var{c}, @var{d}, @var{e}] =} gallery ("dorr", @var{n}, @var{theta})
## @deftypefnx {} {@var{a} =} gallery ("dorr", @dots{})
## Create a diagonally dominant, ill-conditioned, tridiagonal matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("dramadah", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("dramadah", @var{n}, @var{k})
## Create a (0, 1) matrix whose inverse has large integer entries.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("fiedler", @var{c})
## Create a symmetric @nospell{Fiedler} matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("forsythe", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("forsythe", @var{n}, @var{alpha})
## @deftypefnx {} {@var{a} =} gallery ("forsythe", @var{n}, @var{alpha}, @var{lambda})
## Create a @nospell{Forsythe} matrix (a perturbed Jordan block).
##
## @end deftypefn
##
## @deftypefn  {} {@var{f} =} gallery ("frank", @var{n})
## @deftypefnx {} {@var{f} =} gallery ("frank", @var{n}, @var{k})
## Create a Frank matrix (ill-conditioned eigenvalues).
##
## @end deftypefn
##
## @deftypefn {} {@var{c} =} gallery ("gcdmat", @var{n})
## Create a greatest common divisor matrix.
##
## @var{c} is an @var{n}-by-@var{n} matrix whose values correspond to the
## greatest common divisor of its coordinate values, i.e., @var{c}(i,j)
## correspond @code{gcd (i, j)}.
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("gearmat", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("gearmat", @var{n}, @var{i})
## @deftypefnx {} {@var{a} =} gallery ("gearmat", @var{n}, @var{i}, @var{j})
## Create a Gear matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{g} =} gallery ("grcar", @var{n})
## @deftypefnx {} {@var{g} =} gallery ("grcar", @var{n}, @var{k})
## Create a Toeplitz matrix with sensitive eigenvalues.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("hanowa", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("hanowa", @var{n}, @var{d})
## Create a matrix whose eigenvalues lie on a vertical line in the complex
## plane.
##
## @end deftypefn
##
## @deftypefn  {} {@var{v} =} gallery ("house", @var{x})
## @deftypefnx {} {[@var{v}, @var{beta}] =} gallery ("house", @var{x})
## Create a householder matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("integerdata", @var{imax}, [@var{M} @var{N} @dots{}], @var{j})
## @deftypefnx {} {@var{a} =} gallery ("integerdata", @var{imax}, @var{M}, @var{N}, @dots{}, @var{j})
## @deftypefnx {} {@var{a} =} gallery ("integerdata", [@var{imin}, @var{imax}], [@var{M} @var{N} @dots{}], @var{j})
## @deftypefnx {} {@var{a} =} gallery ("integerdata", [@var{imin}, @var{imax}], @var{M}, @var{N}, @dots{}, @var{j})
## @deftypefnx {} {@var{a} =} gallery ("integerdata", @dots{}, "@var{class}")
## Create a matrix with random integers in the range [1, @var{imax}].
## If @var{imin} is given then the integers are in the range
## [@var{imin}, @var{imax}].
##
## The second input is a matrix of dimensions describing the size of the
## output.  The dimensions can also be input as comma-separated arguments.
##
## The input @var{j} is an integer index in the range [0, 2^32-1].  The values
## of the output matrix are always exactly the same (reproducibility) for a
## given size input and @var{j} index.
##
## The final optional argument determines the class of the resulting matrix.
## Possible values for @var{class}: @qcode{"uint8"}, @qcode{"uint16"},
## @qcode{"uint32"}, @qcode{"int8"}, @qcode{"int16"}, int32", @qcode{"single"},
## @qcode{"double"}.  The default is @qcode{"double"}.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("invhess", @var{x})
## @deftypefnx {} {@var{a} =} gallery ("invhess", @var{x}, @var{y})
## Create the inverse of an upper Hessenberg matrix.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("invol", @var{n})
## Create an involutory matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("ipjfact", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("ipjfact", @var{n}, @var{k})
## Create a Hankel matrix with factorial elements.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("jordbloc", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("jordbloc", @var{n}, @var{lambda})
## Create a Jordan block.
##
## @end deftypefn
##
## @deftypefn  {} {@var{u} =} gallery ("kahan", @var{n})
## @deftypefnx {} {@var{u} =} gallery ("kahan", @var{n}, @var{theta})
## @deftypefnx {} {@var{u} =} gallery ("kahan", @var{n}, @var{theta}, @var{pert})
## Create a @nospell{Kahan} matrix (upper trapezoidal).
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("kms", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("kms", @var{n}, @var{rho})
## Create a @nospell{Kac-Murdock-Szego} Toeplitz matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{b} =} gallery ("krylov", @var{a})
## @deftypefnx {} {@var{b} =} gallery ("krylov", @var{a}, @var{x})
## @deftypefnx {} {@var{b} =} gallery ("krylov", @var{a}, @var{x}, @var{j})
## Create a Krylov matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("lauchli", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("lauchli", @var{n}, @var{mu})
## Create a @nospell{Lauchli} matrix (rectangular).
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("lehmer", @var{n})
## Create a @nospell{Lehmer} matrix (symmetric positive definite).
##
## @end deftypefn
##
## @deftypefn {} {@var{t} =} gallery ("lesp", @var{n})
## Create a tridiagonal matrix with real, sensitive eigenvalues.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("lotkin", @var{n})
## Create a @nospell{Lotkin} matrix.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("minij", @var{n})
## Create a symmetric positive definite matrix MIN(i,j).
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("moler", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("moler", @var{n}, @var{alpha})
## Create a @nospell{Moler} matrix (symmetric positive definite).
##
## @end deftypefn
##
## @deftypefn {} {[@var{a}, @var{t}] =} gallery ("neumann", @var{n})
## Create a singular matrix from the discrete @nospell{Neumann} problem
## (sparse).
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("normaldata", [@var{M} @var{N} @dots{}], @var{j})
## @deftypefnx {} {@var{a} =} gallery ("normaldata", @var{M}, @var{N}, @dots{}, @var{j})
## @deftypefnx {} {@var{a} =} gallery ("normaldata", @dots{}, "@var{class}")
## Create a matrix with random samples from the standard normal distribution
## (mean = 0, std = 1).
##
## The first input is a matrix of dimensions describing the size of the output.
## The dimensions can also be input as comma-separated arguments.
##
## The input @var{j} is an integer index in the range [0, 2^32-1].  The values
## of the output matrix are always exactly the same (reproducibility) for a
## given size input and @var{j} index.
##
## The final optional argument determines the class of the resulting matrix.
## Possible values for @var{class}: @qcode{"single"}, @qcode{"double"}.
## The default is @qcode{"double"}.
##
## @end deftypefn
##
## @deftypefn  {} {@var{q} =} gallery ("orthog", @var{n})
## @deftypefnx {} {@var{q} =} gallery ("orthog", @var{n}, @var{k})
## Create orthogonal and nearly orthogonal matrices.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("parter", @var{n})
## Create a @nospell{Parter} matrix (a Toeplitz matrix with singular values
## near pi).
##
## @end deftypefn
##
## @deftypefn  {} {@var{p} =} gallery ("pei", @var{n})
## @deftypefnx {} {@var{p} =} gallery ("pei", @var{n}, @var{alpha})
## Create a Pei matrix.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("poisson", @var{n})
## Create a block tridiagonal matrix from Poisson's equation (sparse).
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("prolate", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("prolate", @var{n}, @var{w})
## Create a prolate matrix (symmetric, ill-conditioned Toeplitz matrix).
##
## @end deftypefn
##
## @deftypefn {} {@var{h} =} gallery ("randhess", @var{x})
## Create a random, orthogonal upper Hessenberg matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("rando", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("rando", @var{n}, @var{k})
## Create a random matrix with elements -1, 0 or 1.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("randsvd", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("randsvd", @var{n}, @var{kappa})
## @deftypefnx {} {@var{a} =} gallery ("randsvd", @var{n}, @var{kappa}, @var{mode})
## @deftypefnx {} {@var{a} =} gallery ("randsvd", @var{n}, @var{kappa}, @var{mode}, @var{kl})
## @deftypefnx {} {@var{a} =} gallery ("randsvd", @var{n}, @var{kappa}, @var{mode}, @var{kl}, @var{ku})
## Create a random matrix with pre-assigned singular values.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("redheff", @var{n})
## Create a zero and ones matrix of @nospell{Redheffer} associated with the
## Riemann hypothesis.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("riemann", @var{n})
## Create a matrix associated with the Riemann hypothesis.
##
## @end deftypefn
##
## @deftypefn {} {@var{a} =} gallery ("ris", @var{n})
## Create a symmetric Hankel matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("smoke", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("smoke", @var{n}, @var{k})
## Create a complex matrix, with a @nospell{"smoke ring"} pseudospectrum.
##
## @end deftypefn
##
## @deftypefn  {} {@var{t} =} gallery ("toeppd", @var{n})
## @deftypefnx {} {@var{t} =} gallery ("toeppd", @var{n}, @var{m})
## @deftypefnx {} {@var{t} =} gallery ("toeppd", @var{n}, @var{m}, @var{w})
## @deftypefnx {} {@var{t} =} gallery ("toeppd", @var{n}, @var{m}, @var{w}, @var{theta})
## Create a symmetric positive definite Toeplitz matrix.
##
## @end deftypefn
##
## @deftypefn  {} {@var{p} =} gallery ("toeppen", @var{n})
## @deftypefnx {} {@var{p} =} gallery ("toeppen", @var{n}, @var{a})
## @deftypefnx {} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}, @var{b})
## @deftypefnx {} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}, @var{b}, @var{c})
## @deftypefnx {} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}, @var{b}, @var{c}, @var{d})
## @deftypefnx {} {@var{p} =} gallery ("toeppen", @var{n}, @var{a}, @var{b}, @var{c}, @var{d}, @var{e})
## Create a pentadiagonal Toeplitz matrix (sparse).
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("tridiag", @var{x}, @var{y}, @var{z})
## @deftypefnx {} {@var{a} =} gallery ("tridiag", @var{n})
## @deftypefnx {} {@var{a} =} gallery ("tridiag", @var{n}, @var{c}, @var{d}, @var{e})
## Create a tridiagonal matrix (sparse).
##
## @end deftypefn
##
## @deftypefn  {} {@var{t} =} gallery ("triw", @var{n})
## @deftypefnx {} {@var{t} =} gallery ("triw", @var{n}, @var{alpha})
## @deftypefnx {} {@var{t} =} gallery ("triw", @var{n}, @var{alpha}, @var{k})
## Create an upper triangular matrix discussed by
## @nospell{Kahan, Golub, and Wilkinson}.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("uniformdata", [@var{M} @var{N} @dots{}], @var{j})
## @deftypefnx {} {@var{a} =} gallery ("uniformdata", @var{M}, @var{N}, @dots{}, @var{j})
## @deftypefnx {} {@var{a} =} gallery ("uniformdata", @dots{}, "@var{class}")
## Create a matrix with random samples from the standard uniform distribution
## (range [0,1]).
##
## The first input is a matrix of dimensions describing the size of the output.
## The dimensions can also be input as comma-separated arguments.
##
## The input @var{j} is an integer index in the range [0, 2^32-1].  The values
## of the output matrix are always exactly the same (reproducibility) for a
## given size input and @var{j} index.
##
## The final optional argument determines the class of the resulting matrix.
## Possible values for @var{class}: @qcode{"single"}, @qcode{"double"}.
## The default is @qcode{"double"}.
##
## @end deftypefn
##
## @deftypefn  {} {@var{a} =} gallery ("wathen", @var{nx}, @var{ny})
## @deftypefnx {} {@var{a} =} gallery ("wathen", @var{nx}, @var{ny}, @var{k})
## Create the @nospell{Wathen} matrix.
##
## @end deftypefn
##
## @deftypefn {} {[@var{a}, @var{b}] =} gallery ("wilk", @var{n})
## Create various specific matrices devised/discussed by Wilkinson.
##
## @end deftypefn

## Code for most of the individual matrices (except binomial, gcdmat,
## integerdata, leslie, normaldata, randcolu, randcorr, randjorth, sampling,
## uniformdata) by Nicholas J. Higham <Nicholas.J.Higham@manchester.ac.uk>
## Adapted for Octave and into single gallery function by Carnë Draug

function varargout = gallery (name, varargin)

  if (nargin < 1)
    print_usage ();
  elseif (! ischar (name))
    error ("gallery: NAME must be a string");
  endif

  ## NOTE: there isn't a lot of input check in the individual functions
  ## that actually build the functions.  This is by design. The original
  ## code by Higham did not perform it and was propagated to Matlab, so
  ## for compatibility, we also don't make it. For example, arguments
  ## that behave as switches, and in theory accepting a value of 0 or 1,
  ## will use a value of 0, for any value other than 1 (only check made
  ## is if the value is equal to 1). It will often also accept string
  ## values instead of numeric. Only input check added was where it
  ## would be causing an error anyway.

  ## we will always want to return at least 1 output
  n_out = nargout;
  if (n_out == 0)
    n_out = 1;
  endif

  switch (lower (name))
    case "binomial"
      error ("gallery: matrix %s not implemented", name);
    case "cauchy"      , [varargout{1:n_out}] = cauchy      (varargin{:});
    case "chebspec"    , [varargout{1:n_out}] = chebspec    (varargin{:});
    case "chebvand"    , [varargout{1:n_out}] = chebvand    (varargin{:});
    case "chow"        , [varargout{1:n_out}] = chow        (varargin{:});
    case "circul"      , [varargout{1:n_out}] = circul      (varargin{:});
    case "clement"     , [varargout{1:n_out}] = clement     (varargin{:});
    case "compar"      , [varargout{1:n_out}] = compar      (varargin{:});
    case "condex"      , [varargout{1:n_out}] = condex      (varargin{:});
    case "cycol"       , [varargout{1:n_out}] = cycol       (varargin{:});
    case "dorr"        , [varargout{1:n_out}] = dorr        (varargin{:});
    case "dramadah"    , [varargout{1:n_out}] = dramadah    (varargin{:});
    case "fiedler"     , [varargout{1:n_out}] = fiedler     (varargin{:});
    case "forsythe"    , [varargout{1:n_out}] = forsythe    (varargin{:});
    case "frank"       , [varargout{1:n_out}] = frank       (varargin{:});
    case "gearmat"     , [varargout{1:n_out}] = gearmat     (varargin{:});
    case "gcdmat"      , [varargout{1:n_out}] = gcdmat      (varargin{:});
    case "grcar"       , [varargout{1:n_out}] = grcar       (varargin{:});
    case "hanowa"      , [varargout{1:n_out}] = hanowa      (varargin{:});
    case "house"       , [varargout{1:n_out}] = house       (varargin{:});
    case "integerdata" , [varargout{1:n_out}] = integerdata (varargin{:});
    case "invhess"     , [varargout{1:n_out}] = invhess     (varargin{:});
    case "invol"       , [varargout{1:n_out}] = invol       (varargin{:});
    case "ipjfact"     , [varargout{1:n_out}] = ipjfact     (varargin{:});
    case "jordbloc"    , [varargout{1:n_out}] = jordbloc    (varargin{:});
    case "kahan"       , [varargout{1:n_out}] = kahan       (varargin{:});
    case "kms"         , [varargout{1:n_out}] = kms         (varargin{:});
    case "krylov"      , [varargout{1:n_out}] = krylov      (varargin{:});
    case "lauchli"     , [varargout{1:n_out}] = lauchli     (varargin{:});
    case "lehmer"      , [varargout{1:n_out}] = lehmer      (varargin{:});
    case "leslie"
      error ("gallery: matrix %s not implemented", name);
    case "lesp"        , [varargout{1:n_out}] = lesp        (varargin{:});
    case "lotkin"      , [varargout{1:n_out}] = lotkin      (varargin{:});
    case "minij"       , [varargout{1:n_out}] = minij       (varargin{:});
    case "moler"       , [varargout{1:n_out}] = moler       (varargin{:});
    case "neumann"     , [varargout{1:n_out}] = neumann     (varargin{:});
    case "normaldata"  , [varargout{1:n_out}] = normaldata  (varargin{:});
    case "orthog"      , [varargout{1:n_out}] = orthog      (varargin{:});
    case "parter"      , [varargout{1:n_out}] = parter      (varargin{:});
    case "pei"         , [varargout{1:n_out}] = pei         (varargin{:});
    case "poisson"     , [varargout{1:n_out}] = poisson     (varargin{:});
    case "prolate"     , [varargout{1:n_out}] = prolate     (varargin{:});
    case "randcolu"
      error ("gallery: matrix %s not implemented", name);
    case "randcorr"
      error ("gallery: matrix %s not implemented", name);
    case "randhess"    , [varargout{1:n_out}] = randhess    (varargin{:});
    case "randjorth"
      error ("gallery: matrix %s not implemented", name);
    case "rando"       , [varargout{1:n_out}] = rando       (varargin{:});
    case "randsvd"     , [varargout{1:n_out}] = randsvd     (varargin{:});
    case "redheff"     , [varargout{1:n_out}] = redheff     (varargin{:});
    case "riemann"     , [varargout{1:n_out}] = riemann     (varargin{:});
    case "ris"         , [varargout{1:n_out}] = ris         (varargin{:});
    case "sampling"
      error ("gallery: matrix %s not implemented", name);
    case "smoke"       , [varargout{1:n_out}] = smoke       (varargin{:});
    case "toeppd"      , [varargout{1:n_out}] = toeppd      (varargin{:});
    case "toeppen"     , [varargout{1:n_out}] = toeppen     (varargin{:});
    case "tridiag"     , [varargout{1:n_out}] = tridiag     (varargin{:});
    case "triw"        , [varargout{1:n_out}] = triw        (varargin{:});
    case "uniformdata" , [varargout{1:n_out}] = uniformdata (varargin{:});
    case "wathen"      , [varargout{1:n_out}] = wathen      (varargin{:});
    case "wilk"        , [varargout{1:n_out}] = wilk        (varargin{:});
    otherwise
      error ("gallery: unknown matrix with NAME %s", name);
  endswitch

endfunction

function C = cauchy (x, y)

  ## CAUCHY  Cauchy matrix.
  ##   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).
  ##   C is totally positive if 0 < X(1) < ... < X(N) and
  ##   0 < Y(1) < ... < Y(N).
  ##
  ##   References:
  ##   N.J. Higham, Accuracy and Stability of Numerical Algorithms,
  ##     Society for Industrial and Applied Mathematics, Philadelphia, PA,
  ##     USA, 1996; sec. 26.1.
  ##   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.
  ##     O. Taussky and M. Marcus, Eigenvalues of finite matrices, in
  ##     Survey of Numerical Analysis, J. Todd, ed., McGraw-Hill, New York,
  ##     pp. 279-313, 1962. (States the totally positive property on p. 295.)

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 or 2 arguments are required for cauchy matrix");
  elseif (! isnumeric (x))
    error ("gallery: X must be numeric for cauchy matrix");
  elseif (nargin == 2 && ! isnumeric (y))
    error ("gallery: Y must be numeric for cauchy matrix");
  endif

  n = numel (x);
  if (isscalar (x) && fix (x) == x)
    n = x;
    x = 1:n;
  elseif (n > 1 && isvector (x))
    ## do nothing
  else
    error ("gallery: X be an integer or a vector for cauchy matrix");
  endif

  if (nargin == 1)
    y = x;
  endif

  ## Ensure x and y are column vectors
  x = x(:);
  y = y(:);
  if (numel (x) != numel (y))
    error ("gallery: X and Y must be vectors of same length for cauchy matrix");
  endif

  C = 1 ./ (x + y.');

endfunction

function C = chebspec (n, k = 0)

  ## CHEBSPEC  Chebyshev spectral differentiation matrix.
  ##   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(REAL(X)) is interesting).
  ##
  ##   References:
  ##   C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral
  ##      Methods in Fluid Dynamics, Springer-Verlag, Berlin, 1988; p. 69.
  ##   L.N. Trefethen and 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 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for chebspec matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for chebspec matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a scalar for chebspec matrix");
  endif

  ## k = 1 case obtained from k = 0 case with one bigger n.
  switch (k)
    case 0, # do nothing
    case 1, n = n + 1;
    otherwise
      error ("gallery: K should be either 0 or 1 for chebspec matrix");
  endswitch

  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;

  ## eye(size(C)) on next line avoids div by zero.
  C = (d * (one./d)') ./ (x*one'-one*x' + eye (size (C)));

  ##  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,:);
    endif
    if (i < n+1)
      C(i,i) = -x(i)/(2*(1-x(i)^2));
    else
      C(n+1,n+1) = -C(1,1);
    endif
  endfor

  if (k == 1)
    C = C(2:n+1,2:n+1);
  endif

endfunction

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 || nargin > 2)
    error ("gallery: 1 or 2 arguments are required for chebvand matrix");
  endif

  ## because the order of the arguments changes if nargin is 1 or 2 ...

  if (nargin == 1)
    p = m;
  endif

  n = numel (p);
  if (! isnumeric (p))
    error ("gallery: P must be numeric for chebvand matrix");
  elseif (isscalar (p) && fix (p) == p)
    n = p;
    p = linspace (0, 1, n);
  elseif (n > 1 && isvector (p))
    ## do nothing
  endif
  p = p(:).'; # Ensure p is a row vector.

  if (nargin == 1)
    m = n;
  elseif (! isnumeric (m) || ! isscalar (m))
    error ("gallery: M must be a scalar for chebvand matrix");
  endif

  C = ones (m, n);
  if (m != 1)
    C(2,:) = p;
    ##      Use Chebyshev polynomial recurrence.
    for i = 3:m
      C(i,:) = 2.*p.*C(i-1,:) - C(i-2,:);
    endfor
  endif

endfunction

function A = chow (n, alpha = 1, delta = 0)

  ## CHOW    Chow matrix - a singular Toeplitz lower Hessenberg matrix.
  ##   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/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 < 1 || nargin > 3)
    error ("gallery: 1 to 3 arguments are required for chow matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for chow matrix");
  elseif (! isnumeric (alpha) || ! isscalar (alpha))
    error ("gallery: ALPHA must be a scalar for chow matrix");
  elseif (! isnumeric (delta) || ! isscalar (delta))
    error ("gallery: DELTA must be a scalar for chow matrix");
  endif

  A = toeplitz (alpha.^(1:n), [alpha 1 zeros(1, n-2)]) + delta * eye (n);

endfunction

function C = circul (v)

  ## CIRCUL  Circulant matrix.
  ##   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.

  if (nargin != 1)
    error ("gallery: 1 argument is required for circul matrix");
  elseif (! isnumeric (v))
    error ("gallery: V must be numeric for circul matrix");
  endif

  n = numel (v);
  if (isscalar (v) && fix (v) == v)
    n = v;
    v = 1:n;
  elseif (n > 1 && isvector (v))
    ## do nothing
  else
    error ("gallery: X must be a scalar or a vector for circul matrix");
  endif

  v = v(:).';   # Make sure v is a row vector
  C = toeplitz ([v(1) v(n:-1:2)], v);

endfunction

function A = clement (n, k = 0)

  ## CLEMENT   Clement matrix - tridiagonal with zero diagonal entries.
  ##   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.
  ##   A. Edelman and E. Kostlan, The road from Kac's matrix to Kac's
  ##      random polynomials. In John~G. Lewis, editor, Proceedings of
  ##      the Fifth SIAM Conference on Applied Linear Algebra Society
  ##      for Industrial and Applied Mathematics, Philadelphia, 1994,
  ##      pp. 503-507.
  ##   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 || nargin > 2)
    error ("gallery: 1 or 2 arguments are required for clement matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for clement matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for clement matrix");
  endif

  n -= 1;
  x = n:-1:1;
  z = 1:n;

  if (k == 0)
    A = diag (x, -1) + diag (z, 1);
  elseif (k == 1)
    y = sqrt (x.*z);
    A = diag (y, -1) + diag (y, 1);
  else
    error ("gallery: K must have a value of 0 or 1 for clement matrix");
  endif

endfunction

function C = compar (A, k = 0)

  ## 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 || nargin > 2)
    error ("gallery: 1 or 2 arguments are required for compar matrix");
  elseif (! isnumeric (A) || ndims (A) != 2)
    error ("gallery: A must be a 2-D matrix for compar matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for compar matrix");
  endif

  [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));
    endfor

  elseif (k == 1)
    C = A';
    for j = 1:p
      C(k,k) = 0;
    endfor
    mx = max (abs (C));
    C  = -mx'*ones (1, n);
    for j = 1:p
      C(j,j) = abs (A(j,j));
    endfor
    if (all (A == tril (A))), C = tril (C); endif
    if (all (A == triu (A))), C = triu (C); endif

  else
    error ("gallery: K must have a value of 0 or 1 for compar matrix");
  endif

endfunction

function A = condex (n, k = 4, theta = 100)

  ## CONDEX   'Counterexamples' to matrix condition number estimators.
  ##   CONDEX(N, K, THETA) is a 'counterexample' 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) (independent of THETA)
  ##       K = 4:   N >= 4,     SONEST (Higham 1988)
  ##   (Note that in practice the K = 4 matrix is not usually a
  ##    counterexample 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
  ##      (Algorithm 674), ACM Trans. Math. Soft., 14 (1988), pp. 381-396.

  if (nargin < 1 || nargin > 3)
    error ("gallery: 1 to 3 arguments are required for condex matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for condex matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for condex matrix");
  elseif (! isnumeric (theta) || ! isscalar (theta))
    error ("gallery: THETA must be a numeric scalar for condex matrix");
  endif

  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 = gallery ("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;

  else
    error ("gallery: unknown estimator K '%d' for condex matrix", k);
  endif

  ## Pad out with identity as necessary.
  m = columns (A);
  if (m < n)
    for i = n:-1:m+1
      A(i,i) = 1;
    endfor
  endif

endfunction

function A = cycol (n, k = max (round (n(end)/4), 1))

  ## CYCOL   Matrix whose columns repeat cyclically.
  ##   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 = RANDN(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).

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 or 2 arguments are required for cycol matrix");
  elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n)
    error ("gallery: N must be a 1 or 2 element integer for cycol matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a scalar for cycol matrix");
  endif

  ## Parameter n specifies dimension: m-by-n
  m = n(1);
  n = n(end);

  A = randn (m, k);
  for i = 2:ceil (n/k)
    A = [A A(:,1:k)];
  endfor
  A = A(:,1:n);

endfunction

function [c, d, e] = dorr (n, theta = 0.01)

  ## DORR  Dorr matrix - diagonally dominant, ill conditioned, tridiagonal.
  ##   [C, D, E] = DORR(N, THETA) returns the vectors defining a row diagonally
  ##   dominant, tridiagonal M-matrix that is ill conditioned for small
  ##   values of the parameter THETA >= 0.
  ##   If only one output parameter is supplied then
  ##   C = FULL(TRIDIAG(C,D,E)), i.e., the matrix iself 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 < 1 || nargin > 2)
    error ("gallery: 1 or 2 arguments are required for dorr matrix");
  elseif (! isscalar (n) || ! isnumeric (n) || fix (n) != n)
    error ("gallery: N must be an integer for dorr matrix");
  elseif (! isscalar (theta) || ! isnumeric (theta))
    error ("gallery: THETA must be a numeric scalar for dorr matrix");
  endif

  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 (m, 1);
  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);
  endif

endfunction

function A = dramadah (n, k = 1)

  ## DRAMADAH  A (0,1) matrix whose inverse has large integer entries.
  ##   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.
  ##
  ##   Another interesting (0,1) matrix:
  ##   K = 3: A has maximal determinant among (0,1) lower Hessenberg
  ##   matrices: det(A) = the n'th Fibonacci number.  A is Toeplitz.
  ##   The eigenvalues have an interesting distribution in the complex
  ##   plane.
  ##
  ##   References:
  ##   R.L. Graham and N.J.A. Sloane, Anti-Hadamard matrices,
  ##      Linear Algebra and Appl., 62 (1984), pp. 113-137.
  ##   L. Ching, The maximum determinant of an nxn lower Hessenberg
  ##      (0,1) matrix, Linear Algebra and Appl., 183 (1993), pp. 147-153.

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for dramadah matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for dramadah matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for dramadah matrix");
  endif

  switch (k)
    case (1)  # Toeplitz
      c = ones (n, 1);
      for i = 2:4:n
        m = min (1, n-i);
        c(i:i+m) = zeros (m+1, 1);
      endfor
      r = zeros (n, 1);
      r(1:4) = [1 1 0 1];
      if (n < 4)
        r = r(1:n);
      endif
      A = toeplitz (c, r);

    case (2)  # Upper triangular and Toeplitz
      c = zeros (n, 1);
      c(1) = 1;
      r = ones (n, 1);
      for i= 3:2:n
        r(i) = 0;
      endfor
      A = toeplitz (c, r);

    case (3)  # Lower Hessenberg
      c = ones (n, 1);
      for i= 2:2:n
        c(i) = 0;
      endfor
      A = toeplitz (c, [1 1 zeros(1,n-2)]);

    otherwise
      error ("gallery: unknown K '%d' for dramadah matrix", k);
  endswitch

endfunction

function A = fiedler (c)

  ## FIEDLER  Fiedler matrix - symmetric.
  ##   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 by 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.

  if (nargin != 1)
    error ("gallery: 1 argument is required for fiedler matrix");
  elseif (! isnumeric (c))
    error ("gallery: C must be numeric for fiedler matrix");
  endif

  n = numel (c);
  if (isscalar (c) && fix (c) == c)
    n = c;
    c = 1:n;
  elseif (n > 1 && isvector (c))
    ## do nothing
  else
    error ("gallery: C must be an integer or a vector for fiedler matrix");
  endif
  c = c(:).';           # Ensure c is a row vector.

  A = abs (c - c.');

endfunction

function A = forsythe (n, alpha = sqrt (eps), lambda = 0)

  ## FORSYTHE  Forsythe matrix - a perturbed Jordan block.
  ##   FORSYTHE(N, ALPHA, LAMBDA) is the N-by-N matrix equal to
  ##   JORDBLOC(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 < 1 || nargin > 3)
    error ("gallery: 1 to 3 arguments are required for forsythe matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for forsythe matrix");
  elseif (! isnumeric (alpha) || ! isscalar (alpha))
    error ("gallery: ALPHA must be a numeric scalar for forsythe matrix");
  elseif (! isnumeric (lambda) || ! isscalar (lambda))
    error ("gallery: LAMBDA must be a numeric scalar for forsythe matrix");
  endif

  A = jordbloc (n, lambda);
  A(n,1) = alpha;

endfunction

function F = frank (n, k = 0)

  ## FRANK   Frank matrix---ill conditioned eigenvalues.
  ##   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).
  ##   F has all positive eigenvalues and they occur in reciprocal pairs
  ##   (so that 1 is an eigenvalue if N is odd).
  ##   The eigenvalues of F may be obtained in terms of the zeros of the
  ##   Hermite polynomials.
  ##   The FLOOR(N/2) smallest eigenvalues of F are ill conditioned,
  ##   the more so for bigger N.
  ##
  ##   DET(FRANK(N)') comes out far from 1 for large N---see Frank (1958)
  ##   and Wilkinson (1960) for discussions.
  ##
  ##   This version incorporates improvements suggested by W. Kahan.
  ##
  ##   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).
  ##   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).
  ##   H. Rutishauser, On test matrices, Programmation en Mathematiques
  ##      Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165,
  ##      1966, pp. 349-365.  Section 9.
  ##   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).
  ##   The next two references give details of the eigensystem, as does
  ##   Rutishauser (see above).
  ##   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 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for frank matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for frank matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for frank matrix");
  endif

  p = n:-1:1;
  F = triu (p(ones (n, 1), :) - diag (ones (n-1, 1), -1), -1);

  switch (k)
    case 0, # do nothing
    case 1, F = F(p,p)';
    otherwise
      error ("gallery: K must have a value of 0 or 1 for frank matrix");
  endswitch

endfunction

function c = gcdmat (n)

  if (nargin != 1)
    error ("gallery: 1 argument is required for gcdmat matrix");
  elseif (! isscalar (n) || ! isnumeric (n) || fix (n) != n)
    error ("gallery: N must be an integer for gcdmat matrix");
  endif
  c = gcd (repmat ((1:n)', [1 n]), repmat (1:n, [n 1]));

endfunction

function A = gearmat (n, i = n, j = -n)

  ## NOTE: this function was named gearm in the original Test Matrix Toolbox
  ## GEARMAT   Gear matrix.
  ##   A = GEARMAT(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.
  ##   GEARMAT(N) is singular.
  ##
  ##   (GEAR is a Simulink function, hence GEARMAT for Gear matrix.)
  ##   Reference:
  ##   C.W. Gear, A simple set of test matrices for eigenvalue programs,
  ##   Math. Comp., 23 (1969), pp. 119-125.

  if (nargin < 1 || nargin > 3)
    error ("gallery: 1 to 3 arguments are required for gearmat matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for gearmat matrix");
  elseif (! isnumeric (i) || ! isscalar (i) || i == 0 || abs (i) > n)
    error ("gallery: I must be a nonzero scalar, and abs (I) <= N for gearmat matrix");
  elseif (! isnumeric (j) || ! isscalar (j) || i == 0 || abs (j) > n)
    error ("gallery: J must be a nonzero scalar, and abs (J) <= N for gearmat matrix");
  endif

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

endfunction

function G = grcar (n, k = 3)

  ## GRCAR     Grcar matrix - a Toeplitz matrix with sensitive eigenvalues.
  ##   GRCAR(N, K) is an N-by-N matrix with -1s on the
  ##   subdiagonal, 1s on the diagonal, and K superdiagonals of 1s.
  ##   The default is K = 3.  The eigenvalues of this matrix form an
  ##   interesting pattern in the complex plane (try PS(GRCAR(32))).
  ##
  ##   References:
  ##   J.F. Grcar, Operator coefficient methods for linear equations,
  ##        Report SAND89-8691, Sandia National Laboratories, Albuquerque,
  ##        New Mexico, 1989 (Appendix 2).
  ##   N.M. Nachtigal, L. Reichel and L.N. Trefethen, A hybrid GMRES
  ##        algorithm for nonsymmetric linear systems, SIAM J. Matrix Anal.
  ##        Appl., 13 (1992), pp. 796-825.

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for grcar matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for grcar matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for grcar matrix");
  endif

  G = tril (triu (ones (n)), k) - diag (ones (n-1, 1), -1);

endfunction

function A = hanowa (n, d = -1)

  ## HANOWA  A matrix whose eigenvalues lie on a vertical line in the complex plane.
  ##   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 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for hanowa matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for hanowa matrix");
  elseif (rem (n, 2) != 0)
    error ("gallery: N must be even for hanowa matrix");
  elseif (! isnumeric (d) || ! isscalar (d))
    error ("gallery: D must be a numeric scalar for hanowa matrix");
  endif

  m = n/2;
  A = [ d*eye(m)  -diag(1:m)
        diag(1:m)  d*eye(m) ];

endfunction

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 1989, 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)|).
  ##
  ##   References:
  ##   G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  ##      Johns Hopkins University Press, Baltimore, Maryland, 1989.
  ##   G.W. Stewart, Introduction to Matrix Computations, Academic Press,
  ##      New York, 1973,
  ##   J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
  ##      Press, 1965.

  if (nargin != 1)
    error ("gallery: 1 argument is required for house matrix");
  elseif (! isnumeric (x) || ! isvector (x))
    error ("gallery: X must be a vector for house matrix");
  endif

  ## must be a column vector
  x = x(:);

  s = norm (x) * (sign (x(1)) + (x(1) == 0)); # Modification for sign (0) == 1.
  v = x;
  if (s == 0)
    ## Quit if x is the zero vector.
    beta = 1;
  else
    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.
  endif

endfunction

function A = integerdata (varargin)

  if (nargin < 3)
    error ("gallery: At least 3 arguments required for integerdata matrix");
  endif

  if (isnumeric (varargin{end}))
    jidx = varargin{end};
    svec = [varargin{:}];
    varargin(end) = [];
  elseif (ischar (varargin{end}))
    if (nargin < 4)
      error (["gallery: CLASS argument requires 4 inputs " ...
              "for integerdata matrix."]);
    endif
    jidx = varargin{end-1};
    svec = [varargin{1:end-1}];
    varargin(end-1) = [];
  else
    error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
            "for integerdata matrix"]);
  endif

  if (! (isnumeric (jidx) && isscalar (jidx)
         && jidx == fix (jidx)
         && jidx >= 0 && jidx <= 0xFFFFFFFF))
    error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
            "for integerdata matrix"]);
  endif

  ## Save and restore random state.  Initialization done so that reproducible
  ## data is available from gallery depending on the jidx and size vector.
  randstate = rand ("state");
  unwind_protect
    rand ("state", svec);
    A = randi (varargin{:});
  unwind_protect_cleanup
    rand ("state", randstate);
  end_unwind_protect

endfunction

function A = invhess (x, y)

  ## INVHESS  Inverse of an upper Hessenberg matrix.
  ##   INVHESS(X, Y), where X is an N-vector and Y an N-1 vector,
  ##   is the matrix whose lower triangle agrees with that of
  ##   ONES(N,1)*X' and whose strict upper triangle agrees with
  ##   that of [1 Y]*ONES(1,N).
  ##   The matrix is nonsingular if X(1) ~= 0 and X(i+1) ~= Y(i)
  ##   for all i, and its inverse is an upper Hessenberg matrix.
  ##   If Y is omitted it defaults to -X(1:N-1).
  ##   Special case: if X is a scalar INVHESS(X) is the same as
  ##   INVHESS(1:X).
  ##
  ##   References:
  ##   F.N. Valvi and V.S. Geroyannis, Analytic inverses and
  ##       determinants for a class of matrices, IMA Journal of Numerical
  ##       Analysis, 7 (1987), pp. 123-128.
  ##   W.-L. Cao and W.J. Stewart, A note on inverses of Hessenberg-like
  ##       matrices, Linear Algebra and Appl., 76 (1986), pp. 233-240.
  ##   Y. Ikebe, On inverses of Hessenberg matrices, Linear Algebra and
  ##       Appl., 24 (1979), pp. 93-97.
  ##   P. Rozsa, On the inverse of band matrices, Integral Equations and
  ##       Operator Theory, 10 (1987), pp. 82-95.

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for invhess matrix");
  elseif (! isnumeric (x))
    error ("gallery: X must be numeric for invhess matrix");
  endif

  if (isscalar (x) && fix (x) == x)
    n = x;
    x = 1:n;
  elseif (! isscalar (x) && isvector (x))
    n = numel (x);
  else
    error ("gallery: X must be an integer scalar, or a vector for invhess matrix");
  endif

  if (nargin < 2)
    y = -x(1:end-1);
  elseif (! isvector (y) || numel (y) != numel (x) -1)
    error ("gallery: Y must be a vector of length -1 than X for invhess matrix");
  endif

  x = x(:);
  y = y(:);

  ##  FIXME: On next line, z = x'; A = z(ones(n,1),:) would be more efficient.
  A = ones (n, 1) * x';
  for j = 2:n
    A(1:j-1,j) = y(1:j-1);
  endfor

endfunction

function A = invol (n)

  ## INVOL   An involutory matrix.
  ##   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.

  if (nargin < 1)
    error ("gallery: 1 argument is required for invol matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for invol matrix");
  endif

  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,:);
  endfor

endfunction

function [A, detA] = ipjfact (n, k = 0)

  ## IPJFACT   A Hankel matrix with factorial elements.
  ##   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.
  ##   The determinant and inverse are known explicitly.
  ##   If a second output argument is present, d = DET(A) is returned:
  ##   [A, d] = IPJFACT(N, K);
  ##
  ##   Suggested by P. R. Graves-Morris.
  ##
  ##   Reference:
  ##   M.J.C. Gover, The explicit inverse of factorial Hankel matrices,
  ##   Dept. of Mathematics, University of Bradford, 1993.

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for ipjfact matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for ipjfact matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for ipjfact matrix");
  endif

  c = cumprod (2:n+1);
  d = cumprod (n+1:2*n) * c(n-1);

  A = hankel (c, d);

  switch (k)
    case 0, # do nothing
    case 1, A = ones (n) ./ A;
    otherwise
      error ("gallery: K must have a value of 0 or 1 for ipjfact matrix");
  endswitch

  if (nargout == 2)
    d = 1;

    if (k == 0)
      for i = 1:n-1
        d *= prod (1:i+1) * prod (1:n-i);
      endfor
      d *= prod (1:n+1);

    elseif (k == 1)
      for i = 0:n-1
        d *= prod (1:i) / prod (1:n+1+i);
      endfor
      if (rem (n*(n-1)/2, 2))
        d = -d;
      endif

    else
      error ("gallery: K must have a value of 0 or 1 for ipjfact matrix");
    endif

    detA = d;
  endif

endfunction

function J = jordbloc (n, lambda = 1)

  ## JORDBLOC  Jordan block.
  ##   JORDBLOC(N, LAMBDA) is the N-by-N Jordan block with eigenvalue
  ##   LAMBDA.  LAMBDA = 1 is the default.

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for jordbloc matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for jordbloc matrix");
  elseif (! isnumeric (lambda) || ! isscalar (lambda))
    error ("gallery: LAMBDA must be a numeric scalar for jordbloc matrix");
  endif

  J = lambda * eye (n) + diag (ones (n-1, 1), 1);

endfunction

function U = kahan (n, theta = 1.2, pert = 25)

  ## KAHAN  Kahan matrix - upper trapezoidal.
  ##   KAHAN(N, THETA) is an upper trapezoidal matrix
  ##   that 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).
  ##   The parameter THETA defaults to 1.2.
  ##   The useful range of THETA is 0 < THETA < PI.
  ##
  ##   To ensure that the QR factorization with column pivoting does not
  ##   interchange columns in the presence of rounding errors, the diagonal
  ##   is perturbed by PERT*EPS*diag( [N:-1:1] ).
  ##   The default is PERT = 25, which ensures no interchanges for KAHAN(N)
  ##   up to at least N = 90 in IEEE arithmetic.
  ##   KAHAN(N, THETA, PERT) uses the given value of PERT.
  ##
  ##   The inverse of KAHAN(N, THETA) is known explicitly: see
  ##   Higham (1987, p. 588), for example.
  ##   The diagonal perturbation was suggested by Christian Bischof.
  ##
  ##   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 < 1 || nargin > 3)
    error ("gallery: 1 to 3 arguments are required for kahan matrix");
  elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n)
    error ("gallery: N must be a 1 or 2 element integer for kahan matrix");
  elseif (! isnumeric (theta) || ! isscalar (theta))
    error ("gallery: THETA must be a numeric scalar for kahan matrix");
  elseif (! isnumeric (pert) || ! isscalar (pert))
    error ("gallery: PERT must be a numeric scalar for kahan matrix");
  endif

  ## Parameter n specifies dimension: r-by-n
  r = n(1);
  n = n(end);

  s = sin (theta);
  c = cos (theta);

  U = eye (n) - c * triu (ones (n), 1);
  U = diag (s.^[0:n-1]) * U + pert*eps* diag ([n:-1:1]);
  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
  endif

endfunction

function A = kms (n, rho = 0.5)

  ## KMS   Kac-Murdock-Szego Toeplitz matrix.
  ##   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 < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for lauchli matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for lauchli matrix");
  elseif (! isscalar (rho))
    error ("gallery: RHO must be a scalar for lauchli matrix");
  endif

  A = (1:n)'*ones (1,n);
  A = abs (A - A');
  A = rho .^ A;
  if (imag (rho))
    A = conj (tril (A,-1)) + triu (A);
  endif

endfunction

function B = krylov (A, x, j)

  ## KRYLOV    Krylov matrix.
  ##   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(RANDN(n)).
  ##
  ##   Reference:
  ##   G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  ##   Johns Hopkins University Press, Baltimore, Maryland, 1989, p. 369.

  if (nargin < 1 || nargin > 3)
    error ("gallery: 1 to 3 arguments are required for krylov matrix");
  elseif (! isnumeric (A) || ! issquare (A) || ndims (A) != 2)
    error ("gallery: A must be a square 2-D matrix for krylov matrix");
  endif

  n = length (A);
  if (isscalar (A))
    n = A;
    A = randn (n);
  endif

  if (nargin < 2)
    x = ones (n, 1);
  elseif (! isvector (x) || numel (x) != n)
    error ("gallery: X must be a vector of length equal to A for krylov matrix");
  endif

  if (nargin < 3)
    j = n;
  elseif (! isnumeric (j) || ! isscalar (j) || fix (j) != j)
    error ("gallery: J must be an integer for krylov matrix");
  endif

  B = ones (n, j);
  B(:,1) = x(:);
  for i = 2:j
    B(:,i) = A*B(:,i-1);
  endfor

endfunction

function A = lauchli (n, mu = sqrt (eps))

  ## LAUCHLI   Lauchli matrix - rectangular.
  ##   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 < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for lauchli matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for lauchli matrix");
  elseif (! isscalar (mu))
    error ("gallery: MU must be a scalar for lauchli matrix");
  endif

  A = [ones(1, n)
       mu*eye(n) ];

endfunction

function A = lehmer (n)

  ## LEHMER  Lehmer matrix - symmetric positive definite.
  ##   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.

  if (nargin < 1)
    error ("gallery: 1 argument is required for lehmer matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for lehmer matrix");
  endif

  A = ones (n, 1) * (1:n);
  A = A./A';
  A = tril (A) + tril (A, -1)';

endfunction

function T = lesp (n)

  ## LESP   A tridiagonal matrix with real, sensitive eigenvalues.
  ##   LESP(N) is an N-by-N matrix whose eigenvalues are real and smoothly
  ##   distributed in the interval approximately [-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,
  ##        Math. Comp., 57 (1991), pp. 221-237.
  ##   L.N. Trefethen, Pseudospectra of matrices, in Numerical Analysis 1991,
  ##        Proceedings of the 14th Dundee Conference,
  ##        D.F. Griffiths and G.A. Watson, eds, Pitman Research Notes in
  ##        Mathematics, volume 260, Longman Scientific and Technical, Essex,
  ##        UK, 1992, pp. 234-266.

  if (nargin < 1)
    error ("gallery: 1 argument is required for lesp matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for lesp matrix");
  endif

  x = 2:n;
  T = full (tridiag (ones (size (x)) ./x, -(2*[x n+1]+1), x));

endfunction

function A = lotkin (n)

  ## LOTKIN  Lotkin matrix.
  ##   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.

  if (nargin < 1)
    error ("gallery: 1 argument is required for lotkin matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for lotkin matrix");
  endif

  A = hilb (n);
  A(1,:) = ones (1, n);

endfunction

function A = minij (n)

  ## MINIJ   Symmetric positive definite matrix MIN(i,j).
  ##   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(N) (Givens' matrix) has tridiagonal inverse and
  ##               eigenvalues .5*sec^2([2r-1)PI/4N], r=1:N.
  ##   (N+1)*ONES(N)-A 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.)

  if (nargin < 1)
    error ("gallery: 1 argument is required for minij matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for minij matrix");
  endif

  A = bsxfun (@min, 1:n, (1:n)');

endfunction

function A = moler (n, alpha = -1)

  ## MOLER   Moler matrix - symmetric positive definite.
  ##   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 (1990) attributes the ALPHA = -1 matrix to Moler.
  ##
  ##   Reference:
  ##   J.C. Nash, Compact Numerical Methods for Computers: Linear
  ##   Algebra and Function Minimisation, second edition, Adam Hilger,
  ##   Bristol, 1990 (Appendix 1).

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for moler matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for moler matrix");
  elseif (! isscalar (alpha))
    error ("gallery: ALPHA must be a scalar for moler matrix");
  endif

  A = triw (n, alpha)' * triw (n, alpha);

endfunction

function [A, T] = neumann (n)

  ## NEUMANN  Singular matrix from the discrete Neumann problem (sparse).
  ##   NEUMANN(N) is the singular, row diagonally dominant matrix resulting
  ##   from discretizing the Neumann problem with the usual five point
  ##   operator on a regular mesh.
  ##   It has a one-dimensional null space with null vector ONES(N,1).
  ##   The dimension N should be a perfect square, or else a 2-vector,
  ##   in which case the dimension of the matrix is N(1)*N(2).
  ##
  ##   Reference:
  ##   R.J. Plemmons, Regular splittings and the discrete Neumann
  ##   problem, Numer. Math., 25 (1976), pp. 153-161.

  if (nargin < 1)
    error ("gallery: 1 argument is required for neumann matrix");
  elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n)
    error ("gallery: N must be a 1 or 2 element integer for neumann matrix");
  endif

  if (isscalar (n))
    m = sqrt (n);
    if (m^2 != n)
      error ("gallery: N must be a perfect square for neumann matrix");
    endif
    n(1) = m;
    n(2) = m;
  endif

  T = tridiag (n(1), -1, 2, -1);
  T(1,2) = -2;
  T(n(1),n(1)-1) = -2;

  A = kron (T, eye (n(2))) + kron (eye (n(2)), T);

endfunction

function A = normaldata (varargin)

  if (nargin < 2)
    error ("gallery: At least 2 arguments required for normaldata matrix");
  endif
  if (isnumeric (varargin{end}))
    jidx = varargin{end};
    svec = [varargin{:}];
    varargin(end) = [];
  elseif (ischar (varargin{end}))
    if (nargin < 3)
      error (["gallery: CLASS argument requires 3 inputs " ...
              "for normaldata matrix."]);
    endif
    jidx = varargin{end-1};
    svec = [varargin{1:end-1}];
    varargin(end-1) = [];
  else
    error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
            "for normaldata matrix"]);
  endif

  if (! (isnumeric (jidx) && isscalar (jidx)
         && jidx == fix (jidx)
         && jidx >= 0 && jidx <= 0xFFFFFFFF))
    error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
            "for normaldata matrix"]);
  endif

  ## Save and restore random state.  Initialization done so that reproducible
  ## data is available from gallery depending on the jidx and size vector.
  randstate = randn ("state");
  unwind_protect
    randn ("state", svec);
    A = randn (varargin{:});
  unwind_protect_cleanup
    randn ("state", randstate);
  end_unwind_protect

endfunction

function Q = orthog (n, k = 1)

  ## 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.
  ##           This is essentially the same matrix as FFT(EYE(N))/SQRT(N)!
  ##   K = 4:  Helmert matrix: a permutation of a lower Hessenberg matrix,
  ##           whose first row is ONES(1:N)/SQRT(N).
  ##   K = 5:  Q(i,j) = SIN( 2*PI*(i-1)*(j-1)/n ) + COS( 2*PI*(i-1)*(j-1)/n ).
  ##           Symmetric matrix arising in the Hartley transform.
  ##   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.
  ##   D. Bini and P. Favati, On a matrix algebra related to the discrete
  ##        Hartley transform, SIAM J. Matrix Anal. Appl., 14 (1993),
  ##        pp. 500-507.

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for orthog matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for orthog matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for orthog matrix");
  endif

  switch (k)
    case (1)
      ## E'vectors second difference matrix
      m = (1:n)'*(1:n) * (pi/(n+1));
      Q = sin (m) * sqrt (2/(n+1));

    case (2)
      m = (1:n)'*(1:n) * (2*pi/(2*n+1));
      Q = sin (m) * (2/ sqrt (2*n+1));

    case (3)
      ## Vandermonde based on roots of unity
      m = 0:n-1;
      Q = exp (m'*m*2*pi* sqrt (-1) / n) / sqrt (n);

    case (4)
      ## Helmert matrix
      Q = tril (ones (n));
      Q(1,2:n) = ones (1, n-1);
      for i = 2:n
        Q(i,i) = -(i-1);
      endfor
      Q = diag (sqrt ([n 1:n-1] .* [1:n])) \ Q;

    case (5)
      ## Hartley matrix
      m = (0:n-1)'*(0:n-1) * (2*pi/n);
      Q = (cos (m) + sin (m)) / sqrt (n);

    case (-1)
      ##  extrema of T(n-1)
      m = (0:n-1)'*(0:n-1) * (pi/(n-1));
      Q = cos (m);

    case (-2)
      ## zeros of T(n)
      m = (0:n-1)'*(.5:n-.5) * (pi/n);
      Q = cos (m);

    otherwise
      error ("gallery: unknown K '%d' for orthog matrix", k);
  endswitch

endfunction

function A = parter (n)

  ## PARTER    Parter matrix - a Toeplitz matrix with singular values near PI.
  ##   PARTER(N) is the matrix with (i,j) element 1/(i-j+0.5).
  ##   It is a Cauchy matrix and a Toeplitz matrix.
  ##
  ##   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.

  if (nargin < 1)
    error ("gallery: 1 argument is required for parter matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for parter matrix");
  endif

  A = cauchy ((1:n) + 0.5, -(1:n));

endfunction

function P = pei (n, alpha = 1)

  ## PEI    Pei matrix.
  ##   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 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for pei matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for pei matrix");
  elseif (! isnumeric (alpha) || ! isscalar (alpha))
    error ("gallery: ALPHA must be a scalar for pei matrix");
  endif

  P = alpha * eye (n) + ones (n);

endfunction

function A = poisson (n)

  ## POISSON   Block tridiagonal matrix from Poisson's equation (sparse).
  ##   POISSON(N) is the block tridiagonal matrix of order N^2
  ##   resulting from discretizing Poisson's equation with the
  ##   5-point operator on an N-by-N mesh.
  ##
  ##   Reference:
  ##   G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  ##   Johns Hopkins University Press, Baltimore, Maryland, 1989
  ##   (Section 4.5.4).

  if (nargin < 1)
    error ("gallery: 1 argument is required for poisson matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for poisson matrix");
  endif

  S = tridiag (n, -1, 2, -1);
  I = speye (n);
  A = kron (I, S) + kron (S, I);

endfunction

function A = prolate (n, w = 0.25)

  ## PROLATE   Prolate matrix - symmetric, ill-conditioned Toeplitz matrix.
  ##   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:
  ##   J.M. Varah. The Prolate matrix. Linear Algebra and Appl.,
  ##   187:269--278, 1993.

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for prolate matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for prolate matrix");
  elseif (! isnumeric (w) || ! isscalar (w))
    error ("gallery: W must be a scalar for prolate matrix");
  endif

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

endfunction

function H = randhess (x)

  ## NOTE: this function was named ohess in the original Test Matrix Toolbox
  ## RANDHESS  Random, orthogonal upper Hessenberg matrix.
  ##   H = RANDHESS(N) is an N-by-N real, random, orthogonal
  ##   upper Hessenberg matrix.
  ##   Alternatively, H = RANDHESS(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 (nargin < 1)
    error ("gallery: 1 argument is required for randhess matrix");
  elseif (! isnumeric (x) || ! isreal (x))
    error ("gallery: N or X must be numeric real values for randhess matrix");
  endif

  if (isscalar (x))
    n = x;
    x = rand (n-1, 1) * 2*pi;
    H = eye (n);
    H(n,n) = sign (randn);
  elseif (isvector (x))
    n = numel (x);
    H = eye (n);
    H(n,n) = sign (x(n)) + (x(n) == 0); # Second term ensures H(n,n) nonzero.
  else
    error ("gallery: N or X must be a scalar or a vector for randhess matrix");
  endif

  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,:) ];
  endfor

endfunction

function A = rando (n, k = 1)

  ## RANDO   Random matrix with elements -1, 0 or 1.
  ##   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 < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for rando matrix");
  elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n)
    error ("gallery: N must be an integer for rando matrix");
  elseif (! isnumeric (k) || ! isscalar (k))
    error ("gallery: K must be a numeric scalar for smoke matrix");
  endif

  ## Parameter n specifies dimension: m-by-n.
  m = n(1);
  n = n(end);

  switch (k)
    case (1), A =   floor (  rand(m, n) + 0.5);     # {0, 1}
    case (2), A = 2*floor (  rand(m, n) + 0.5) -1;  # {-1, 1}
    case (3), A =   round (3*rand(m, n) - 1.5);     # {-1, 0, 1}
    otherwise
      error ("gallery: unknown K '%d' for smoke matrix", k);
  endswitch

endfunction

function A = randsvd (n, kappa = sqrt (1/eps), mode = 3, kl = max (n) -1,
                      ku = kl)

  ## RANDSVD  Random matrix 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.
  ##
  ##   Reference:
  ##   N.J. Higham, Accuracy and Stability of Numerical Algorithms,
  ##      Society for Industrial and Applied Mathematics, Philadelphia, PA,
  ##      USA, 1996; sec. 26.3.
  ##
  ##   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 < 1 || nargin > 5)
    error ("gallery: 1 to 5 arguments are required for randsvd matrix");
  elseif (! isnumeric (n) || all (numel (n) != [1 2]) || fix (n) != n)
    error ("gallery: N must be a 1 or 2 element integer vector for randsvd matrix");
  elseif (! isnumeric (kappa) || ! isscalar (kappa))
    error ("gallery: KAPPA must be a numeric scalar for randsvd matrix");
  elseif (abs (kappa) < 1)
    error ("gallery: KAPPA must larger than or equal to 1 for randsvd matrix");
  elseif (! isnumeric (mode) || ! isscalar (mode))
    error ("gallery: MODE must be a numeric scalar for randsvd matrix");
  elseif (! isnumeric (kl) || ! isscalar (kl))
    error ("gallery: KL must be a numeric scalar for randsvd matrix");
  elseif (! isnumeric (ku) || ! isscalar (ku))
    error ("gallery: KU must be a numeric scalar for randsvd matrix");
  endif

  posdef = 0;
  if (kappa < 0)
    posdef = 1;
    kappa  = -kappa;
  endif

  ## Parameter n specifies dimension: m-by-n.
  m = n(1);
  n = n(end);
  p = min ([m n]);

  ## If A will be a vector
  if (p == 1)
    A = randn (m, n);
    A /= norm (A);
    return;
  endif

  ##  Set up vector sigma of singular values.
  switch (abs (mode))
    case (1)
      sigma = ones (p, 1) ./ kappa;
      sigma(1) = 1;
    case (2)
      sigma = ones (p, 1);
      sigma(p) = 1 / kappa;
    case (3)
      factor = kappa^(-1/(p-1));
      sigma  = factor.^[0:p-1];
    case (4)
      sigma = ones (p, 1) - (0:p-1)'/(p-1)*(1-1/kappa);
    case (5)
      ## In this case cond (A) <= kappa.
      rand ("uniform");
      sigma = exp (-rand (p, 1) * log (kappa));
    otherwise
      error ("gallery: unknown MODE '%d' for randsvd matrix", mode);
  endswitch

  ##  Convert to diagonal matrix of singular values.
  if (mode < 0)
    sigma = sigma(p:-1:1);
  endif
  sigma = diag (sigma);

  if (posdef)
    ## handle case where KAPPA was negative
    Q = qmult (p);
    A = Q' * sigma * Q;
    A = (A + A') / 2;  # Ensure matrix is symmetric.
    return;
  endif

  if (m != n)
    ## Expand to m-by-n diagonal matrix
    sigma(m, n) = 0;
  endif

  if (kl == 0 && ku == 0)
    ## Diagonal matrix requested - nothing more to do.
    A = sigma;
  else
    ##  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);
    endif
  endif

endfunction

function A = redheff (n)

  ## REDHEFF    A (0,1) matrix of Redheffer associated with the Riemann hypothesis.
  ##   A = REDHEFF(N) is an N-by-N matrix of 0s and 1s defined by
  ##       A(i,j) = 1 if j = 1 or if i divides j,
  ##       A(i,j) = 0 otherwise.
  ##   It has N - FLOOR(LOG2(N)) - 1 eigenvalues equal to 1,
  ##   a real eigenvalue (the spectral radius) approximately SQRT(N),
  ##   a negative eigenvalue approximately -SQRT(N),
  ##   and the remaining eigenvalues are provably ``small''.
  ##   Barrett and Jarvis (1992) conjecture that
  ##     ``the small eigenvalues all lie inside the unit circle
  ##       ABS(Z) = 1'',
  ##   and a proof of this conjecture, together with a proof that some
  ##   eigenvalue tends to zero as N tends to infinity, would yield
  ##   a new proof of the prime number theorem.
  ##   The Riemann hypothesis is true if and only if
  ##   DET(A) = O( N^(1/2+epsilon) ) for every epsilon > 0
  ##                                     ('!' denotes factorial).
  ##   See also RIEMANN.
  ##
  ##   Reference:
  ##   W.W. Barrett and T.J. Jarvis,
  ##   Spectral Properties of a Matrix of Redheffer,
  ##   Linear Algebra and Appl., 162 (1992), pp. 673-683.

  if (nargin < 1)
    error ("gallery: 1 argument is required for redheff matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for redheff matrix");
  endif

  i = (1:n)' * ones (1, n);
  A = ! rem (i', i);
  A(:,1) = ones (n, 1);

endfunction

function A = riemann (n)

  ## RIEMANN    A matrix associated with the Riemann hypothesis.
  ##   A = RIEMANN(N) is an N-by-N matrix for which the
  ##   Riemann hypothesis 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.
  ##
  ##   See also REDHEFF.
  ##
  ##   Reference:
  ##   F. Roesler, Riemann's hypothesis as an eigenvalue problem,
  ##   Linear Algebra and Appl., 81 (1986), pp. 153-198.

  if (nargin < 1)
    error ("gallery: 1 argument is required for riemann matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for riemann matrix");
  endif

  n += 1;
  i = (2:n)' * ones (1, n-1);
  j = i';
  A = i .* (! rem (j, i)) - ones (n-1);

endfunction

function A = ris (n)

  ## NOTE: this function was named dingdong in the original Test Matrix Toolbox
  ## RIS       Dingdong matrix - a symmetric Hankel matrix.
  ##   A = RIS(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, second edition, Adam Hilger,
  ##   Bristol, 1990 (Appendix 1).

  if (nargin < 1)
    error ("gallery: 1 argument is required for ris matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for ris matrix");
  endif

  p = -2*(1:n) + (n+1.5);
  A = cauchy (p);

endfunction

function A = smoke (n, k = 0)

  ## SMOKE     Smoke matrix - complex, with a 'smoke ring' pseudospectrum.
  ##   SMOKE(N) is an N-by-N matrix with 1s on the
  ##   superdiagonal, 1 in the (N,1) position, and powers of
  ##   roots of unity along the diagonal.
  ##   SMOKE(N, 1) is the same except for a zero (N,1) element.
  ##   The eigenvalues of SMOKE(N, 1) are the N'th roots of unity;
  ##   those of SMOKE(N) are the N'th roots of unity times 2^(1/N).
  ##
  ##   Try PS(SMOKE(32)).  For SMOKE(N, 1) the pseudospectrum looks
  ##   like a sausage folded back on itself.
  ##   GERSH(SMOKE(N, 1)) is interesting.
  ##
  ##   Reference:
  ##   L. Reichel and L.N. Trefethen, Eigenvalues and pseudo-eigenvalues of
  ##   Toeplitz matrices, Linear Algebra and Appl., 162-164:153-185, 1992.

  if (nargin < 1 || nargin > 2)
    error ("gallery: 1 to 2 arguments are required for smoke matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be an integer for smoke matrix");
  elseif (! isnumeric (n) || ! isscalar (n))
    error ("gallery: K must be a numeric scalar for smoke matrix");
  endif

  w = exp (2*pi*i/n);
  A = diag ( [w.^(1:n-1) 1] ) + diag (ones (n-1,1), 1);

  switch (k)
    case 0, A(n,1) = 1;
    case 1, # do nothing
    otherwise,
      error ("gallery: K must have a value of 0 or 1 for smoke matrix");
  endswitch

endfunction

function T = toeppd (n, m = n, w = rand (m,1), theta = rand (m,1))

  ## NOTE: this function was named pdtoep in the original Test Matrix Toolbox
  ## TOEPPD   Symmetric positive definite Toeplitz matrix.
  ##   TOEPPD(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).
  ##
  ##   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 < 1 || nargin > 4)
    error ("gallery: 1 to 4 arguments are required for toeppd matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be a numeric integer for toeppd matrix");
  elseif (! isnumeric (m) || ! isscalar (m) || fix (m) != m)
    error ("gallery: M must be a numeric integer for toeppd matrix");
  elseif (numel (w) != m || numel (theta) != m)
    error ("gallery: W and THETA must be vectors of length M for toeppd matrix");
  endif

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

  for i = 1:m
    T += w(i) * cos (theta(i)*E);
  endfor

endfunction

function P = toeppen (n, a = 1, b = -10, c = 0, d = 10, e = 1)

  ## NOTE: this function was named pentoep in the original Test Matrix Toolbox
  ## TOEPPEN   Pentadiagonal Toeplitz matrix (sparse).
  ##   P = TOEPPEN(N, A, B, C, D, E) is the N-by-N pentadiagonal
  ##   Toeplitz matrix with diagonals composed of the numbers
  ##   A =: P(3,1), B =: P(2,1), C =: P(1,1), D =: P(1,2), E =: P(1,3).
  ##   Default: (A,B,C,D,E) = (1,-10,0,10,1) (a matrix of Rutishauser).
  ##             This matrix has eigenvalues lying approximately on
  ##             the line segment 2*cos(2*t) + 20*i*sin(t).
  ##
  ##   Interesting plots are
  ##   PS(FULL(TOEPPEN(32,0,1,0,0,1/4)))  - 'triangle'
  ##   PS(FULL(TOEPPEN(32,0,1/2,0,0,1)))  - 'propeller'
  ##   PS(FULL(TOEPPEN(32,0,1/2,1,1,1)))  - 'fish'
  ##
  ##   References:
  ##   R.M. Beam and R.F. Warming, The asymptotic spectra of
  ##      banded Toeplitz and quasi-Toeplitz matrices, SIAM J. Sci.
  ##      Comput. 14 (4), 1993, pp. 971-1006.
  ##   H. Rutishauser, On test matrices, Programmation en Mathematiques
  ##      Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165,
  ##      1966, pp. 349-365.

  if (nargin < 1 || nargin > 6)
    error ("gallery: 1 to 6 arguments are required for toeppen matrix");
  elseif (! isnumeric (n) || ! isscalar (n) || fix (n) != n)
    error ("gallery: N must be a numeric integer for toeppen matrix");
  elseif (any (! cellfun ("isnumeric", {a b c d e}))
          || any (cellfun ("numel", {a b c d e}) != 1))
    error ("gallery: A, B, C, D and E must be numeric scalars for toeppen matrix");
  endif

  P = spdiags ([a*ones(n,1) b*ones(n,1) c*ones(n,1) d*ones(n,1) e*ones(n,1)],
                -2:2, n, n);

endfunction

function T = tridiag (n, x = -1, y = 2, z = -1)

  ## TRIDIAG  Tridiagonal matrix (sparse).
  ##   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 1977)
  ##            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).
  ##
  ##   References:
  ##   J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  ##     Birkhauser, Basel, and Academic Press, New York, 1977, p. 155.
  ##   D.E. Rutherford, Some continuant determinants arising in physics and
  ##     chemistry---II, Proc. Royal Soc. Edin., 63, A (1952), pp. 232-241.

  if (nargin != 1 && nargin != 3 && nargin != 4)
    error ("gallery: 1, 3, or 4 arguments are required for tridiag matrix");
  elseif (nargin == 3)
    z = y;
    y = x;
    x = n;
  endif

  ## Force column vectors
  x = x(:);
  y = y(:);
  z = z(:);

  if (isscalar (x) && isscalar (y) && isscalar (z))
    x *= ones (n-1, 1);
    z *= ones (n-1, 1);
    y *= ones (n,   1);
  elseif (numel (y) != numel (x) + 1)
    error ("gallery: X must have one element less than Y for tridiag matrix");
  elseif (numel (y) != numel (z) + 1)
    error ("gallery: Z must have one element less than Y for tridiag matrix");
  endif

  ##  T = diag (x, -1) + diag (y) + diag (z, 1);  # For non-sparse matrix.
  n = numel (y);
  T = spdiags ([[x;0] y [0;z]], -1:1, n, n);

endfunction

function t = triw (n, alpha = -1, k = n(end) - 1)

  ## TRIW   Upper triangular matrix discussed by Wilkinson and others.
  ##   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 Kahan, Golub and Wilkinson.
  ##
  ##   Ostrowski (1954) shows that
  ##     COND(TRIW(N,2)) = COT(PI/(4*N))^2,
  ##   and for large ABS(ALPHA),
  ##     COND(TRIW(N,ALPHA)) is approximately ABS(ALPHA)^N*SIN(PI/(4*N-2)).
  ##
  ##   Adding -2^(2-N) to the (N,1) element makes TRIW(N) singular,
  ##   as does adding -2^(1-N) to all elements in the first column.
  ##
  ##   References:
  ##   G.H. Golub and J.H. Wilkinson, Ill-conditioned eigensystems and the
  ##      computation of the Jordan canonical form, SIAM Review,
  ##      18(4), 1976, pp. 578-619.
  ##   W. Kahan, Numerical linear algebra, Canadian Math. Bulletin,
  ##      9 (1966), pp. 757-801.
  ##   A.M. Ostrowski, On the spectrum of a one-parametric family of
  ##      matrices, J. Reine Angew. Math., 193 (3/4), 1954, pp. 143-160.
  ##   J.H. Wilkinson, Singular-value decomposition---basic aspects,
  ##      in D.A.H. Jacobs, ed., Numerical Software---Needs and Availability,
  ##      Academic Press, London, 1978, pp. 109-135.

  if (nargin < 1 || nargin > 3)
    error ("gallery: 1 to 3 arguments are required for triw matrix");
  elseif (! isnumeric (n) || all (numel (n) != [1 2]))
    error ("gallery: N must be a 1 or 2 elements vector for triw matrix");
  elseif (! isscalar (alpha))
    error ("gallery: ALPHA must be a scalar for triw matrix");
  elseif (! isscalar (k) || ! isnumeric (k) || fix (k) != k || k < 0)
    error ("gallery: K must be a numeric integer >= 0 for triw matrix");
  endif

  m = n(1);              # Parameter n specifies dimension: m-by-n.
  n = n(end);

  t = tril (eye (m, n) + alpha * triu (ones (m, n), 1), k);

endfunction

function A = uniformdata (varargin)

  if (nargin < 2)
    error ("gallery: At least 2 arguments required for uniformdata matrix");
  endif
  if (isnumeric (varargin{end}))
    jidx = varargin{end};
    svec = [varargin{:}];
    varargin(end) = [];
  elseif (ischar (varargin{end}))
    if (nargin < 3)
      error (["gallery: CLASS argument requires 3 inputs " ...
              "for uniformdata matrix."]);
    endif
    jidx = varargin{end-1};
    svec = [varargin{1:end-1}];
    varargin(end-1) = [];
  else
    error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
            "for uniformdata matrix"]);
  endif

  if (! (isnumeric (jidx) && isscalar (jidx)
         && jidx == fix (jidx)
         && jidx >= 0 && jidx <= 0xFFFFFFFF))
    error (["gallery: J must be an integer in the range [0, 2^32-1] " ...
            "for uniformdata matrix"]);
  endif

  ## Save and restore random state.  Initialization done so that reproducible
  ## data is available from gallery depending on the jidx and size vector.
  randstate = rand ("state");
  unwind_protect
    rand ("state", svec);
    A = rand (varargin{:});
  unwind_protect_cleanup
    rand ("state", randstate);
  end_unwind_protect

endfunction

function A = wathen (nx, ny, k = 0)

  ## WATHEN returns the Wathen matrix.
  ##
  ## Discussion:
  ##
  ##   The Wathen matrix is a finite element matrix which is sparse.
  ##
  ##   The entries of the matrix depend in part on a physical quantity
  ##   related to density.  That density is here assigned random values between
  ##   0 and 100.
  ##
  ##   A = WATHEN ( NX, NY ) is a sparse random N-by-N finite element matrix
  ##   where N = 3*NX*NY + 2*NX + 2*NY + 1.
  ##
  ##   A is the consistent mass matrix for a regular NX-by-NY
  ##   grid of 8-node (serendipity) elements in 2 space dimensions.
  ##
  ##   Here is an illustration for NX = 3, NX = 2:
  ##
  ##    23-24-25-26-27-28-29
  ##     |     |     |     |
  ##    19    20    21    22
  ##     |     |     |     |
  ##    12-13-14-15-16-17-18
  ##     |     |     |     |
  ##     8     9    10    11
  ##     |     |     |     |
  ##     1--2--3--4--5--6--7
  ##
  ##   For this example, the total number of nodes is, as expected,
  ##
  ##     N = 3 * 3 * 2 + 2 * 2 + 2 * 3 + 1 = 29.
  ##
  ##   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).
  ##
  ##   A = WATHEN ( NX, NY, 1 ) returns the diagonally scaled matrix.
  ##
  ## Modified:
  ##
  ##   17 September 2007
  ##
  ## Author:
  ##
  ##   Nicholas Higham
  ##
  ## Reference:
  ##
  ##   Nicholas Higham,
  ##   Algorithm 694: A Collection of Test Matrices in MATLAB,
  ##   ACM Transactions on Mathematical Software,
  ##   Volume 17, Number 3, September 1991, pages 289-305.
  ##
  ##   Andrew Wathen,
  ##   Realistic eigenvalue bounds for the Galerkin mass matrix,
  ##   IMA Journal of Numerical Analysis,
  ##   Volume 7, 1987, pages 449-457.
  ##
  ## Parameters:
  ##
  ##   Input, integer NX, NY, the number of elements in the X and Y directions
  ##   of the finite element grid.  NX and NY must each be at least 1.
  ##
  ##   Optional input, integer K, is used to request that the diagonally scaled
  ##   version of the matrix be returned.  This happens if K is specified with
  ##   the value 1.
  ##
  ##   Output, sparse real A(N,N), the matrix.  The dimension N is determined by
  ##   NX and NY, as described above.  A is stored in the MATLAB sparse matrix
  ##   format.

  if (nargin < 2 || nargin > 3)
    error ("gallery: 2 or 3 arguments are required for wathen matrix");
  elseif (! isnumeric (nx) || ! isscalar (nx) || nx < 1)
    error ("gallery: NX must be a positive scalar for wathen matrix");
  elseif (! isnumeric (ny) || ! isscalar (ny) || ny < 1)
    error ("gallery: NY must be a positive scalar for wathen matrix");
  elseif (! isscalar (k))
    error ("gallery: K must be a scalar for wathen matrix");
  endif

  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 = sparse (n, n);

  rho = 100 * rand (nx, ny);

  for j = 1:ny
    for i = 1:nx
      ##
      ##   For the element (I,J), determine the indices of the 8 nodes.
      ##
      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);
        endfor
      endfor

    endfor
  endfor

  ## If requested, return A with diagonal scaling.
  if (k)
    A = diag (diag (A)) \ A;
  endif

endfunction

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 (nargin < 1)
    error ("gallery: 1 argument is required for wilk matrix");
  elseif (! isnumeric (n) || ! isscalar (n))
    error ("gallery: N must be a numeric scalar for wilk matrix");
  endif

  if (n == 3)
    ## Wilkinson (1961) p.323.
    A = [ 1e-10   0.9  -0.4
          0       0.9  -0.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)
    ## Wilkinson (1965), p.308.
    E = diag (ones (n-1, 1), 1);
    m = (n-1)/2;
    A = diag (abs (-m:m)) + E + E';

  else
    error ("gallery: unknown N '%d' for wilk matrix", n);
  endif

endfunction

## NOTE: bandred is part of the Test Matrix Toolbox and is used by randsvd()
function A = bandred (A, kl, ku)

  ## BANDRED  Band reduction by two-sided unitary transformations.
  ##   B = BANDRED(A, KL, KU) is a matrix unitarily 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.
  ##
  ##   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!
  ##
  ##   Reference:
  ##   G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  ##   Johns Hopkins University Press, Baltimore, Maryland, 1989.
  ##   Section 5.4.3.

  ##  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 = false;
  if (ku == 0)
    flip = true;
    A = A';
    [ku, kl] = deal (kl, ku);
  endif

  [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);
    endif

    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);
    endif
  endfor

  if (flip)
    A = A';
  endif

endfunction

## NOTE: qmult is part of the Test Matrix Toolbox and is used by randsvd()
function B = qmult (A)

  ## QMULT  Pre-multiply by random orthogonal matrix.
  ##   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)).
  ##
  ##   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 (isscalar (A))
    n = A;
    A = eye (n);
  endif

  d = zeros (n);

  for k = n-1:-1:1
    ## Generate random Householder transformation.
    x = randn (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);
  endfor

  ## Tidy up signs
  for i = 1:n-1
    A(i,:) = d(i)*A(i,:);
  endfor
  A(n,:) = A(n,:) * sign (randn);
  B = A;

endfunction


## BIST testing for just a few functions to verify that the main gallery
## dispatch function works.
%!assert (gallery ("clement", 3), [0 1 0; 2 0 2; 0 1 0])
%!assert (gallery ("invhess", 2), [1 -1; 1 2])

## Test input validation of main dispatch function only
%!error <Invalid call> gallery ()
%!error <NAME must be a string> gallery (123)
%!error <matrix binomial not implemented> gallery ("binomial")
%!error <unknown matrix with NAME foobar> gallery ("foobar")

## BIST testing for individual gallery functions
%!assert (gallery ("minij", 4), [1 1 1 1; 1 2 2 2; 1 2 3 3; 1 2 3 4])
%!assert (gallery ("minij", 1), 1)
%!assert (gallery ("minij", 0), [])
%!assert (gallery ("minij", -1), [])

%!test
%! exp = 1 ./ [
%!   2  3  4  5  6
%!   3  4  5  6  7
%!   4  5  6  7  8
%!   5  6  7  8  9
%!   6  7  8  9  10];
%! assert (gallery ("cauchy", 5), exp);
%! assert (gallery ("cauchy", 1:5), exp);
%! assert (gallery ("cauchy", 1:5, 1:5), exp);
%!
%! exp = 1 ./ [
%!   1  2  3  4  5
%!   2  3  4  5  6
%!   3  4  5  6  7
%!   4  5  6  7  8
%!   5  6  7  8  9];
%! assert (gallery ("cauchy", 0:4, 1:5), exp);
%! assert (gallery ("cauchy", 1:5, 0:4), exp);
%! assert (gallery ("cauchy", 1:5, 4:-1:0), fliplr (exp));
%!
%! exp = 1 ./ [
%!  -1  0  1  2  3
%!   0  1  2  3  4
%!   1  2  3  4  5
%!   2  3  4  5  6
%!   3  4  5  6  7];
%! assert (gallery ("cauchy", 1:5, -2:2), exp);
%!
%! exp = 1 ./ [
%!    8  18  -4  2
%!   13  23   1  7
%!    9  19  -3  3
%!   15  25   3  9];
%! assert (gallery ("cauchy", [-2 3 -1 5], [10 20 -2 4]), exp);
%! assert (gallery ("cauchy", [-2 3 -1 5], [10 20 -2 4]'), exp);
%! assert (gallery ("cauchy", [-2 3 -1 5]', [10 20 -2 4]), exp);

%!assert (size (gallery ("chebspec", 5)), [5 5])
%!assert (size (gallery ("chebspec", 5, 1)), [5 5])
%!assert (size (gallery ("chebspec", 5, 0)), [5 5])

%!assert (size (gallery ("chebvand", 7)), [7 7])
%!assert (size (gallery ("chebvand", 1:7)), [7 7])
%!assert (size (gallery ("chebvand", 5, 7)), [5 7])

%!assert (size (gallery ("chow", 5)), [5 5])
%!assert (size (gallery ("chow", 5, 6)), [5 5])
%!assert (size (gallery ("chow", 5, 6, 7)), [5 5])

%!assert (gallery ("circul", 3), [1 2 3; 3 1 2; 2 3 1])
%!assert (gallery ("circul", [1 3 6]), [1 3 6; 6 1 3; 3 6 1])

%!assert (size (gallery ("clement", 5)), [5 5])
%!assert (size (gallery ("clement", 5, 1)), [5 5])
%!assert (size (gallery ("clement", 5, 0)), [5 5])

%!assert (size (gallery ("compar", ones (5))), [5 5])
%!assert (size (gallery ("compar", ones (5), 0)), [5 5])
%!assert (size (gallery ("compar", ones (5), 1)), [5 5])

%!assert (size (gallery ("condex", 4)), [4 4])
%!assert (size (gallery ("condex", 4, 1)), [4 4])
%!assert (size (gallery ("condex", 4, 1, 50)), [4 4])

%!assert (size (gallery ("cycol", [4 5])), [4 5])
%!assert (size (gallery ("cycol", [4 5], 1)), [4 5])
%!assert (size (gallery ("cycol", 4)), [4 4])
%!assert (size (gallery ("cycol", 4, 1)), [4 4])

%!assert (size (gallery ("dorr", 4)), [4 4])
%!assert (cellfun (@rows, nthargout (1:3, @gallery, "dorr", 4)), [3 4 3])

%!assert (size (gallery ("dramadah", 5)), [5 5])
%!assert (size (gallery ("dramadah", 5, 2)), [5 5])

%!test
%! exp = [
%!   0   1   2   3   4
%!   1   0   1   2   3
%!   2   1   0   1   2
%!   3   2   1   0   1
%!   4   3   2   1   0];
%! assert (gallery ("fiedler", 5), exp);
%! assert (gallery ("fiedler", 1:5), exp);
%! assert (gallery ("fiedler", -2:2), exp);
%! assert (gallery ("fiedler", 2:5), exp(1:4,1:4));

%!assert (size (gallery ("forsythe", 5)), [5 5])
%!assert (size (gallery ("forsythe", 5, 1, 0.5)), [5 5])
%!assert (size (gallery ("forsythe", 5, 4, 7)), [5 5])

%!assert (size (gallery ("frank", 5)), [5 5])
%!assert (size (gallery ("frank", 5, 1)), [5 5])

%!assert (size (gallery ("gcdmat", 5)), [5 5])

%!assert (size (gallery ("gearmat", 5)), [5 5])
%!assert (size (gallery ("gearmat", 5, 4)), [5 5])
%!assert (size (gallery ("gearmat", 5, 4, 3)), [5 5])

%!assert (size (gallery ("grcar", 5)), [5 5])
%!assert (size (gallery ("grcar", 5, 2)), [5 5])

%!error <N must be even> gallery ("hanowa", 5)
%!assert (size (gallery ("hanowa", 6, 5)), [6 6])
%!assert (size (gallery ("hanowa", 6, 5)), [6 6])

%!assert (size (gallery ("house", [1:5]')), [5 1])
%!assert (cellfun (@rows, nthargout (1:2, @gallery, "house", [1:5]')), [5 1])

%!assert (size (gallery ("integerdata", 5, [3 2], 5)), [3 2])
%!assert (size (gallery ("integerdata", 5, [3 2 6], 5)), [3 2 6])

%!assert (size (gallery ("invhess", 1:4, 1:3)), [4 4])

%!assert (size (gallery ("invol", 4)), [4 4])

%!assert (size (gallery ("ipjfact", 4)), [4 4])
%!assert (size (gallery ("ipjfact", 4, 0)), [4 4])
%!assert (size (gallery ("ipjfact", 4, 1)), [4 4])

%!assert (size (gallery ("jordbloc", 4)), [4 4])
%!assert (size (gallery ("jordbloc", 4, 1)), [4 4])
%!assert (size (gallery ("jordbloc", 4, 3)), [4 4])

%!assert (size (gallery ("kahan", 4)), [4 4])
%!assert (size (gallery ("kahan", [4 5])), [4 5])
%!assert (size (gallery ("kahan", [4 5], 1)), [4 5])
%!assert (size (gallery ("kahan", [4 5], 1, 30)), [4 5])

%!assert (size (gallery ("kms", 5)), [5 5])
%!assert (size (gallery ("kms", 5, 0.2)), [5 5])

%!assert (size (gallery ("krylov", 4)), [4 4])
%!assert (size (gallery ("krylov", ones (4))), [4 4])
%!assert (size (gallery ("krylov", ones (4), [.2 .3 .4 .5], 3)), [4 3])

%!assert (size (gallery ("lauchli", 5)), [6 5])
%!assert (size (gallery ("lauchli", 5, 3)), [6 5])

%!assert (size (gallery ("lehmer", 5)), [5 5])

%!assert (size (gallery ("lesp", 5)), [5 5])

%!assert (size (gallery ("lotkin", 5)), [5 5])

%!assert (size (gallery ("minij", 5)), [5 5])

%!assert (size (gallery ("moler", 5)), [5 5])
%!assert (size (gallery ("moler", 5, 0.2)), [5 5])

%!assert (size (gallery ("neumann", 4)), [4 4])

%!assert (size (gallery ("normaldata", [5 4 6], 3)), [5 4 6])

%!assert (size (gallery ("orthog", 5)), [5 5])
%!assert (size (gallery ("orthog", 5, 2)), [5 5])
%!assert (size (gallery ("orthog", 5, -2)), [5 5])

%!assert (size (gallery ("parter", 5)), [5 5])

%!assert (size (gallery ("pei", 5)), [5 5])
%!assert (size (gallery ("pei", 5, 4)), [5 5])

%!assert (size (gallery ("poisson", 1)), [1 1])
%!assert (size (gallery ("poisson", 4)), [16 16])
%!assert (size (gallery ("poisson", 5)), [25 25])

%!assert (size (gallery ("prolate", 5)), [5 5])
%!assert (size (gallery ("prolate", 5, 0.5)), [5 5])

%!assert (size (gallery ("randhess", 5)), [5 5])
%!assert (size (gallery ("randhess", 2:5)), [4 4])

%!assert (size (gallery ("rando", 5)), [5 5])
%!assert (size (gallery ("rando", 5, 2)), [5 5])

%!assert (size (gallery ("randsvd", 5)), [5 5])
%!assert (size (gallery ("randsvd", [5 3])), [5 3])

%!assert (size (gallery ("redheff", 5)), [5 5])
%!assert (size (gallery ("riemann", 5)), [5 5])
%!assert (size (gallery ("ris", 5)), [5 5])

%!assert (size (gallery ("smoke", 5)), [5 5])
%!assert (size (gallery ("smoke", 5, 1)), [5 5])
%!assert (gallery ("smoke", 5, 1)(5, 1), 0)

%!assert (size (gallery ("toeppd", 5)), [5 5])

%!assert (size (gallery ("toeppen", 5)), [5 5])

%!assert (size (gallery ("tridiag", 5)), [5 5])
%!assert (size (gallery ("tridiag", 1:4, 1:5, 1:4)), [5 5])
%!assert (gallery ("tridiag", 5), gallery ("tridiag", 5, -1, 2, -1))

%!assert (size (gallery ("triw", 5)), [5 5])

%!assert (size (gallery ("uniformdata", [5 3 4], 3)), [5 3 4])

%!assert (size (gallery ("wathen", 2, 3)), [29 29])

%!assert (cellfun (@rows, nthargout (1:2, @gallery, "wilk", 3)), [3 3])
%!assert (cellfun (@rows, nthargout (1:2, @gallery, "wilk", 4)), [4 4])
%!assert (size (gallery ("wilk", 5)), [5 5])
%!assert (size (gallery ("wilk", 21)), [21 21])