%%% ====================================================================
%%% BibTeX-file{
%%%     author          = "Nicholas J. Higham",
%%%     version         = "1.00",
%%%     date            = "October 24, 2005",
%%%     time            = "",
%%%     filename        = "u-manchester-mccm.bib",
%%%     address         = "Department of Mathematics
%%%                        University of Manchester
%%%                        Manchester M13 9PL
%%%                        England",
%%%     telephone       = "+44 (0)161 275 5800",
%%%     FAX             = "+44 (0)161 275 5819",
%%%     checksum        = ""
%%%     email           = "higham at ma.man.ac.uk (Internet)",
%%%     codetable       = "ISO/ASCII",
%%%     keywords        = "bibliography, Manchester, numerical
%%%                        analysis",
%%%     supported       = "yes",
%%%     docstring       = "BibTeX bibliography for the Numerical
%%%                        Analysis Reports of the Manchester Centre for
%%%                        Computational Mathematics
%%%                        (University of Manchester and UMIST)."
%%%     }
%%% ====================================================================

%%% The URL (Uniform Resource Locator) of this file is
%%%
%%%   URL = "http://www.maths.man.ac.uk/~nareports/u-manchester-mccm.bib"
%%%
%%% This is narep.bib, the index for the Numerical Analysis Reports of the
%%% Manchester Centre for Computational Mathematics
%%% (University of Manchester and UMIST).
%%% It is in BibTeX form, with additional fields.
%%% Entries are in reverse chronological order.

@String{inst-MCCM               = "Manchester Centre for
                                   Computational Mathematics"}
@String{inst-MCCM:adr           = "Manchester, England"}
%%% Numerical Analysis Reports: ISSN 1360-1725

@String{t-NAREP                 = "Numerical Analysis Report"}

@String{inst-U-MANCHESTER       = "University of Manchester"}
@String{inst-U-MANCHESTER:adr   = "Manchester M13 9PL, England"}

@String{pub-LONGMAN             = "Longman Scientific and Technical"}
@String{pub-LONGMAN:adr         = "Essex, UK"}

@String{pub-OUP                 = "Oxford University Press"}
@String{pub-OUP:adr             = "Oxford, UK"}

@String{pub-SV                  = "Spring{\-}er-Ver{\-}lag"}
@String{pub-SV:adr              = "Berlin, Germany~/ Heidelberg,
                                   Germany~/ London, UK~/ etc."}
%%% ====================================================================

@techreport{Guo:2005:SNM,
  author = "Chun-Hua Guo and Nicholas J. Higham",
  title = "A {Schur--Newton} Method for the Matrix $p$th Root and its Inverse",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 475",
  pages = "17",
  month = oct,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep475.pdf",
  abstract = "Newton's method for the inverse matrix $p$th root, $A^{-1/p}$, has
              the attraction that it involves only matrix multiplication. We
              show that if the starting matrix is $c^{-1}I$ for $c\in\R^+$ then
              the iteration converges quadratically to $A^{-1/p}$ if the
              eigenvalues of $A$ lie in a wedge-shaped convex set containing the
              disc $\{\, z: |z-c^p| < c^p\,\}$. We derive an optimal choice of
              $c$ for the case where $A$ has real, positive eigenvalues. An
              application is described to roots of transition matrices from
              Markov models, in which for certain problems the convergence
              condition is satisfied with $c=1$. Although the basic Newton
              iteration is numerically unstable, a coupled version is stable and
              a simple modification of it provides a new coupled iteration for
              the matrix $p$th root. For general matrices we develop a hybrid
              algorithm that computes a Schur decomposition, takes square roots
              of the upper (quasi)triangular factor, and applies the coupled
              Newton iteration to a matrix for which fast convergence is
              guaranteed. The new algorithm can be used to compute either
              $A^{1/p}$ or $A^{-1/p}$, and for large $p$ that are not highly
              composite it is more efficient than the method of Smith based
              entirely on the Schur decomposition.",
 keywords = "Matrix $p$th root, principal $p$th root, matrix logarithm, inverse,
             Newton's method, preprocessing, Schur decomposition, numerical
             stability, convergence, Markov model, transition matrix"
}

@techreport{Elman:05:IFI,
  author = "Howard Elman and Alison Ramage and David Silvester",
  title = "IFISS: a {Matlab} Toolbox for Modelling Incompressible Flow",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 474",
  month = oct,
  year = 2005,
  note = "Submitted to ACM Transactions on Mathematical Software",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep474.pdf",
  abstract = "IFISS is a graphical Matlab package for the interactive numerical study
            of incompressible flow problems. It includes algorithms for discretisation by
           mixed finite element methods and a posteriori error estimation of the
           computed solutions. The package can also be used as
           a computational laboratory for experimenting with state-of-the-art
          preconditioned iterative solvers for the discrete linear equation systems that
          arise in incompressible flow modelling. A unique feature of the package is its
          comprehensive nature; for each
          problem addressed, it enables the study of both discretisation and iterative
          solution algorithms as well as the interaction between the two and the
          resulting effect on overall efficiency."
}

@techreport{Tisseur:05:SCN,
  author = "Fran\c{c}oise Tisseur and Stef Graillat",
  title = "Structured Condition Numbers and Backward Errors
           in Scalar Product Spaces",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 473",
  pages = "21",
  month = sep,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep473.pdf",
  abstract = "We investigate the effect of structure-preserving perturbations on the
          solution to a linear system, matrix inversion, and distance to
          singularity. Particular attention is paid to linear and nonlinear
          structures that form Lie algebras, Jordan algebras and automorphism
          groups of a scalar product. These include complex symmetric,
          pseudo-symmetric, persymmetric, skew-symmetric, Hamiltonian, unitary,
          complex orthogonal and symplectic matrices. We show that under  mild
          assumptions on the scalar product, there is little or no difference
          between  structured and unstructured condition numbers and distance to
          singularity for matrices in Lie and Jordan algebras. Hence, for these
          classes of matrices, the usual unstructured perturbation analysis is
          sufficient. We show this is not true in general for structures in
          automorphism groups. Bounds and computable expressions for the
          structured condition numbers for a linear system and matrix inversion
          are derived for these nonlinear structures.
              Structured backward errors for the approximate solution of linear
          systems are also considered. Conditions are given for the structured
          backward error to be finite. We prove that for Lie and Jordan
          algebras, whenever the structured backward error is finite, it is
          within a small factor of the unstructured one. The same conclusion
          holds for orthogonal and unitary structures but cannot easily be
          extended to other matrix groups. For some nonlinear structures such as
          complex orthogonal, we show that the structured backward error can be
          rewritten in terms of a challenging Procrustes problem.
          This work extends and unifies earlier analyses.",
  keywords =  "structured matrices, normwise structured perturbations, structured
         linear systems, condition number, backward error, distance to
         singularity, Lie algebra, Jordan algebra, automorphism group, scalar
         product, bilinear form, sesquilinear form, orthosymmetric."
  }

@techreport{MCCM:2005:AR4,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 2004",
  type =         t-NAREP,
  number =       "No. 472",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "25",
  month =        sep,
  year =         "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep472.pdf"
  }

@techreport{Bhatia:2005:IDM,
  author = "Rajendra Bhatia",
  title = "Infinitely Divisible Matrices",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 471",
  pages = "36",
  month = aug,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep471.pdf"
}

@techreport{Lancaster:2005:NWS,
  author = "Peter Lancaster and Panayiotis Psarrakos",
  title = "A Note on Weak and Strong Linearizations of
           Regular Matrix Polynomials",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 470",
  pages = "6",
  month = jun,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep470.pdf"
}

@techreport{Berhanu:05:PBH,
  author = "Michael Berhanu",
  title = "Perturbation Bounds for Hyperbolic Matrix Factorizations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 469",
  pages = "24",
  month = jun,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep469.pdf",
  abstract = "Several matrix factorizations depend on orthogonal factors,
              matrices that preserve the Euclidean scalar product. Some of these
              factorizations can be extended and generalized to
              $(J,\Jt)$-orthogonal factors, that is, matrices that satisfy
              $H^TJH=\Jt$, where $J$ and $\Jt$ are diagonal with diagonal
              elements $\pm1$. The purpose of this work is to analyze the
              perturbation of matrix factorizations that have a
              $(J,\Jt)$-orthogonal or orthogonal factor and to give first order
              perturbation bounds. For each factorization analyzed, we give the
              sharpest possible first order bound,  which yields a condition
              number. The cost of computing these condition numbers is high. It
              is usually equivalent to computing the $2$-norm of an $n^2\times
              n^2$ matrix for a problem of size $n$. Thus, we also propose less
              sharp bounds that are less expensive to compute.",
  keywords = "Condition number, hyperbolic matrix factorization, manifold"
}

@techreport{Higham:05:IPL,
  author = "Nicholas J. Higham",
  title = "An Interview with {Peter Lancaster}",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 468",
  pages = "10",
  month = jun,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep468.pdf",
}

@techreport{Tisseur:05:SEC,
  author = "Michael Karow and Daniel Kressner and Fran\c{c}oise Tisseur",
  title = "Structured Eigenvalue Condition Numbers",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 467",
  pages = "16",
  month = apr,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep467.pdf",
  abstract = "This paper investigates the effect of structure-preserving
           perturbations on the eigenvalues of linearly and nonlinearly
           structured eigenvalue problems. Particular attention is paid to
           structures that form Jordan algebras, Lie algebras, and automorphism
           groups of a scalar product. Bounds and computable expressions for
           structured eigenvalue condition numbers are derived for these classes
           of matrices, which include complex symmetric, pseudo symmetric,
           persymmetric, skew-symmetric, Hamiltonian, symplectic, and orthogonal
           matrices. In particular we show that under mild assumptions on the
           scalar product, the structured and unstructured eigenvalue condition
           numbers are equal for structures in Jordan algebras. For Lie
           algebras, the effect on the condition number of incorporating
           structure varies greatly with the structure. We identify  Lie
           algebras for which structure does not affect the eigenvalue condition
           number.",
  keywords = "Structured eigenvalue problem, condition number, Jordan algebra,
           Lie algebra, automorphism group, symplectic, perplectic,
           pseudo-orthogonal, pseudo-unitary, complex symmetric, persymmetric,
           perskew-symmetric, Hamiltonian, skew-Hamiltonian, structure
           preservation."
}

@techreport{Mackey:2005:PPE,
  author = "D. Steven Mackey and Niloufer Mackey and Christian Mehl and
            Volker Mehrmann",
  title = "Palindromic Polynomial Eigenvalue Problems: {Good} Vibrations from
           Good Linearizations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type =         t-NAREP,
  number = "No. 466",
  pages = 26,
  month = apr,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep466.pdf",
  abstract = "Palindromic polynomial eigenvalue problems and related classes of
              structured eigenvalue problems are considered. These structures
              generalize the concepts of symplectic and Hamiltonian matrices to
              matrix polynomials. We discuss several applications where these
              matrix polynomials arise, and show how linearizations can be
              derived that re ect the structure of all these structured matrix
              polynomials and therefore preserve symmetries in the spectrum.",
   keywords = "nonlinear eigenvalue problem, palindromic matrix polynomial, even
               matrix polynomial, odd matrix polynomial, Cayley transformation,
               linearization, Hamiltonian matrix, symplectic matrix"
   }

@techreport{Higham:2005:CLM,
  author = "Nicholas J. Higham and D. Steven Mackey and Fran\c{c}oise Tisseur",
  title = "The Conditioning of Linearizations of Matrix Polynomials",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 465",
  pages = "25",
  month = apr,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep465.pdf",
  abstract = "The standard way of solving the polynomial eigenvalue problem of
              degree $m$ in $n\times n$ matrices is to ``linearize'' to a pencil
              in $mn\times mn$ matrices and solve the generalized eigenvalue
              problem. For a given polynomial, $P$, infinitely many
              linearizations exist and they can have widely varying eigenvalue
              condition numbers. We investigate the conditioning of
              linearizations from a vector space $\mathbb{DL}(P)$ of pencils
              recently identified and studied by Mackey, Mackey, Mehl, and
              Mehrmann. We look for the best conditioned linearization and
              compare the conditioning with that of the original polynomial. Two
              particular pencils are shown always to be almost optimal over
              linearizations in $\mathbb{DL}(P)$ for eigenvalues of modulus
              greater than or less than 1, respectively, provided that the
              problem is not too badly scaled and that the pencils are
              linearizations. Moreover, under this scaling assumption, these
              pencils are shown to be about as well conditioned as the original
              polynomial. For quadratic eigenvalue problems that are not too
              heavily damped, a simple scaling is shown to convert the problem
              to one that is well scaled. We also analyze the eigenvalue
              conditioning of the widely used first and second companion
              linearizations. The conditioning of the first companion
              linearization relative to that of $P$ is shown to depend on the
              coefficient matrix norms, the eigenvalue, and the left \ev s of
              the linearization and of $P$. The companion form is found to be
              potentially much more ill conditioned than $P$, but if the 2-norms
              of the coefficient matrices are all approximately 1 then the
              companion form and $P$ are guaranteed to have similar condition
              numbers. Analogous results hold for the second companion form. Our
              results are phrased in terms of both the standard relative
              condition number and the condition number of Dedieu and Tisseur
              for the problem in homogeneous form, this latter condition number
              having the advantage of applying to zero and infinite
              eigenvalues.",
   keywords = "matrix polynomial, matrix pencil, linearization, companion form,
               condition number, homogeneous form, quadratic eigenvalue problem,
               vector space, scaling"
}

@techreport{Mackey:2005:VSL,
  author = "D. Steven Mackey and Niloufer Mackey and Christian Mehl and
            Volker Mehrmann",
  title = "Vector Spaces of Linearizations for Matrix Polynomials",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type =         t-NAREP,
  number = "No. 464",
  pages = "32",
  month = apr,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep464.pdf",
  abstract = "The classical approach to investigating polynomial eigenvalue
              problems is linearization, where the polynomial is converted into
              a larger matrix pencil with the same eigenvalues. For any
              polynomial there are innitely many linearizations with widely
              varying properties, but in practice the companion forms are
              typically used. However, these companion forms are not always
              entirely satisfactory, and linearizations with special properties
              may sometimes be required. In this paper we develop a systematic
              approach to generating large classes of linearizations for matrix
              polynomials. Given a polynomial P, we show how to simply construct
              two vector spaces of pencils that generalize the companion forms
              of P, and prove that almost all of these pencils are
              linearizations for P. Eigenvectors of these pencils are shown to
              be closely related to those of P. A distinguished subspace is then
              isolated, and the special properties of these pencils are
              investigated. These spaces of pencils provide a convenient arena
              in which to look for structured linearizations of structured
              polynomials, as well as to try to optimize the conditioning of
              linearizations, issues to be addressed in further work.",
   keywords = "matrix polynomial, matrix pencil, linearization, strong
               linearization, shifted sum, companion form"
}

@techreport{Lancaster:2005:CFS,
  author = "Peter Lancaster and Leiba Rodman",
  title = "Canonical Forms for Symmetric/Skew-symmetric Real Matrix Pairs
           under Strict Equivalence and Congruence",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 463",
  pages = "77",
  month = mar,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep463.pdf",
  abstract = "A systematic development is made of the simultaneous reduction of
              pairs of quadratic forms over the reals, one of which is
              skew-symmetric and the other is either symmetric or
              skew-symmetric. These reductions are by strict equivalence and by
              congruence, over the reals or over the complex numbers, and
              essentially complete proofs are presented.  Some closely related
              results which can be derived from the canonical forms of pairs of
              symmetric/skew-symmetric real forms are also included. They
              concern simultaneously neutral subspaces, Hamiltonian and
              skew-Hamiltonian matrices, and canonical structures of real
              matrices which are selfadjoint or skew-adjoint in a regular
              skew-symmetric indefinite inner product, and real matrices which
              are skew-adjoint in a regular symmetric indefinite inner product.
              The paper is largely expository, and continues the comprehensive
              account of the reduction of pairs of matrices started in a
              preceding paper.",
   keywords = ""
}

@techreport{Lancaster:2005:ISP,
  author = "Peter Lancaster",
  title = "Inverse Spectral Problems for Semi-simple Damped Vibrating Systems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 462",
  pages = "24",
  month = mar,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep462.pdf",
  abstract = "Computational schemes are investigated for the solution of inverse
              spectral problems for $n\times n$ real systems of the form
              $L(\la)=M\la^2 +D\la +K$. Thus, admissible sets of data concerning
              systems of eigenvalues and eigenvectors are examined and
              procedures for generating associated (isospectral) families of
              systems are developed. The analysis includes symmetric systems,
              systems with mixed real/non-real spectrum, systems with  positive
              definite coefficients, and hyperbolic systems (with real
              spectrum). A one-to-one correspondence between Jordan pairs and
              structure preserving similarities is clarified. An examination of
              complex symmetric matrices is included.",
   keywords = ""
}

@techreport{Hargreaves:2005:EAM,
  author = "Gareth I. Hargreaves and Nicholas J. Higham",
  title = "Efficient Algorithms for the Matrix Cosine and Sine",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "461",
  pages = "17",
  month = feb,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep461.pdf",
  abstract = "Several improvements are made to an algorithm of Higham and Smith
              for computing the matrix cosine. The original algorithm scales the
              matrix by a power of 2 to bring the $\infty$-norm to 1 or less,
              evaluates the [8/8] Pad\'e approximant, then uses the double-angle
              formula $\cos(2A) = 2\cos^2A - I$ to recover the cosine of the
              original matrix. The first improvement is to phrase truncation
              error bounds in terms of $\norm{A^2}^{1/2}$ instead of the (no
              smaller and potentially much larger quantity) $\norm{A}$. The
              second is to choose the degree of the Pad\'e approximant to
              minimize the computational cost subject to achieving a desired
              truncation error. A third improvement is to use an absolute,
              rather than relative, error criterion in the choice of Pad\'e
              approximant; this allows the use of higher degree approximants
              without worsening an a priori error bound. Our theory and
              experiments show that each of these modifications brings a
              reduction in computational cost. Moreover, because the
              modifications tend to reduce the number of double-angle steps they
              usually result in a more accurate computed cosine in floating
              point arithmetic. We also derive an algorithm for computing both
              $\cos(A)$ and $\sin(A)$, by adapting the ideas developed for the
              cosine and intertwining the cosine and sine double angle
              recurrences.",
  keywords = "matrix function, matrix cosine, matrix sine, matrix exponential,
              Taylor series, Pad\'e approximation, Pad\'e approximant,
              double-angle formula, rounding error analysis, Schur--Parlett
              method, MATLAB "
}

@techreport{Lancaster:2005:CFH,
  author = "Peter Lancaster and Leiba Rodman",
  title = "Canonical Forms for Hermitian Matrix Pairs
           under Strict Equivalence and Congruence",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 460",
  pages = "45",
  month = jan,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep460.pdf",
  abstract = "After a brief historical review and an account of the canonical
              forms attributed to Jordan and Kronecker, a systematic development
              is made of the simultaneous reduction of pairs of quadratic forms
              over the complex numbers and over the reals. These reductions are
              by strict equivalence and by congruence and essentially complete
              proofs are presented.  Some closely related results which can be
              derived from the canonical forms are also included. They concern
              simultaneous diagonalization, a new criterion for the existence of
              positive definite linear combinations of a pair of hermitian
              matrices, and the canonical structures of matrices which are
              selfadjoint in an indefinite inner product. ",
   keywords = "",
   note = "To appear in SIAM Review."
}

@techreport{Lancaster:2005:IVS2,
  author = "Uwe Prells and Peter Lancaster",
  title = "Isospectral Vibrating Systems. {Part 2}: Structure Preserving
           Transformations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 459",
  pages = "23",
  month = jan,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep459.pdf",
  abstract = "The study of inverse problems for $n\times n$ systems of the form
              $L(\la):=M\la^2 +D\la +K$ is continued. In this paper it is
              assumed that one vibrating system is specified and the objective
              is to generate isospectral families of systems, i.e. systems which
              reproduce precisely the eigenvalues of the given system together
              with their multiplicities. Two central ideas are developed and
              used, namely, standard triples of matrices, and structure
              preserving transformations.",
  keywords = ""
}

@techreport{Lancaster:2005:IVS1,
  author = "Peter Lancaster",
  title = "Isospectral Vibrating Systems. {Part 1}: The Spectral Method",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 458",
  pages = "17",
  month = jan,
  year = "2005",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep458.pdf",
  abstract = "A study is made of inverse problems for $n\times n$ systems of the
              form $L(\la)=M\la^2 +D\la +K$.  This paper concerns the
              determination of systems in an equivalence class defined by a
              fixed $2n\times 2n$ admissible Jordan matrix, i.e. a class of
              isospectral systems. Constructive methods are obtained for complex
              or real systems with no symmetry constraints. It is also shown how
              isospectral families of complex hermitian matrices can be formed.
              The case of real symmetric matrices is more difficult. Some
              partial solutions are obtained but, in this case, the theory
              remains incomplete. Examples are given.",
  keywords = "",
  note = "To appear in Linear Algebra Appl."
}

@techreport{Mathias:2004:DGE,
  author = "Roy Mathias and Chi-Kwong Li",
  title = "The Definite Generalized Eigenvalue Problem:
           {A} New Perturbation Theory",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 457",
  pages = 51,
  month = oct,
  year = 2004,
  URL =     "http://www.maths.man.ac.uk/~nareports/narep457.pdf",
  abstract = "Let $(A,B)$ be a definite pair of $n \times n$ Hermitian matrices.
              That is, $|x^*Ax| + |x^*Bx| \neq 0$ for all non-zero vectors
              $x \in \IC^n$. It is possible to find an $n \times n$ non-singular
              matrix $X$ with unit columns such that
              $$
	      X^{*} (A + iB) X = {\rm diag}(\alpha_1 + i \beta_1, \dots, \alpha_n + i
	                      \beta_n)
	      $$
	      where $\alpha_j$ and $\beta_j$ are real numbers.  We call the pairs
	      $(\alpha_j, \beta_j)$ {\it normalized generalized eigenvalues} of  the
	      definite pair $(A, B)$.
	      These pairs have not been studied previously.  We rework the perturbation
	      theory for  the eigenvalues and eigenvectors  of the definite
	      generalized eigenvalue problem $\beta Ax = \alpha Bx$ in terms of
	      these normalized generalized eigenvalues and show that they play
	      a crucial rule in obtaining the best possible perturbation  bounds.
              In particular, in existing perturbation bounds, one can replace most
	      instances of the Crawford number
	      \[
	      c(A,B) = \min \{ |x^{*} (A + iB)x|:  x \in \IC^n, x^{*}x = 1 \}
	      \]
	      with the larger quantity
	      \[
	      d_{\min}  = \min\{ |\alpha_j + i\beta_j|: j = 1, \ldots, n \}.
	      \]
              This results in bounds that can be stronger by an arbitrarily large
	      factor. We also give a new measure of the separation of the $j$th
	      eigenvalue from the $k$th:
	      \[
	      |(\alpha_j + i \beta_j) \sin (\arg(\alpha_j + i \beta_j)
	              - \arg(\alpha_k + i \beta_k))|  .
	      \]
	      This asymmetric measure is entirely new, and again
	      results in bounds that can be arbitrarily stronger  than
	      the existing bounds. We show that all but one of our bounds are attainable.
	      We also show that the Crawford number is the infimum of the distance
	      from a definite pencil, {\it a fortiori} diagonalizable, to a
	      non-diagonalizable pair."
      }

@techreport{Shardlow:2004:GED,
  author = "Shardlow, Tony and Yan, Yubin",
  title = "Geometric Ergodicity for Dissipative Particle Dynamics",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 456",
  pages = 29,
  month = oct,
  year = 2004,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep456.pdf",
  abstract = "Dissipative particle dynamics is a system of stochastic
              differential equations modelling complex fluid flows. We consider
              the problem of $N$ particles evolving on the one dimensional
              periodic domain of length $L$ and, if the density of particles is
              large, prove geometric convergence to a unique invariant measure.
              The proof uses minorization and drift arguments, but allows
              elements of the drift and diffusion matrix to have compact support
              where hypoellipticity arguments are not directly available."
      }

@techreport{Shardlow:2004:MES,
  author = "Shardlow, Tony",
  title = "Modified Equations for Stochastic Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 455",
  pages = 12,
  month = oct,
  year = 2004,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep455.pdf",
  abstract = "We describe a backward error analysis for stochastic differential
              equations with respect to weak convergence. Modified equations are
              provided for forward and backward Euler approximations to Ito SDEs
              with additive noise, and extensions to other types of equations
              and approximation discussed."
      }

@techreport{Bini:2004:AMP,
  author = "Dario A. Bini and Nicholas J. Higham and Beatrice Meini",
  title = "Algorithms for the Matrix {$p$}th Root",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 454",
  pages = 26,
  month = aug,
  year = 2004,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep454.pdf",
  abstract = "New theoretical results are presented about the principal matrix
              $p$th root. In particular, we show that the $p$th root is related
              to the matrix sign function and to the Wiener-Hopf factorization,
              and that it can be expressed as an integral over the unit circle.
              These results are used in the design and analysis of several new
              algorithms for the numerical computation of the $p$th root. We
              also analyze the convergence and numerical stability properties of
              Newton's method for the inverse $p$th root. Preliminary
              computational experiments are presented to compare the methods.",
  keywords = "Matrix $p$th root, matrix sign function, Wiener-Hopf
              factorization, Newton's method, Graeffe iteration, cyclic
              reduction, Laurent polynomial"
      }

@techreport{Powell:2004:PFP,
  author = "Catherine Powell",
  title = "Parameter-free $\mathbf{H(div)}$ preconditioning for mixed
           finite element formulation of diffusion problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 453",
  month = aug,
  year = 2004,
  URL =  "http://www.maths.man.ac.uk/~nareports/narep453.ps.gz",
  abstract = "In this paper, we implement the lowest-order Raviart-Thomas scheme
              for generalised diffusion problems and analyse the robustness of
              so-called $H(div)$ preconditioning for the resulting saddle-point
              systems. Properties of the exact version of the suggested
              block-diagonal preconditioner were discussed in some previous
              work,  `Optimal preconditioning for Raviart-Thomas mixed
              formulation of second-order elliptic problems', SIAM J. Matrix
              Anal. Appl., 25, 3(2003). Here, new bounds  are established for
              the eigenvalue spectrum of the  preconditioned system matrix when
              a single V-cycle of any available multigrid scheme is used to
              approximately invert the matrix corresponding to the discrete
              weighted $H(div)$ operator."
}

@techreport{Higham:2004:SSM,
  author = "Nicholas J. Higham",
  title = "The Scaling and Squaring Method for the Matrix Exponential
           Revisited",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 452",
  pages = "15",
  month = jul,
  year = "2004",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep452.pdf",
  abstract = "The scaling and squaring method is the most widely used method for
              computing the matrix exponential, not least because it is the
              method implemented in MATLAB's \texttt{expm} function. The method
              scales the matrix by a power of 2 to reduce the norm to order 1,
              computes a Pad\'e approximant to the matrix exponential, then
              repeatedly squares to undo the effect of the scaling. We give a
              new backward error analysis of the method (in exact arithmetic)
              that employs sharp bounds for the truncation errors and leads to
              an implementation of essentially optimal efficiency. We also give
              new rounding error analysis that shows the computed Pad\'e
              approximant of the scaled matrix to be highly accurate. For IEEE
              double precision arithmetic the best choice of degree of Pad\'e
              approximant turns out to be 13, rather than the 6 or 8 used by
              previous authors. Our implementation of the scaling and squaring
              method always requires at least two fewer matrix multiplications
              than \texttt{expm} when the matrix norm exceeds 1, which can
              amount to a 37\% saving in the number of multiplications, and it
              is typically more accurate owing to the fewer required squarings.
              We also investigate a different scaling and squaring algorithm
              proposed by Najfeld and Havel that employs a Pad\'e approximation
              to the function $x \coth(x)$. This method is found to be
              essentially a variation of the standard one with weaker supporting
              error analysis.",
  keywords = "matrix function, matrix exponential, Pad\'e approximation, matrix
              polynomial evaluation, scaling and squaring method, MATLAB, expm,
              backward error analysis, performance profile"
      }

@techreport{Roberts:2004:SVO,
  author = "Leonid E. Shaikhet and Jason A. Roberts",
  title = "Stochastic {Volterra} Integro-Differential Equations:
           Stability and Numerical Methods",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 450",
  pages = 38,
  month = jun,
  year = 2004,
  url = "http://www.maths.man.ac.uk/~nareports/narep450.pdf",
  abstract = "We consider the reliability of some numerical methods in
              preserving the stability properties of the linear stochastic
              functional differential equation $$\dot{x}(t)=\alpha
              x(t)+\beta\int^t_0x(s)ds+\sigma x(t-\tau)\dot{W}(t),$$ where
              $\alpha, \beta, \sigma, \tau \geq 0$ are real constants, and
              $W(t)$ is a standard Wiener process. We adopt the shorthand
              notation of $\dot{x}(t)$ to represent the differential $dx(t)$
              etc. Our choice of test equation is a stochastic perturbation of
              the classical deterministic Brunner \& Lambert test equation for
              $\sigma=0$ and so our investigation may be thought of as an
              extension of their work. Sufficient conditions for the asymptotic
              mean square stability of solutions to both the differential
              equation and discrete analogues are derived using the general
              method of Lyapunov functionals construction proposed by
              Kolmanovskii \& Shaikhet which has previously been successfully
              employed for deterministic and stochastic differential and
              difference equations with delay. The areas of the regions of
              asymptotic stability for each $\theta$ method, indicated by the
              sufficient conditions for the discrete system, are shown to be
              equal and we show that an upper bound can be put on the time-step
              parameter for the numerical method fo which the system is
              asymptotically mean-square stable. We illustrate our results by
              means of numerical experiments and various stability diagrams. We
              examine the extent to which the continuous system can tolerate
              stochastic perturbations before losing its stability properties
              and we illustrate how one may accurately choose a numerical method
              to preserve the stability properties of the original problem in
              the numerical solution. Our numerical experiments also indicate
              that the quality of the sufficient conditions is very high."
}

@techreport{MCCM:2004:AR3,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 2003",
  type =         t-NAREP,
  number =       "No. 449",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "25",
  month =        may,
  year =         "2004",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep449.pdf",
  }

@techreport{Baker:2004:HAL,
  author = "Christopher T. H. Baker and Evelyn Buckwar",
  title = "A note on {Halanay}-type stability results for the
           {$\theta$}--{Maruyama} method for stochastic delay differential
            equations",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 448",
  pages = 13,
  month = apr,
  year = 2004,
  url = "http://www.maths.man.ac.uk/~nareports/narep448.pdf",
  abstract = "Using an approach that has its origins in work of Halanay, we
              consider stability in  mean square of numerical solutions obtained
              from the {\em $\theta$--Maruyama discretization} of a test
              stochastic delay differential equation \[ dX(t) = \{f(t) -\alpha
            	X(t) + \beta X(t-\tau)\} {\rd}t + \{g(t) + \eta \ X(t) + \mu
              X(t-\tau) \} \ {\rd}W(t) \] interpreted in the It\^o sense, where
              $W(t)$ denotes a Wiener process (white noise). Our objective in
              this report is limited to  demonstrating that we may use
              techniques advanced in a recent report by Baker and Buckwar to
              obtain criteria for asymptotic and exponential stability, in mean
              square, for the solutions of a recurrence \[ {\widetilde X}_{n+1}
              - {\widetilde X}_n = \theta\{f_{n+1} -\A h {\widetilde X}_{n+1} +
              \B h {\widetilde X}_{n+1-N}\} + (1-\theta)\{f_n -\A h{\widetilde
              X}_{n} + \B h {\widetilde X}_{n-N}\}\] \[  + (g_n + \eta \sqrt{h}
              \, {\widetilde X}_n + \mu \sqrt{h} \, {\widetilde X}_{n-N_{}})\
              \xi_n \] where $\xi_n \in N(0,1)$ ($\xi_n$ is normally distributed
              with zero mean and variance unity).  This is the $\theta$-Maruyama
              formula; recognizable cases arise from taking $\theta = 0$,
              $\theta = \half$ or $\theta = 1$; for convenience we assume $
              \theta \in [0,1] $."
}

@techreport{Hargreaves:2004:CCN,
   author = "Gareth I. Hargreaves",
   title = "Computing the Condition Number of Tridiagonal and
            Diagonal-Plus-Semiseparable Matrices in Linear Time",
   institution = inst-MCCM,
   address = inst-MCCM:adr,
   type = t-NAREP,
   number = "No. 447",
   pages = 19,
   month = apr,
   year = 2004,
   url = "http://www.maths.man.ac.uk/~nareports/narep447.pdf",
   abstract = "For an $n \times n$ tridiagonal matrix we exploit the structure of
   its  QR factorization to devise two new algorithms for computing the 1-norm
   condition number in $O(n)$ operations.  The algorithms avoid
   underflow and overflow, and are simpler than existing algorithms since
   tests are not required for degenerate cases.  An error analysis of the
   first algorithm is given, while the second algorithm is shown to be
   competitive in speed with existing algorithms.  We then turn our
   attention to an $n \times n$ diagonal-plus-semiseparable matrix, $A$,
   for which several algorithms have recently been developed to solve
   $Ax=b$ in $O(n)$ operations.  We again exploit the QR factorization of
   the matrix to present an algorithm that computes the 1-norm condition
   number in $O(n)$ operations. ",
}

@techreport{Higham:2004:FPM,
  author = "Nicholas J. Higham and D. Steven Mackey and Niloufer Mackey and
            Fran\c{c}oise Tisseur",
  title = "Functions Preserving Matrix Groups
           and Iterations for the Matrix Square Root",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "446",
  pages = 29,
  month = mar,
  year = "2004",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep446.pdf",
  abstract = "For any matrix automorphism group $\G$ associated with a bilinear
              or sesquilinear form, Mackey, Mackey, and Tisseur have recently
              shown that the matrix sign decomposition factors of $A\in\G$ also
              lie in $\G$; moreover, the polar factors of $A$ lie in $\G$ if the
              matrix of the underlying form is unitary. Groups satisfying the
              latter condition include the complex orthogonal, real and complex
              symplectic, and pseudo-orthogonal groups. This work is concerned
              with exploiting the structure of $\G$ when computing the polar and
              matrix sign decompositions of matrices in $\G$. We give sufficient
              conditions for a matrix iteration to preserve the group structure
              and show that a family of globally convergent rational
              Pad\'e-based iterations of Kenney and Laub satisfy these
              conditions. The well-known scaled Newton iteration for computing
              the unitary polar factor does not preserve group structure, but we
              show that the approach of the iterates to the group is precisely
              tethered to the approach to unitarity, and that this forces a
              different and exploitable structure in the iterates. A similar
              relation holds for the Newton iteration for the matrix sign
              function. We also prove that the number of iterations needed for
              convergence of the structure-preserving methods can be precisely
              predicted by running an associated scalar iteration. Numerical
              experiments are given to compare the cubically and quintically
              converging iterations with Newton's method and to test stopping
              criteria. The overall conclusion is that the structure-preserving
              iterations and the scaled Newton iteration are all of practical
              interest, and which iteration is to be preferred is
              problem-dependent.",
  keywords = "automorphism group, bilinear form, sesquilinear form, scalar
              product, adjoint, Fr\'echet derivative, stability analysis,
              perplectic matrix, pseudo-orthogonal matrix, generalized polar
              decomposition, matrix sign function, matrix $p$th root, matrix
              square root, structure preservation, matrix iteration, Newton
              iteration"
}

@techreport{Lancaster:2004:OPM,
  author = "Peter Lancaster and Panayiotis Psarrakos",
  title = "On the Pseudospectra of Matrix Polynomials",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 445",
  pages = "17",
  month = feb,
  year = "2004",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep445.pdf",
  abstract = "The pseudospectra of a matrix polynomial $P(\la)$ are sets of
              complex numbers that are eigenvalues of matrix polynomials which
              are near to $P(\la)$, i.e., their coefficients are within some
              fixed magnitude of the coefficients of $P(\la)$. Pseudospectra
              provide important insights into the sensitivity of eigenvalues
              under perturbations, and have several applications. First,
              qualitative properties concerning boundedness and connected
              components of pseudospectra are obtained. Then an accurate
              continuation algorithm for the numerical determination of the
              boundary of pseudospectra of matrix polynomials is devised and
              illustrated. This algorithm is based on a prediction-correction
              scheme.",
keywords = "matrix polynomial, eigenvalue, singular value, perturbation,
             $\eps$-pseudospectum, boundary, stability"
}

@techreport{Baker:2004:INP,
  author = "Christopher T.H.Baker and Eugene I. Parmuzin",
  title = "Identification of the initial function for delay differential
           equations: Part III: Nonlinear DDEs \& computational results.",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 444",
  pages = 27,
  month = mar,
  year = 2004,
  url = "http://www.maths.man.ac.uk/~nareports/narep444.pdf",
  abstract = "In this report we consider the ``data assimilation problem'' for
              nonlinear delay differential equations. This problem consists of
              finding an initial function that gives rise to a solution of a
              given nonlinear delay differential equation, which is a close fit
              to observed data. A r\^ole for adjoint equations and fundamental
              solutions in the nonlinear case is established. The
              ``Pseudo-Newton'' method to solve the ``data assimilation
              problem'' for nonlinear DDE is presented. This is an extension of
              results obtain in Nareport No.431 for linear delay differential
              equations."
}

@techreport{Baker:2004:IDP,
  author = "Christopher T. H. Baker and Eugene I. Parmuzin",
  title = "Identification of the initial function for delay differential
           equations: Part II: A discrete analogue of the data assimilation
           problem \& computational results.",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 443",
  pages = 29,
  month = mar,
  year = 2004,
  url = "http://www.maths.man.ac.uk/~nareports/narep443.pdf",
  abstract = "In this report we consider a difference analogue of the ``data
             assimilation problem'' for delay differential equations (DDEs)
             presented in the first part (Numerical Analysis Report No. 431).
             The problem consists of finding an initial function that gives rise
             to a solution of a discretized DDE, which is a close fit to
             observed data. A r\^ole for adjoint equations and fundamental
             solutions in the discrete case is established, and related discrete
             integral equations (summation equations) are obtained. The final
             part of  the report consists of a report of numerical experiments
             that  demonstrate the performance of the algorithm."
}

@techreport{Lucas:2004:LSC,
  author = "Craig Lucas",
  title = "{LAPACK}-Style Codes for Level 2 and 3
           Pivoted {Cholesky} Factorizations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 442",
  pages = "31",
  month = feb,
  year = "2004",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep442.pdf",
  abstract = "Fortran 77 codes exist in LAPACK for computing the Cholesky
              factorization of a symmetric positive definite matrix using
              Level 2 and 3 BLAS. In LINPACK there is a Level 1 BLAS
              routine for computing the Cholesky factorization with complete
              pivoting of a symmetric positive semi-definite matrix. We present
              two  new algorithms and Fortran 77 LAPACK-style codes for
              computing this pivoted factorization: one using Level 2 BLAS
              and one using Level 3 BLAS. We show that on modern machines the
              new codes can be many times faster than the LINPACK code. Also,
              with a new  stopping criterion they provide more reliable rank
              detection and can have a smaller normwise backward error.",
  keywords = "pivoted Cholesky factorization, semi-definite matrices,
              LAPACK, Level 3 BLAS."
}

@techreport{Davies:2004:SCM,
  author = "Philip I. Davies",
  title = "Structured Conditioning of Matrix Functions",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 441",
  pages = "33",
  month = jan,
  year = "2004",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep441.pdf",
  abstract = "The existing theory of conditioning for matrix functions
              $f(X) \colon \C^{\nbyn} \to \C^{\nbyn}$ does not cater for structure
              in the matrix $X$. An extension of this theory is presented in which
              when $X$ has structure, all perturbations of $X$ are required to
              have the same structure. We consider two classes of structured
              matrices, those comprising the Jordan algebra $\J$ and the Lie
              algebra $\L$ associated with a nondegenerate bilinear or sesquilinear
              form on $\R^n$ or $\C^n$. Examples of such classes are the symmetric,
              skew-symmetric, Hamiltonian and skew-Hamiltonian matrices. Structured
              condition numbers are defined for these two classes. Under certain
              conditions on the underlying scalar product, explicit representations
              are given for the structured condition numbers. Comparisons between
              the unstructured and structured condition numbers are then made.
              We show that when the underlying scalar product is a sesquilinear
              form, there is no difference between the values of the two condition
              numbers for (i) all functions of $X \in \J$,  and (ii) odd and even
              functions of $X \in \L$. When the underlying scalar product is a
              bilinear form then equality is not guaranteed in all these cases.
              Where equality is not guaranteed, bounds are obtained for the ratio
              of the unstructured and structured condition numbers.",
  keywords = "Matrix functions, Fr\'echet derivative, Condition numbers, Bilinear
              forms, Sesquilinear forms, Structured Matrices, Jordan algebra,
              Symmetric matrices, Lie algebra, Skew-symmetric matrices"
}

@techreport{Higham:2003:NSB,
  author = "Nicholas J. Higham",
  title = "The Numerical Stability of Barycentric {Lagrange} Interpolation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 440",
  pages = "10",
  month = dec,
  year = "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep440.pdf",
  abstract = "The Lagrange representation of the interpolating polynomial can be
              rewritten in two more computationally attractive forms: a modified
              Lagrange form and a barycentric form. We give an error analysis of
              the evaluation of the interpolating polynomial using these two
              forms. The modified Lagrange formula is shown to be backward
              stable. The barycentric formula has a less favourable error
              analysis, but is forward stable for any set of interpolating
              points with a small Lebesgue constant. Therefore the barycentric
              formula can be significantly less accurate than the modified
              Lagrange formula only for a poor choice of interpolating points.
              This analysis provides further weight to the argument of Berrut
              and Trefethen that barycentric Lagrange interpolation should be
              the polynomial interpolation method of choice.",
  keywords = "polynomial interpolation, Lagrange interpolation, barycentric
              formula, rounding error analysis, backward error, forward error,
              Lebesgue constant"
}

@techreport{Buckwar:2003:WAS,
  author = "Buckwar, Evelyn and Shardlow, Tony",
  title = "Weak Approximation of Stochastic Delay Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 439",
  month = dec,
  year = "2003",
  URL = "http://www.maths.man.ac.uk/~nareports/narep439.pdf"
   }

@techreport{Shardlow:2003:NWE,
  author = "Shardlow, Tony",
  title = "Nucleation of Waves in Excitable Media by Noise",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 438",
  month = dec,
  year = "2003",
  URL = "http://www.maths.man.ac.uk/~nareports/narep438.pdf"
   }

@techreport{Shardlow:2003:NSS,
  author = "Shardlow, Tony",
  title = "Numerical Simulation of Stochastic {PDE}s for Excitable Media",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 437",
  month = dec,
  year = "2003",
  URL = "http://www.maths.man.ac.uk/~nareports/narep437.pdf"
   }

@techreport{Davies:2003:CFM,
  author = "Philip I. Davies and Nicholas J. Higham",
  title = "Computing {$f(A)b$} for Matrix Functions {$f$}",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 436",
  pages = "10",
  month = nov,
  year = "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep436.pdf",
  abstract = "For matrix functions $f$ we investigate how to compute a
              matrix-vector product $f(A)b$ without explicitly
              computing $f(A)$. A general method is described that
              applies quadrature to the matrix version of the Cauchy
              integral theorem. Methods specific to the logarithm,
              based on quadrature, and fractional matrix powers, based
              on solution of an ordinary differential equation initial
              value problem, are also described."
  }

@techreport{Ford:2003:SGB,
  author = "Judith M. Ford and Neville J. Ford and John Wheeler",
  title = "Simulation of Grain Boundary Creep: Analysis of Some New Numerical
           Techniques",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 435",
  month = oct,
  year = "2003",
  URL = "http://www.maths.man.ac.uk/~nareports/narep435.pdf"
   }

@techreport{Ford:2003:SLS,
  author = "Judith M. Ford and Eugene E. Tyrtyshnikov",
  title = "Solving Linear Systems using Wavelet Compression Combined with
           {Kronecker} Product Approximation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 434",
  month = oct,
  pages = "11",
  year = "2003",
  URL = "http://www.maths.man.ac.uk/~nareports/narep434.pdf",
  note = "Submitted to special edition of Numerical Algorithms, Proceedings of
          Quatrième séminaire sur l'algorithmique numérique appliquée aux
          problèmes industriels, Calais, 2003"
}

@techreport{Alam:2003:OCN,
  author = "Rafikul Alam",
  title = "On the Construction of Nearest Defective Matrices to a Normal
           Matrix",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 433",
  pages = "6",
  month = sep,
  year = "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep433.pdf",
  abstract = "Let $A$ be an $n$-by-$n$ normal matrix with $n$ distinct
              eigenvalues. We describe a simple procedure to construct a
              defective matrix whose distance from $A,$ among all the defective
              matrices, is the shortest. Our construction requires only an
              appropriate pair of eigenvalues of $A$ and their corresponding
              eigenvectors. We also show that the nearest defective matrices
              influence the evolution of spectral projections of $A$.",
  keywords = "Defective matrix, multiple eigenvalue,
              pseudospectra, eigenvector, spectral projection"
}

@techreport{Mackey:2003:SFS,
  author = "D. Steven Mackey and Niloufer Mackey and Fran\c{c}oise Tisseur",
  title =  "Structured Factorizations in Scalar Product Spaces",
  type =         t-NAREP,
  number =       "No. 432",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "",
  month =        "",
  year = "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep432.pdf",
  note = "In preparation."
}

@techreport{Baker:2004:ICP,
  author = "Christopher T. H. Baker and Eugene I. Parmuzin",
  title = "Identification of the initial function for delay differential
           equations: Part I: The continuous problem \& an integral equation
           analysis.",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 431",
  pages = 25,
  month = jan,
  year = 2004,
  url = "http://www.maths.man.ac.uk/~nareports/narep431.pdf",
  abstract = "In this report, which is the first of three related reports, we
              consider the ``data assimilation problem'' for delay differential
              equations. This problem consists of finding an initial function
              that gives rise to a solution which is a close fit to observed
              data. A r\^ole for adjoint equations and fundamental solutions is
              established, and related integral equations are obtained."
}

@techreport{Baker:2004:GTV,
  author = "Christopher T. H. Baker and Eugene I. Parmuzin",
  title = "A Guided Tour of Variation of Parameters Formulae for
           Continuous and Discretized {DDE}s.",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 430",
  pages = 21,
  month = jan,
  year = 2004,
  url = "http://www.maths.man.ac.uk/~nareports/narep430.pdf",
  abstract = "In this report we discuss variation of parameters formulae ($i$)
              for systems of linear delay differential equations (DDEs), $$
              \frac{\dss d y(t)}{\dss dt} -A(t) y(t)   -B(t) y(t-\tau) = f(t),
              \; \mbox{for}\  t\in [0,T], $$ subject to an initial condition
              $y(t) =\varphi(t) \quad \mbox{for} \ t \in [-\tau,0]$,  and ($ii$)
              for the corresponding discretized equations that result from the
              application of Euler methods. The explicit Euler method for the
              above DDE, using a step $h = \tau / N$ with $T = Kh$ and $N,K \in
              \Nset$, yields the equations $$ {\widetilde{y}}_{n+1} -
              {\widetilde{y}}_{n} -A_n {\widetilde{y}}_n -B_n
              {\widetilde{y}}_{n-N} =f_n, \quad \mbox{for}\ n \in \{0,1, \cdots,
              K-1 \}, $$ subject to the initial condition $ {\widetilde{y}}_n =
              \varphi_n \; \mbox{for} \ -N \leq n \leq 0, $ where $B_n = B(nh)$,
              $A_n = A(nh)$ $f_n = f(nh)$ and $\varphi_n = \varphi(nh)$. Of
              interest is the form of the solution to the DDE or its discretized
              form, in a representation that depends upon the inhomogeneous term
              $f$ or the initial function $\varphi$. Various ways to derive
              variation of parameters formulae will be pursued. In the
              continuous case, the fundamental solution can be obtained as a
              solution of a DDE, a related integral equation, by seeking a
              generalized Green's function or by considering adjoint
              equations.(A number of results related to the continuous case are
              presented in the literature and these are recalled in an
              appendix.)  We look for analogous results in the discretized case,
              drawing on parallels in linear algebra."
}

@techreport{MCCM:2003:AR2,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 2002",
  type =         t-NAREP,
  number =       "No. 429",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "19",
  month =        jun,
  year =         "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep429.pdf",
  }

@techreport{Tisseur:2003:EAM,
  author = "Dario A. Bini and Luca Gemignani and Fran\c{c}oise Tisseur",
  title = "The {Ehrlich-Aberth} Method for the Nonsymmetric
           Tridiagonal Eigenvalue Problem",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 428",
  pages = "23",
  month = jun,
  year = "2003",
  URL =  "ftp://ftp.ma.man.ac.uk/pub/narep/narep428.pdf",
  abstract = "An algorithm based on the Ehrlich-Aberth iteration is presented
              for the computation of the zeros of $p(\l)=\det(T-\l I)$, where
              $T$ is an irreducible tridiagonal matrix. The algorithm requires
              the evaluation of $p(\l)/p'(\l)=-1/\trace(T-\l I)^{-1}$, which is
              done here by exploiting the QR factorization of $T-\l I$ and the
              semiseparable structure of $(T-\l I)^{-1}$. Two choices of the
              initial approximations are considered; the most effective relies
              on a divide-and-conquer strategy, and some results motivating this
              strategy are given. A Fortran 95 module implementing the algorithm
              is provided and numerical experiments that confirm the
              effectiveness and the robustness of the approach are presented. In
              particular, comparisons with the LAPACK subroutine \texttt{dhseqr}
              show that our algorithm is faster for large dimensions.",
  keywords = "nonsymmetric eigenvalue problem,
              symmetric indefinite generalized eigenvalue problem,
              tridiagonal matrix, root finder, QR decomposition,
              divide and conquer."
    }

@techreport{Mackey:SPA:2003,
  author = "D. Steven Mackey and Niloufer Mackey and Daniel M. Dunlavy",
  title = "Structure Preserving Algorithms for Perplectic Eigenproblems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 427",
  pages = "30",
  month = may,
  year = "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep427.pdf"
}

@techreport{Higham:2003:CPD,
  author = "Nicholas J. Higham and D. Steven Mackey and Niloufer Mackey and
            Fran\c{c}oise Tisseur",
  title = "Computing the Polar Decomposition and the Matrix Sign Decomposition
           in Matrix Groups",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 426",
  pages = "15",
  month = apr,
  year = "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep426.pdf",
  abstract = "For any matrix automorphism group $\G$ associated with a
              bilinear or sesquilinear form, Mackey, Mackey, and Tisseur have
              recently shown that the matrix sign decomposition factors of
              $A\in\G$ also lie in $\G$; moreover, the polar factors of $A$
              lie in $\G$ if the matrix of the underlying form is unitary.
              Groups satisfying the latter condition include the complex
              orthogonal, real and complex symplectic, and pseudo-orthogonal
              groups. This work is concerned with exploiting the structure of
              $\G$ when computing the polar and matrix sign decompositions of
              matrices in $\G$. We give sufficient conditions for a matrix
              iteration to preserve the group structure and show that a family
              of globally convergent rational Pad\'e-based iterations of
              Kenney and Laub satisfy these conditions. The well-known scaled
              Newton iteration for computing the unitary polar factor does not
              preserve group structure, but we show that the approach of the
              iterates to the group is precisely tethered to the approach to
              unitarity, and that this forces a different and exploitable
              structure in the iterates. A similar relation holds for the
              Newton iteration for the matrix sign function. We also prove
              that the number of iterations needed for convergence of the
              structure-preserving methods can be precisely predicted by
              running an associated scalar iteration. Numerical experiments
              are given to compare the cubically and quintically converging
              iterations with Newton's method and to test stopping criteria.
              The overall conclusion is that the structure-preserving
              iterations and the scaled Newton iteration are all of practical
              interest, and which iteration is to be preferred is
              problem-dependent.",
  keywords = "automorphism group, bilinear form, sesquilinear form, adjoint,
              complex orthogonal matrix, symplectic matrix, perplectic matrix,
              pseudo-orthogonal matrix, polar decomposition, matrix sign
              decomposition, structure preservation, matrix iteration, Newton
              iteration, convergence tests"
}

@techreport{Baker:2003:MDC,
  author = "Christopher T. H. Baker and Gennadii A. Bocharov  and Christopher
            A. H. Paul and Fathalla A. Rihan",
  title = "Models with Delays for Cell Population Dynamics:
           Identification, Selection and Analysis. Part I: Computational
           Modelling with Functional Differential Equations:
           Identification, Selection, and Sensitivity",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 425",
  pages = "28",
  month = feb,
  year = "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep425.pdf",
  abstract = "This is the first in a series of reports in which we discuss aspects
        of the identification of, discrimination between, and sensitivity of,
        mathematical models, and we illustrate with some examples.
        Mathematical models based upon certain types of differential
        equations, functional differential equations, or systems of
        such equations, are often employed to represent the dynamics
        of natural, in particular biological, phenomena.  We present
        some of the principles underlying the choice of a methodology
        for the computational development of quantitatively consistent
        models, using scientifically meaningful parameters, based on
        observational data.
        Exploiting the link between a (weighted) least-squares fit,
        information theory, and maximum likelihood, we can select an
        appropriate objective function which is minimized over a set
        of parameters. To compare the best fit models from amongst a
        hierarchy of possible models we can employ indicators that can
        be regarded as taking parsimony into account. A sensitivity
        analysis provides feedback on the covariances of the
        parameters in relation to the best fit and on whether we need
        to include all the parameters in the model.  We attempt to
        formulate a methodology for assessing the adequacy of
        mathematical models, founded on this approach.
        Much of our thinking is influenced by discussions that are
        dispersed amongst the existing literature. However, we believe
        that the value of our material is more than the sum of that of
        its separate parts. Additionally, in a number of places we
        extend or supplement or reformulate existing results. The
        models that we consider generate infinite-dimensional
        dynamical systems.",
  keywords = " Cell growth, model selection, neutral delay
         differential equation, parameter estimation, first \& second
         order sensitivity analysis, nonlinear bias, objective
         functions, methodologies."
}

@techreport{Higham:2003:SCS,
  author = "Nicholas J. Higham and Mihail Konstantinov and Volker Mehrmann
            and Petko Petkov",
  title = "Sensitivity of Computational Control Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 424",
  pages = "32",
  month = mar,
  year = "2003",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep424.pdf",
  abstract = "It is well-known that many factors contribute to the accurate and
              efficient numerical solution of mathematical problems such as
              those arising in computational control system design. In simple
              terms these are the arithmetic of the machine on which the
              calculations are carried out, sensitivity (or conditioning) of the
              mathematical model to small changes of the data and the numerical
              stability of the algorithms. It happens quite often that these
              concepts are confused. We define these concepts and demonstrate
              some of the subtleties that often lead to confusion. In particular
              we demonstrate with several examples what may happen when a
              problem is modularized, i.e., split into subproblems for which
              computational modules are available. For three classical problems
              in computational control, pole placement, linear quadratic control
              and optimal $H_\infty$ control, we then discuss the conditioning
              of the problems and point out sources of difficulties. We give
              some ill-conditioned examples for which even numerically stable
              methods fail. We also stress the need for condition and error
              estimators that supplement the numerical algorithm and inform the
              user about potential or actual difficulties, and we explain what
              can be done to avoid these difficulties.",
  keywords = "numerical stability, machine
               arithmetic, pole placement, linear quadratic control,
               algebraic Riccati equation, $H_\infty$ control"
}

@techreport{Mackey:2003:ODS,
  author = "D. Steven Mackey and Niloufer Mackey",
  title = "On the Determinant of Symplectic Matrices",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 422",
  pages = "12",
  month = feb,
  year = "2003",
  URL =  "http://www.maths.man.ac.uk/~nareports/narep422.pdf",
  abstract = "A collection of new and old proofs showing that the determinant of
              any symplectic matrix is $+1$ is presented. Structured
              factorizations of symplectic matrices play a key role in several
              arguments. A constructive derivation of the symplectic analogue of
              the Cartan-Dieudonn\'{e} theorem is one of the new proofs in this
              essay.",
  keywords = "symplectic, determinant, bilinear form, skew-symmetric,
              structure-preserving factorizations, symmetries, transvections,
              $\G$-reflector, Pfaffian."
   }

@techreport{Tisseur:04:SFSPS,
  author = "D. Steven Mackey and Niloufer Mackey and Fran\c{c}oise Tisseur",
  title = "Structured Factorizations in Scalar Product Spaces",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 421",
  pages = "30",
  month = Nov,
  year = "2004",
  URL =  "http://www.maths.man.ac.uk/~nareports/narep421.pdf",
  abstract = "
             Let $A$ belong to an automorphism group, Lie algebra or Jordan
             algebra of a scalar product. When $A$ is factored, to what extent
             do the factors inherit structure from $A$? We answer this question
             for the principal matrix square root, the matrix sign
             decomposition, and the polar decomposition. For general $A$, we
             give a simple derivation and characterization of a particular
             generalized polar decomposition, and we relate it to other such
             decompositions in the literature. Finally, we study
             eigendecompositions and structured singular value decompositions,
             considering in particular the structure in eigenvalues,
             eigenvectors and singular values that persists across a wide range
             of scalar products.
             A key feature of our analysis is the identification of two
             particular classes of scalar products, termed unitary and
             orthosymmetric, which serve to unify assumptions for the existence
             of %structured factors and structured factorizations. A variety of
             different characterizations of these scalar product classes is
             given.",
  keywords =  "automorphism group, Lie group, Lie algebra, Jordan algebra,
             bilinear form, sesquilinear form, scalar product, indefinite inner
             product, orthosymmetric, adjoint, factorization, symplectic,
             Hamiltonian, pseudo-orthogonal, polar decomposition, matrix sign
             function, matrix square root, generalized polar decomposition,
             eigenvalues, eigenvectors, singular values, structure
             preservation."
    }

@techreport{Tisseur:2003:GRS,
  author = "D. Steven Mackey and Niloufer Mackey
            and Fran\c{c}oise Tisseur",
  title = "$\mathbb{G}$-Reflectors in Scalar Product Spaces",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 420",
  pages = "24",
  month = feb,
  year = "2003",
  URL =  "http://www.maths.man.ac.uk/~nareports/narep420.ps.gz",
  abstract = " We characterize the analogues of Householder reflectors in matrix
               groups associated with scalar products. Examples of such groups
               include the symplectic and pseudo-unitary groups. We also
               precisely delimit the mapping capabilities of these Householder
               analogues: given a matrix group $\G$ and vectors $x, y$, we give
               necessary and sufficient conditions for the existence of a
               Householder-like analogue $G \in \G$ such that $Gx = y$. When $G$
               exists, we show how it can be constructed from $x$ and $y$.",
  keywords =  "scalar product, bilinear, sesquilinear, isotropic, Householder
               reflection, hyperbolic transformation, symmetries, transvections,
               symplectic, pseudo-unitary, structure-preserving."
}

@techreport{Tisseur:2003:STS,
  author = "D. Steven Mackey and Niloufer Mackey and Fran\c{c}oise Tisseur",
  title = "Structured Tools for Structured Matrices",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 419",
  pages = "36",
  month = feb,
  year = "2003",
  URL =  "http://www.maths.man.ac.uk/~nareports/narep419.pdf",
  abstract = "We present an extensive and unified collection of
              structure-preserving transformations, organized for easy
              reference. The structures we work with arise in the context of a
              non-degenerate bilinear or sesquilinear form on $\R^n$ or $\C^n$.
              We construct a variety of transformations belonging to the
              automorphism groups of these forms, that imitate the action of
              Givens rotations, Householder reflectors, and Gauss
              transformations. We also describe transformations for performing
              structured scaling actions. The matrix groups considered in this
              paper are the complex orthogonal, real, complex and conjugate
              symplectic, real perplectic, real and complex pseudo-orthogonal,
              and pseudo-unitary groups. In addition to deriving new
              transformations, this paper collects and unifies existing
              structure-preserving tools.",
  keywords = "structured matrices, matrix groups, Givens rotations, Householder
              reflections, complex orthogonal, symplectic, complex symplectic,
              conjugate symplectic, real perplectic, pseudo-orthogonal, complex
              pseudo-orthogonal, pseudo-unitary, scalar product, bilinear form,
              sesquilinear form, automorphism groups, Jordan algebra, Lie
              algebra."
}

@techreport{Ford:2003:CKP,
author = "Judith M. Ford and Eugene E. Tyrtyshnikov",
title = "Combining {Kronecker} Product Approximation with Discrete Wavelet
         Transforms to Solve Dense, Function-Related Linear Systems",
institution = inst-MCCM,
address = inst-MCCM:adr,
type =    t-NAREP,
number = "No. 418",
month = jan,
year = 2003,
URL = "http://www.maths.man.ac.uk/~nareports/narep418.pdf"
}

@techreport{Baker:2003:NDE,
  author = "Christopher T. H. Baker and Christopher A. H. Paul",
  title = "Piecewise Continuous Solutions of Neutral Delay
           Differential Equations",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 417",
  pages = 35,
  month = Dec,
  year = 2003,
  url =  "http://www.maths.man.ac.uk/~nareports/narep417.pdf",
  abstract = "Examples of neutral delay differential equations (NDDEs) are \[
              \frac{d}{dt} u(t) = F\big(t,u(t),u(t-\tau),\frac{d}{dt}
              u(t-\tau)\big) \quad (t \ge t_0), \] for $\tau > 0$, and (a form
              frequently called ``Hale's form'') \[ \frac{d}{dt}
              \{u(t)-g\big(t,u(t),u(t-\tau)\big)\} = f\big(t,u(t),u(t-\tau)\big)
              \quad (t \ge t_0) \] which is an NDDE when $g$ does not vanish
              identically.  A particular solution of such an equation is defined
              by initial data on $[-\tau,0]$, and we observe in this report
              that, when the problem is suitably interpreted, the solutions of
              an NDDE need not be continuous. This fact, largely overlooked in
              the numerical analysis literature (many of the texts on the
              analysis of NDDEs admit a function as a solution only if it is
              continuous), presents some difficulties when considering numerical
              methods for approximating solutions. Our discussion of NDDEs
              points to the need to be precise in the definition of the term
              ``solution'', and to consider questions of \emph{existence},
              \emph{uniqueness} and \emph{smoothness} of any solutions of NDDEs.
              Such questions are a precursor to addressing the problem of
              constructing numerical methods. Numerical codes for solving
              differential equations are generally based on the assumption that
              a solution is continuous, and our discussion will indicate
              measures that may be taken if one seeks a discontinuous solution
              of an NDDE. First we show that the difficulties can sometimes be
              alleviated by \emph{perturbing the initial function} to yield a
              continuously differentiable solution. Alternatively,
              \emph{tracking discontinuities} and converting the problem to a
              \emph{delay differential algebraic equation} (DDAE) -- followed
              either by direct solution of the DDAE or by treatment after
              regularizing the DDAE to a singularly perturbed problem -- are
              discussed in the Appendices. The investigation of the simple NDDE
              that opens our discussion is followed by a number of results for
              general NDDEs."
}

@techreport{Hargreaves:2002:IAM,
  author = "Gareth I. Hargreaves",
  title = "Interval Analysis in {MATLAB}",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type =    t-NAREP,
  number = "No. 416",
  pages = 49,
  month = dec,
  year = 2002,
  url = "http://www.maths.man.ac.uk/~nareports/narep416.pdf",
  abstract = "The introduction of fast and efficient software for interval
              arithmetic, such as the MATLAB toolbox INTLAB, has resulted in the
              increased popularity of the use of interval analysis.  We give an
              introduction to interval arithmetic and explain how it is
              implemented in the toolbox INTLAB.  A tutorial is provided for
              those who wish to learn how to use INTLAB.  We then focus on the
              interval versions of some important problems in numerical
              analysis.  A variety of techniques for solving interval linear
              systems of equations are discussed, and these are then tested to
              compare timings and accuracy.  We consider univariate and
              multivariate interval nonlinear systems and describe algorithms
              that enclose all the roots. Finally, we give an application of
              interval analysis.  Interval arithmetic is used to take account of
              rounding errors in the computation of Viswanath's constant, the
              rate at which a random Fibonacci sequence increases.",
  keywords = "interval analysis, interval arithmetic, INTLAB, MATLAB, tutorial,
              interval linear equations, interval nonlinear equations,
              Viswanath's constant"
}

@techreport{Powell:2002:BBP,
  author = "Catherine Powell and David Silvester",
  title = "Black-Box Preconditioning for Mixed Formulation of Self-Adjoint
           Elliptic {PDE}s",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 415",
  month = nov,
  year = 2002,
  note = "Submitted to the conference proceedings of CISC 2002 for publication
          in Springer Lecture Notes in Computer Science and Engineering",
  URL =  "http://www.maths.man.ac.uk/~nareports/narep415.ps.gz",
  abstract = "Mixed finite element approximation of self-adjoint elliptic PDEs
              leads to symmetric indefinite linear systems of equations.
              Preconditioning strategies commonly focus on reduced symmetric
              positive definite systems and require nested iteration. This
              defficiency is avoided if preconditioned \textsc{minres} is
              applied to the full indefinite system.  We outline such a
              preconditioning strategy, the key building block for which is a
              fast solver for a scalar diffusion operator based on black-box
              algebraic multigrid. Numerical results are presented for the
              Stokes equations arising in incompressible flow modelling and a
              variable diffusion equation that arises in modelling potential
              flow. We prove that the eigenvalues of the preconditioned matrices
              are contained in intervals that are bounded independently of the
              discretisation parameter and the PDE coefficients."
}

@techreport{Agyingi:2003:VOP,
  author = "Ephraim Agyingi and Christopher T. H. Baker",
  title = "Variation of Parameters Formulae for {Volterra} Integral Equations",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 414",
  pages = 20,
  month = Nov,
  year = 2003,
  url = "http://www.maths.man.ac.uk/~nareports/narep414.pdf",
  abstract = "In this report, the authors obtain variation of parameters
               formulae for nonlinear Volterra integral equations and
               integro-differential equations. These are related to a variation
               of parameters formula for ordinary differential equations that is
               due to V.M.~Alekseev.  The authors' results are obtained through
               the use of an embedding technique. Some doubts have been cast on
               previous results in the literature. Amongst other things, the
               authors show that a result on which much of the numerical
               analysis for Volterra equations has been founded in the past can
               be obtained by a new approach. "
}

@techreport{Bojanczyk:2002:CIL,
  author = "Adam Bojanczyk and Nicholas J. Higham and Harikrishna Patel",
  title = "The Constrained Indefinite Least Squares Problem: {Theory} and
           Algorithms",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 413",
  pages = "11",
  month = oct,
  year = "2002",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep413.pdf",
  abstract = "We present theory and algorithms for the equality constrained
              indefinite least squares problem, which requires minimization of
              an indefinite quadratic form subject to a linear equality
              constraint. A generalized hyperbolic QR factorization is
              introduced and used in the derivation of perturbation bounds and
              to construct a numerical method. An alternative method is obtained
              by employing a generalized QR factorization in combination with a
              Cholesky factorization. Rounding error analysis is given to show
              that both methods have satisfactory numerical stability properties
              and numerical experiments are given for illustration. This work
              builds on recent work on the unconstrained indefinite least
              squares problem by Chandrasekaran, Gu and Sayed and by the
              present authors.",
   keywords = "equality constrained indefinite least squares problem,
               $J$-orthogonal matrix, hyperbolic rotation, hyperbolic QR
               factorization, generalized hyperbolic QR factorization, rounding
               error analysis, forward stability, perturbation theory, Cholesky
               factorization"
}

@techreport{Ford:2002:IDP,
  author = "Judith M. Ford",
  title = "An Improved {DWT}-Based Preconditioner for Dense Matrix Problems",
  institution = inst-MCCM,
  address =     inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 412",
  month = sep,
  year = 2002,
  URL = "http://www.maths.man.ac.uk/~nareports/narep412.ps.gz",
  abstract = "We present a new preconditioning method for a class of dense
              matrices, based on sparse approximation in a discrete wavelet
              basis. We first prove theoretical results that enable us to reduce
              the cost of an existing wavelet-based preconditioner and then
              incorporate those ideas into the design of our new preconditioner.
              We demonstrate the effectiveness of our methods with numerical
              examples drawn from physical problems including the
              two-dimensional Helmholtz equation for which our preconditioner
              gives savings of about 40\% compared with existing discrete
              wavelet preconditioners and of up to 87\% compared with direct
              solution."
}

@techreport{Higham:2002:CMC,
  author = "Nicholas J. Higham and Matthew I. Smith",
  title = "Computing the Matrix Cosine",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 411",
  pages = 15,
  month = sep,
  year = 2002,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep411.pdf",
  abstract = "An algorithm is developed for computing the matrix cosine,
              building on a proposal of Serbin and Blalock. The algorithm scales
              the matrix by a power of 2 to make the $\infty$-norm less than or
              equal to 1, evaluates a Pad\'e approximant, and then uses the
              double angle formula $\cos(2A) = 2\cos(A)^2-I$ to recover the
              cosine of the original matrix. In addition, argument reduction and
              balancing is used initially to decrease the norm. We give
              truncation and rounding error analyses to show that an [8,8]
              Pad\'e approximant produces the cosine of the scaled matrix
              correct to machine accuracy in IEEE double precision arithmetic,
              and we show that this Pad\'e approximant can be more efficiently
              evaluated than a corresponding Taylor series approximation. We
              also provide error analysis to bound the propagation of errors in
              the double angle recurrence. Numerical experiments show that our
              algorithm is competitive in accuracy with the Schur--Parlett
              method of Davies and Higham, which is designed for general matrix
              functions, and it is substantially less expensive than that method
              for matrices of $\infty$-norm of order 1. The dominant
              computational kernels in the algorithm are matrix multiplication
              and solution of a linear system with multiple right-hand sides, so
              the algorithm is well suited to modern computer architectures.",
  keywords = "matrix function, matrix cosine, matrix exponential, Taylor series,
              Pad\'e approximation, double angle formula, rounding error
              analysis, Schur--Parlett method, MATLAB"
}

@techreport{Higham:2002:MCT,
  author = "Nicholas J. Higham",
  title = "The {Matrix Computation Toolbox} for {MATLAB} (Version 1.0)",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 410",
  pages = "19",
  month = aug,
  year = "2002",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep410.pdf",
  abstract = "The Matrix Computation Toolbox is a collection of MATLAB M-files
              containing functions for constructing test matrices,
              computing matrix factorizations, visualizing matrices,
              and carrying out direct search optimization.",
  keywords = "test matrix, MATLAB, matrix factorizations, Gaussian elimination,
              complete pivoting, rook pivoting, symmetric indefinite matrix,
              Gram-Schmidt orthogonalization, polar decomposition, matrix sign
              function, complete orthogonal decomposition,
              generalized QR factorization, Gauss-Jordan elimination,
              Gershgorin disks, pseudospectrum, field of values, visualization,
              direct search optimization, Strassen's method, vec-permutation
              matrix, Vandermonde matrix"
}

@techreport{Tisseur:2002:TDRSIP,
  author = "Fran\c{c}oise Tisseur",
  title = "Tridiagonal-Diagonal Reduction of Symmetric Indefinite Pairs",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 409",
  pages = "20",
  month = sep,
  year = "2002",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep409.pdf",
  abstract = "We consider the reduction of a symmetric indefinite
              matrix pair $(A,B)$, with $B$ nonsingular, to
              tridiagonal-diagonal form by congruence transformations.
              This is an important reduction in solving polynomial
              eigenvalue problems with symmetric coefficient matrices
              and in frequency response computations.  First the pair
              is reduced to symmetric-diagonal form.  Then we describe
              three methods reducing the symmetric-diagonal pair to
              tridiagonal-diagonal form.  Two of them employ more
              stable versions of Brebner and Grad's pseudosymmetric
              Givens and pseudosymmetric Householder reductions, while
              the third is new and based on a combination of
              Householder reflectors and hyperbolic rotations.  We
              prove an optimality condition for the transformations
              used in the third reduction.  We present numerical
              experiments that compare the different approaches and
              show improvements over Brebner and Grad's reductions.",
  keywords = "symmetric indefinite generalized eigenvalue problem,
              tridiagonalization, hyperbolic rotation, unified
              rotation, hyperbolic Householder reflector."
}

@techreport{Higham:2002:JOM,
  author = "Nicholas J. Higham",
  title = "{$J$}-Orthogonal Matrices: {Properties} and Generation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 408",
  pages = "15",
  month = sep,
  year = 2002,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep408.pdf",
  abstract = "A real, square matrix $Q$ is $J$-orthogonal if $Q^TJQ = J$, where
              the signature matrix $J = \diag(\pm 1)$. $J$-orthogonal matrices
              arise in the analysis and numerical solution of various matrix
              problems involving indefinite inner products, including, in
              particular, the downdating of Cholesky factorizations. We present
              techniques and tools useful in the analysis, application and
              construction of these matrices, giving a self-contained treatment
              that provides new insights. First, we define and explore the
              properties of the exchange operator, which maps $J$-orthogonal
              matrices to orthogonal matrices and vice-versa. Then we show how
              the exchange operator can be used to obtain a hyperbolic CS
              decomposition of a $J$-orthogonal matrix directly from the usual
              CS decomposition of an orthogonal matrix. We employ the
              decomposition to derive an algorithm for constructing random
              $J$-orthogonal matrices with specified norm and condition number.
              We also give a short proof of the fact that $J$-orthogonal
              matrices are optimally scaled under two-sided diagonal scalings.
              We introduce the indefinite polar decomposition and investigate
              two iterations for computing the $J$-orthogonal polar factor: a
              Newton iteration involving only matrix inversion and a Schulz
              iteration involving only matrix multiplication. We show that these
              iterations can be used to $J$-orthogonalize a matrix that is not
              too far from being $J$-orthogonal.",
  keywords = "$J$-orthogonal matrix, exchange operator, gyration operator,
               sweep operator, principal pivot transform, hyperbolic CS
               decomposition, two-sided scaling, indefinite least squares
               problem, hyperbolic QR factorization, Newton's method, Schulz
               iteration"
}

@techreport{Tisseur:2002:STTSM,
  author = "Seamus. D. Garvey and Fran\c{c}oise Tisseur
            and Mike I. Friswell and John E. T. Penny and
            Uwe Prells",
  title = "Simultaneous Tridiagonalization of Two Symmetric Matrices",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 407",
  pages = "18",
  month = aug,
  year = "2002",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep407.pdf",
  abstract = "We show how to simultaneously reduce a pair of symmetric
              matrices to tridiagonal form by congruence
              transformations.  No assumptions are made on the
              nonsingularity or definiteness of the two matrices.  The
              reduction follows a strategy similar to the one used for
              the tridiagonalization of a single symmetric matrix via
              Householder reflectors.  Two algorithms are proposed,
              one using non-orthogonal rank-one modifications of the
              identity matrix and the other, more costly but more
              stable, using a combination of Householder reflectors
              and non-orthogonal rank-one modifications of the
              identity matrix with minimal condition numbers.  Each of
              these tridiagonalization processes requires $O(n^3)$
              arithmetic operations and respects the symmetry of the
              problem.  We illustrate and compare the two algorithms
              with some numerical experiments.",
keywords = "Symmetric matrices,
            generalized eigenvalue problem,
            tridiagonalization,
            symmetric quadratic eigenvalue problem."
}

@techreport{Mihajlovic:2002:EPS,
  author = "Milan Mihajlovic and David Silvester",
  title = "Efficient Parallel Solvers for the Biharmonic Equation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 406",
  month = jul,
  year = 2002,
  note = "Submitted to Parallel Computing",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep406.ps.gz",
  abstract = "We examine the convergence characteristics and performance of
       parallelised Krylov subspace solvers applied to the linear algebraic
       systems that arise from low-order mixed finite element approximation of
       the biharmonic problem. Our strategy results in preconditioned systems
       that have nearly optimal eigenvalue distribution, which consists of a
       tightly clustered set together with a small number of outliers. We
       implement the preconditioner operator in a ``black-box'' fashion using
       publicly available parallelised sparse direct solvers and multigrid
       solvers for the discrete Dirichlet Laplacian. We present convergence and
       timing results that demonstrate the efficiency and scalability of our
       strategy when implemented on contemporary computer architectures."
}

@techreport{Davies:2002:USV,
  author = "Philip I. Davies and Matthew I. Smith",
  title = "Updating the Singular Value Decomposition",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 405",
  pages = "30",
  month = aug,
  year = "2002",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep405.ps.gz",
  abstract = "The spectral decomposition of a symmetric matrix $A$
              with small off-diagonal and distinct diagonal elements
              can be approximated using a direct scheme of R. Davies and Modi
              [\emph{Linear Algebra and Appl.,} 77:61-74, 1986].
              In this paper a generalization of this method
              for computing the singular value decomposition of
              close-to-diagonal $A\in\R^{m \times n}$ is presented.
              When $A$ has repeated or ``close'' singular values it is
              possible to apply the direct method to split the problem
              in two with one part containing the well-separated singular
              values and one requiring the computation of the ``close''
              singular values."
}

@techreport{Davies:2002:SPA,
  author = "Philip I. Davies and Nicholas J. Higham",
  title = "A {Schur--Parlett} Algorithm for Computing Matrix Functions",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 404",
  pages = "22",
  month = jul,
  year = "2002",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep404.pdf",
  abstract = "An algorithm for computing matrix functions is presented. It
              employs a Schur decomposition with reordering and blocking
              followed by the block form of a recurrence of Parlett, with
              functions of the nontrivial diagonal blocks evaluated via a Taylor
              series. A parameter is used to balance the conflicting
              requirements of producing small diagonal blocks and keeping the
              separations of the blocks large. The algorithm is intended
              primarily for functions having a Taylor series with an infinite
              radius of convergence, but it can be adapted for certain other
              functions, such as the logarithm. Novel features introduced here
              include a convergence test that avoids premature termination of
              the Taylor series evaluation and an algorithm for reordering and
              blocking the Schur form. Numerical experiments show that the
              algorithm is competitive with existing special-purpose algorithms
              for the matrix exponential, logarithm, and cosine. Nevertheless,
              the algorithm can be numerically unstable with the default choice
              of its blocking parameter (or in certain cases for all choices),
              and we explain why determining the optimal parameter appears to be
              a very difficult problem. A MATLAB implementation is available
              that is much more reliable than the existing MATLAB function
              \texttt{funm}."
  }

@techreport{Song:PTD:2002,
  author = "Y. Song and C. T. H. Baker",
  title = "Perturbation Theory for Discrete {Volterra} Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 403",
  pages = "18",
  month = jun,
  year = "2002",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep403.ps.gz",
  abstract = "Fixed  point theory is used to investigate nonlinear discrete
          Volterra equations that are perturbed versions of linear equations.
          Sufficient conditions are established ($i$) to ensure that stability
          (in a sense that is defined) of the solutions of the linear equation
          implies a corresponding stability of the zero solution of the
          nonlinear equation and ($ii$) to ensure the existence of
          asymptotically periodic solutions",
    keywords = "Fixed point theorem, perturbation, discrete Volterra equations,
               periodic, asymptotically periodic."
}

@techreport{Song:SDV:2002,
  author = "Y. Song and C. T. H. Baker",
  title = "Stability in Discrete {Volterra} Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 402",
  pages = "19",
  month = jun,
  year = "2002",
  URL =         "ftp://ftp.ma.man.ac.uk/pub/narep/narep402.ps.gz",
  abstract = "A linearized stability theory for discrete Volterra
        equations of the type \[
        x(n)=f(n)+\sum^n_{j=0}K\big(n, j, x(n)\big),  \mbox{($n\ge 0$)}
        \] is developed: Various types of stability are defined and
        studied by using the fundamental and resolvent matrices that
        are associated with linear equations of the form \[ x(n)=f(n)
        + \sum^n_{j=0}B(n, j)x(j),  \mbox{($n\ge 0$)},  \] and several
        necessary and sufficient conditions for stability are obtained
        for solutions of this linear equation. In conclusion,  a result
        concerning the stability of a solution of the nonlinear
        equation is provided under certain assumptions.
        The results parallel corresponding results for classical
        Volterra integral equations of the second kind",
  keywords = "Stability,  asymptotic stability,  discrete Volterra equations,
              resolvent matrix,  fundamental matrix."
}

@techreport{MCCM:2002:AR1,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 2001",
  type =         t-NAREP,
  number =       "No. 401",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "21",
  month =        may,
  year =         "2002",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep401.pdf",
  }

@techreport{Baker:2002:PSD,
  author = "C. T. H. Baker and Y. Song",
  title = "Periodic Solutions of Discrete {Volterra} Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 400",
  pages = "21",
  month = may,
  year = "2002",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep400.ps.gz",
  abstract =
        "In this paper we investigate periodic solutions of linear
         and nonlinear discrete Volterra equations of convolution,
         or non-convolution type with unbounded memory. For linear
         discrete Volterra equations,  we establish  Fredholm's
         alternative theorem for convolution equations and we show
         that, under certain conditions, only one particular bounded
         initial function can generate a  unique periodic solution
         and this unique periodic solution attracts all other
         solutions with bounded initial function.  In other words,
         we show that all solutions of linear discrete Volterra
         equations are asymptotically periodic under certain conditions.",
   keywords =
         "Periodic, asymptotically periodic, discrete Volterra
          equations, resolvent matrices, Fredholm's alternative."
}

@techreport{Powell:2002:OPR,
  author = "Catherine Powell and David Silvester",
  title = "Optimal Preconditioning for {Raviart-Thomas} Mixed Formulation
           of Second-Order Elliptic Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 399",
  month = mar,
  year = 2002,
  note = "Submitted to the SIAM Journal on Matrix Analysis and Applications",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep399.ps.gz",
  abstract = "We evaluate two preconditioning strategies for the indefinite
              linear system obtained from Raviart-Thomas mixed finite element
              formulation of a second-order elliptic problem with variable
              coefficients. In contrast to other approaches, our emphasis is on
              linear algebra; we establish the optimality of the
              preconditioning using basic properties of the finite element
              matrices. The underlying saddle-point problem is well-posed in
              two function spaces, $H(div)\times L^{2}$ and $L^{2} \times
              H^{1}$, leading to the possibility of two distinct types of
              preconditioner. For homogeneous Dirichlet boundary conditions,
              the discrete problems are identical. This motivates the use of
              Raviart-Thomas approximation in both frameworks yielding a
              non-conforming method in the second case. We construct two
              block-diagonal preconditioners and establish a new uniform
              eigenvalue bound. Trials of preconditioned MINRES illustrate that
              both preconditioners are robust with respect to the mesh
              parameter, and with respect to jumps and anisotropies in the
              diffusion coefficients.  "
}

@techreport{Baker:2002:DVO,
  author = "C. T. H. Baker and Y. Song",
  title = "Discrete Volterra Operators, Fixed Point Theorems \& their
           Application",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 398",
  pages = "23",
  month = mar,
  year = "2002",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep398.ps.gz",
  abstract =

        "The existence of solutions of nonlinear discrete Volterra
        equations is established.
           We define discrete Volterra operators on normed spaces of
        infinite sequences of finite-dimensional vectors.  and present
        some of their basic properties (continuity, boundedness, and
        representation). The treatment relies upon the use of
        coordinate functions, and the existence results are obtained
        using fixed point theorems for discrete Volterra operators on
        infinite-dimensional spaces based on fixed point theorems of
        Schauder, Rothe, and Altman, and Banach's contraction mapping
        theorem, for finite-dimensional spaces.",
   keywords =
        "Fixed point theorem, discrete Volterra equations, discrete
        Volterra operators."
}

@techreport{Bojanczyk:2002:SIL,
  author = "Adam Bojanczyk and Nicholas J. Higham and Harikrishna Patel",
  title = "Solving the Indefinite Least Squares Problem by Hyperbolic {QR}
           Factorization",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 397",
  pages = "19",
  month = jan,
  year = "2002",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep397.pdf",
  abstract = "The indefinite least squares (ILS) problem is to minimize a
              certain type of indefinite quadratic form. We develop
              perturbation theory for the problem and identify a condition
              number. We describe and analyze a method for solving the ILS
              problem based on hyperbolic QR factorization. This method is
              faster than one recently proposed by Chandrasekaran, Gu and Sayed
              that employs both QR and Cholesky factorizations. We give a
              rounding error analysis of the new method and use the
              perturbation theory to show that under a reasonable assumption
              the method is forward stable. Our analysis is quite general and
              sheds some light on the stability properties of hyperbolic
              transformations. Our numerical experiments suggest that the new
              method is in practice just as accurate as the method of
              Chandrasekaran, Gu and Sayed.",
   keywords = "indefinite least squares problem, downdating, hyperbolic
               rotation, hyperbolic QR factorization, rounding error analysis,
               forward stability, perturbation theory, condition number"
}

@techreport{Baker:2002:BSD,
  author = "C. T. H. Baker and J. M. Ford and N. J. Ford",
  title = "Bifurcations in Stochastic Delay Differential Equations---{I}:
           {Numerical} Investigation ",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 396",
  pages = "18",
  month = mar,
  year = "2002",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep396.ps.gz",
  abstract =
        "We consider {\em stochastic delay differential equations} of
        the form \[ dY(t)=F(t,Y(t),Y(t-\tau))dt + G(t,Y(t),Y(t-\tau))
        dW(t) \quad (t \geq t_0), \] interpreted in the It\^{o} sense
        where $W(t)$ is a standard {\em Wiener process}, with $ Y(t) =
        \Phi(t)$ for $t \in [-\tau,0]$ (where $\tau > 0$ is the
        constant lag, or `time-lag').  We are interested in {\em
        bifurcations} (that is, changes in the qualitative behaviour
        of the solutions to these equations) and we draw on insights
        from the deterministic counterpart \[
        y'(t)=f(t,y(t),y(t-\tau)), \quad \mbox{ for $t \geq t_0$}, \]
        with $y(t) = \phi(t)$ for $ t \in [t_0-\tau,t_0]$ and constant
        lag $\tau$.  For these deterministic problems we can exploit
        known theory and numerical results that enable us to predict
        where changes occur in the behaviour of the solutions (exact
        and approximate) of the equation.
        For those unfamiliar with the area, rather diverse components
        of the mathematical background are necessary to understand the
        questions of interest.  In this introductory report we first
        review some deterministic results and some basic elements of
        the stochastic analysis that ($i$) suggests lines of
        investigation for the stochastic case and ($ii$) are expected
        to facilitate the theoretical investigation of the stochastic
        problem. We then give numerical evidence that the presence of
        noise, corresponding to the change from the deterministic to
        the stochastic form, may change the qualitative behaviour of
        the solutions in an interesting way.",
   keywords =
        "Bifurcation and stability, stochastic and deterministic delay
        equations with fixed lag, Euler-type methods with fixed step,
        numerical experiments."
}

@techreport{Tisseur:2001:CBE,
  author = "Fran\c{c}oise Tisseur",
  title = "A Chart of Backward Errors and Condition Numbers for
           Singly and Doubly Structured Eigenvalue Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 395",
  pages = "23",
  month = dec,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep395.ps.gz",
  abstract = "We present a chart of structured backward errors for
              approximate eigenpairs of singly and doubly structured
              eigenvalue problems.  We aim to give, wherever possible,
              formulae that are inexpensive to compute so that they
              can be used routinely in practice.  We identify a number
              of problems for which the structured backward error is
              within a factor $\sqrt{2}$ of the unstructured backward
              error. This paper collects, unifies and extends existing
              work on this subject.",
}

@techreport{Tisseur:2001:IGT,
  author = "Jean-Pierre Dedieu and Myong-Hi Kim and Michael Shub
            and Fran\c{c}oise Tisseur",
  title = "Implicit Gamma Theorems {(I)}: {Pseudoroots} and Pseudospectra",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 394",
  pages = "25",
  month = Nov,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep394.ps.gz",
  abstract = "Let $g: \E \rightarrow \F $ be an analytic function
              between two Hilbert spaces $\E$ and $\F$. We study the
              set $g(B(x,\eps ))\subset \F$, the image under $g$ of
              the closed ball about $x\in\E$ with radius $\eps$. When
              $g(x)$ expresses the solution of an equation depending
              on $x$ then the elements of $g(B(x,\eps ))$ are
              $\eps$-pseudosolutions . Our aim is to investigate the
              size of the set $g(B(x,\eps ))$. We derive upper and
              lower bounds of the following form $$ g(x) + Dg(x)\left
              ( B(0, c_1 \eps ) \right ) \subseteq g(B(x,\eps ))
              \subseteq g(x) + Dg(x)\left ( B(0, c_2 \eps ) \right ),
              $$ where $Dg(x)$ denotes the derivative of $g$ at
              $x$. We consider both the case where $g$ is given
              explicitly and the case where $g$ is given
              implicitly. We apply our results to the implicit
              function associated with the evaluation map, namely the
              solution map, and to the polynomial eigenvalue
              problem. Our results are stated in terms of an invariant
              $\gamma$ which has been extensively used by various
              authors in the study of Newton's method. The main tool
              used here is an implicit $\gamma$ theorem, which
              estimates the $\gamma$ of an implicit function in terms
              of the $\gamma$ of the function defining it.",
}

@techreport{Cheng:2001:ILBA,
  author = "Sheung Hun Cheng and Nicholas J. Higham",
  title = "Implementation for {LAPACK} of a Block Algorithm for Matrix 1-Norm
           Estimation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 393",
  pages = "19",
  month = aug,
  year = "2001",
  note = "LAPACK Working Note 152",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep393.ps.gz",
  abstract = "We describe double precision and complex*16 Fortran~77
              implementations, in LAPACK style, of a block matrix 1-norm
              estimator of Higham and Tisseur. This estimator differs from that
              underlying the existing LAPACK code, \texttt{xLACON}, in that it
              iterates with a matrix with $t$ columns, where $t\ge 1$ is a
              parameter, rather than with a vector, and so the basic
              computational kernel is level 3 BLAS operations. Our experiments
              with random matrices on a Sun SPARCStation Ultra-5 show that with
              $t=2$ or $4$ the new code offers better estimates than
              \texttt{xLACON} with a similar execution time. Moreover, with
              $t>2$, estimates exact over $95\%$ and $75\%$ of the time are
              achieved for the real and complex version respectively, with
              execution time growing much slower than $t$. We recommend this
              new code be included as an auxiliary routine in LAPACK to
              complement the existing LAPACK routine \texttt{xLACON}, upon
              which the various drivers should still be based for compatibility
              reasons."
        }

@techreport{Smith:2001:ASAC,
  author =       "Matthew I. Smith",
  title =        "A {Schur} Algorithm for Computing Matrix $p$th Roots",
  institution =   inst-MCCM,
  address =       inst-MCCM:adr,
  type =          t-NAREP,
  number =       "No. 392",
  month =         jul,
  year =         "2001",
  pages =        "18",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep392.ps.gz",
  abstract =     "Any nonsingular matrix has $p$th roots,
             and they can be classified into two groups, those that are
             functions of $A$ and those that are not. Matrix $p$th roots can be
             computed by a specialized version of Newton's method, but we show
             that this iteration has poor convergence and stability properties
             in general. We present a Schur algorithm for computing a matrix
             $p$th root that generalizes methods of Bj{\"{o}}rck and Hammarling
             [\emph{Linear Algebra and Appl.,} 52/53:127-140, 1983] and Higham
             [\emph{Linear Algebra and Appl.,} 88/89:405-430, 1987] for the
             square root. Our algorithm forms a Schur decomposition of $A$ and
             computes a $p$th root of the (quasi)triangular factor by a
             recursion. The backward error associated with our Schur method is
             examined, and we find that the method has excellent numerical
             stability.",
   keywords = "matrix $p$th root, Schur algorithm, Newton's method,
               commutativity, stability, real arithmetic, rounding error
               analysis"
}

@techreport{MCCM:2001:AR0,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 2000",
  type = t-NAREP,
  number =       "No. 391",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "23",
  month =        mar,
  year =         "2001",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep391.pdf",
  }

@techreport{FordJ:2001:PFWT,
  author = "J. M. Ford and K. Chen and N. J. Ford",
  title = "Parallel Implementation of Fast Wavelet Transforms",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 389",
  pages = "15",
  month = feb,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep389.ps.gz",
  abstract = "We look at the use of parallel techniques for fast wavelet transforms, focusing
	particularly on how to make best use of any given number of available  processors.
	Efficient algorithms for parallel  FWT have been established in the case where the
	lengths of the vectors to be transformed and the number of processors are both powers of 2.
	We show that a na{\"}{\i}ve extension to arbitrary vector length and number of processors
	may give poor performance. We present new algorithms for the distribution of matrix
	elements between processors that overcome these difficulties. We have implemented our
	algorithms on a multi-processor Sun Sparc machine, and we present analytical and
	experimental results demonstrating significant improvements in performance of FWT
	calculation.",
  keywords = "parallel, FWT, implementation, MPI"
        }

@techreport{FordJ:2001:SPFW,
  author = "J. M. Ford and K. Chen and N. J. Ford",
  title = "Small-Scale Parallel Implementation of Fast Wavelet Transforms",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 388",
  pages = "7",
  month = feb,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep388.ps.gz",
  abstract = "We look at the use of parallel techniques for Fast Wavelet Transforms, focusing
	particularly on how to make best use of a very limited (and, perhaps, variable)
	number of available  processors. Efficient algorithms for parallel  FWT have been
	established in the case where the lengths of the vectors to be transformed and the
	number of processors are both powers of 2. We show that a na{\"}{\i}ve extension to
	arbitrary vector length and number of processors may give poor performance or, in
	some cases, be impossible. We present new algorithms for the distribution of matrix
	elements between processors that overcome these difficulties. We have implemented our
	algorithms on a local area network of inexpensive PCs using MPI, and present
	experimental results demonstrating significant improvements in performance of FWT
	calculation.",
  keywords = "PC, parallel, FWT, MPI"
        }

@techreport{Ford:2001:NAFDE,
  author = "N. J. Ford and A. C. Simpson",
  title = "Numerical and Analytical Treatment of Differential Equations of
           Fractional Order",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 387",
  pages = "7",
  month = feb,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep387.ps.gz",
  abstract = "We begin with a brief survey of differential equations of fractional order and
	their applications. We focus on key issues in the exact and numerical solution
	of the equations. We summarise our recent results: We describe the importance
	of the appropriate choice of	initial conditions. We show how the behaviour of
	solutions changes with small changes in the fractional order. We discuss the use
	of {\em systems} of fractional order equations as a tool to solve equations of
	higher order. We give results that show how the computational effort required to
	solve an equation of fractional order can be reduced, either through careful
	application of the  {\em short memory}  principle or through the use of schemes
	of varying step-length.",
  keywords = "short memory, systems, fractional order"
        }

@techreport{Ford:2001:AFDG1,
  author = "N. J. Ford and A. C. Simpson",
  title = "The Approximate Solution of Fractional Differential Equations of Order
 	   Greater Than 1",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 386",
  pages = "6",
  month = feb,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep386.ps.gz",
  abstract = "Differential equations of fractional order arise, for example, in the generalised
	viscoelastic constitutive equations, the generalised Basset force, and
        in fractional Brownian motion. Numerical methods for solution of
        equations of fixed order $\alpha$ are well established, particularly
        for $0<\alpha<1$, however many of the analytical results depend on the
        assumption that the equation has a particularly simple form. Recent
        work (by Diethelm and Ford and others) has considered approaches for
        solving more complicated and multi-term equations. In this paper we
        present an alternative approach to numerical solution of order greater
        than one and we conclude the paper by comparing the numerical results
        for our proposed method with known numerical and analytical solutions
        for prototype nonlinear problems.",
  keywords = "systems, fractional order"
        }

@techreport{Ford:2001:SVA,
  author = "N. J. Ford and A. C. Simpson",
  title = "The Numerical Solution of Fractional Differential Equations: Speed
		  Versus Accuracy",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 385",
  pages = "12",
  month = feb,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep385.ps.gz",
  abstract = "This paper is concerned with the development of efficient algorithms for the
              approximate solution of fractional differential equations of the
              form $D^\alpha y(t)=f(t,y(t)), \mathbb{R}^+-\mathbb{N}$. We
              briefly review standard numerical techniques for the solution of
              ($\dagger$) and we consider how the computational cost may be
              reduced by taking into account the structure of the calculations
              to be undertaken. We analyse the {\it fixed memory principle} and
	present an alternative {\it nested mesh} variant that gives a good
              approximation to the true solution at reasonable computational
              cost. We conclude with some numerical examples.",
   keywords = "short memory, finite memory, logarithmic memory"
   }

@techreport{Edwards:2001:BSDE,
  author = "J. T. Edwards and N. J. Ford",
  title = "Boundedness and Stability of Difference Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 384",
  pages = "14",
  month = feb,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep384.ps.gz",
  abstract = "This paper is concerned with the qualitative behaviour of solutions to
		difference equations.	We focus on boundedness and stability of solutions
		and we present a unified theory that applies both to autonomous and
		non-autonomous equations and to nonlinear equations as well as linear
		equations. Our presentation brings together new, established, and
		hard-to-find results from the literature and provides a theory that is both
		memorable and easy to apply. We show how the theoretical results given here
		relate to some of those in the established literature and by means of simple
		examples we indicate how the use of Lipschitz constants in this way can provide
		useful insights into the qualitative behaviour of solutions to some nonlinear
		problems including those arising in numerical analysis.",
  keywords = "boundedness, stability, Lipschitz"
        }

@techreport{Ford:2001:SHQ,
  author = "N. J. Ford and P. D. Crofts and R. H. T. Edwards",
  title = "Sensitivity of Hospital Clinic Queues to Patient Non-Attendance",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "383",
  month = feb,
  pages = "6",
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep383.ps.gz",
  abstract = "The recently introduced Patients' Charter offers patients in the United
	Kingdom a guaranteed minimum quality of service when they require hospital
	treatment. For example, the following guarantees are offered: 1. No
	patient should be required to spend more than one year	awaiting an offer
	of an appointment with a Consultant 2. Patients should not be expected to
	wait in clinic for more than 30 minutes without seeing a doctor. Health
	economists are active in devising strategies and systems to deal with these
	enhanced patient expectations. In this paper, we describe the work of an
	interdisciplinary team, funded by the UK National Health Service Research
	and Development which investigated computer simulation of clinic queues to
	provide further insight into the problem and to produce computer software
	that can assist in clinic optimisation.",
   keywords = "optimisation, modelling, queues"
       }

@techreport{Ford:2001:NDSS,
  author = "Neville J Ford and Sjoerd M Verduyn Lunel",
  title = "Numerical approximation of delay differential equations with small
		  solutions",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "382",
  month = feb,
  pages = "6",
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep382.ps.gz",
  abstract = "We are concerned with the detection of small solutions of delay differential equations.
	We consider, as prototype, an equation of the form $\dot{y}(t)=b(t)y(t-\tau)$ for
	which conditions on the function $b$ that govern the existence (or otherwise) of small
	solutions are well-known. We apply simple numerical methods to the problem and
	investigate the characteristic values of the solution operator. We show how these
	characteristic roots may be used to detect the existence of small solutions.",
  keywords = "small solutions, floquet theory, characteristic values"
        }

@techreport{Ford:2001:CSS,
  author = "Neville J. Ford and Sjoerd M. Verduyn Lunel",
  title = "Characterising Small Solutions in Delay Differential Equations
	   Through Numerical Approximations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "381",
  month = feb,
  pages = "15",
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep381.ps.gz",
  abstract = "The existence of small solutions, i.e., solutions that decay
              faster than any exponential, for linear time-dependent
              differential delay equations with bounded coefficients depends on
              specific properties of the coefficients. Although small solutions
              do not occur in the finite dimensional approximations of the
              differential delay equation we show that the existence of small
              solutions for differential delay equations can be predicted from
              the behaviour of the spectrum of the finite dimensional
              approximations.",
  keywords = "small solutions, floquet theory"
      }

@techreport{Edwards:2001:IDBP,
  author = "J. T. Edwards, J. A. Roberts and N. J. Ford",
  title = "The numerical Solution of an Integro-Differential Equation Close to
	   Bifurcation Points",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 380",
  pages = "6",
  month = feb,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep380.ps.gz",
  abstract = "In this paper we  investigate the qualitative behaviour of numerical
	approximations to a Volterra integro-differential equation. We consider (as
	prototype) a linear problem with fading memory kernel of the form
	$$y'(t)=-\int^t_0e^{-\lambda(t-s)}y(s)ds, \qquad \qquad y(0)=1$$ and we consider
	the performance of simple numerical schemes applied to solve the equation. We
	are concerned with the preservation (or otherwise) of qualitative properties of
	the analytical solution in the numerical approximation.   We outline the known
	stability behaviour and derive the values of $\lambda$ at which the true solution
	bifurcates. We give the corresponding analysis for the discrete schemes and we show
	that, as the step size of the numerical scheme decreases, the bifurcation points
	tend towards those of the continuous scheme. We illustrate our results with some
	numerical examples.",
  keywords = "bifurcations, numerical approximations, integro-differential equations"
}

@techreport{Ford:2001:FESO,
  author = "N. J. Ford and K. Diethelm",
  title = "The numerical Solution of Linear and Non-Linear Fractional
 	  Differential Equations Involving Fractional Derivatives of
		  several orders",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 379",
  month = feb,
  pages = 14,
  year = 2001,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep379.ps.gz",
  abstract = "We consider fractional differential equations of the very general
	form $y^{(\alpha)}(t)= f(t,y(t),y^{(\beta_1)}(t),y^{(\beta_2)}(t),
	\ldots,y^{(\beta_n)}(t))$ with $\alpha>\beta_1>\beta_2>\ldots>\beta_n$
	and $\alpha-\beta_1\le1$,$\beta_j-\beta_{j+1}\le1$, $0\le\beta_n\le1$,
	and where the derivatives are understood in Caputo's sense.For equations of
	this type, we discuss the questions of existence and uniqueness of solutions,
	and we investigate how the solutions depend on the given data. We propose
	convergent and stable numerical methods based on a nearly equivalent system
	of fractional differential equations of order not exceeding $\beta_n$. We
	give particular emphasis to the practically important linear case and also
	present results for nonlinear problems.",
  keywords = "milti-term equation, fractional order, fractional Adams method"
}

@techreport{Ford:2001:NBT,
  author = "N. J. Ford and K. Diethelm",
  title = "Numerical Solution of the {Bagley Torvik} Equation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 378",
  month = feb,
  pages = 12,
  year = 2001,
  note = "Submitted to BIT",
  URL = "http://www.maths.man.ac.uk/~nareports/narep378.ps.gz",
  abstract = "We consider the numerical solution of the Bagley-Torvik
	equation $$ A y''(t) + B D_*^{3/2}y(t) + C y(t) = f(t)$$
	as a prototype fractional differential equation with two
	derivatives. Approximate solutions have recently been proposed in
	the book and papers of Podlubny in which the solution obtained with
	approximate methods is compared to the exact solution. In this paper
	we consider the reformulation of the Bagley-Torvik equation as a
	system of fractional differential equations of order $1/2$. This
	allows us to propose numerical methods for its solution which are
	consistent and stable and have arbitrarily high order.",
  keywords = "Bagley Torvik, multi-term, fractional order"
}

@techreport{Ford:2001:AFDE,
  author = "N. J. Ford and K. Diethelm",
  title = "Analysis of Fractional Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 377",
  month = feb,
  pages = 17,
  year = 2001,
  note = "To appear in Journal of Mathematical Analysis and Applications",
  URL = "http://www.maths.man.ac.uk/~nareports/narep377.ps.gz",
  abstract = "We discuss existence, uniqueness and {\it structural stability}
	of solutions of nonlinear differential equations of fractional
	order. The differential operators are taken in the
	Riemann-Liouville sense and the initial conditions are specified
	according to Caputo's suggestion, thus allowing for interpretation
	in a physically meaningful way. We investigate in particular the
	dependence of the solution on the order of the differential
	equation and on the initial condition and we relate our results
	to the selection of appropriate numerical schemes for the solution
	of fractional differential equations.",
  keywords = "fractional order, stability, existence, uniqueness, numerical schemes"
}

@techreport{Bolton:2001:MCLS,
  author = "P. Bolton and J. Stratakis and R. W. Thatcher",
  title = "Mass Conservation in Least Squares Methods for {Stokes} Flow",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 376",
  pages = "22",
  month = jun,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep376.pdf"
        }

@techreport{Basebi:2001:ASMM,
  author =      "Thebe Basebi and Ruth M. Thomas",
  title =       "A Study of Moving Mesh Methods Applied to a Thin Flame
                 Propagating in a Detonator Delay Element",
  type = t-NAREP,
  number =      "No. 375",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        oct,
  year =        "2001",
  URL =         "http://www.maths.man.ac.uk/~nareports/narep375.ps.gz",
  abstract =    "We study the application of moving mesh methods to a
                 one-dimensional (time dependent) detonator delay element
                 problem. We consider moving mesh methods based on the
                 equidistribution principle derived by
                 Huang $et$ $al.$ (1994). Adaptive mesh methods have been widely
                 used recently to solve time dependent partial differential
                 equations having large solution gradients. Significant
                 improvements in accuracy and efficiency are achieved by
                 adapting the nodes (mesh points) so that they are
                 concentrated about areas of large solution variations.
                 Each system of equations for the moving mesh methods is
                 solved in conjunction with the detonator problem. In this
                 paper, the system of ordinary differential equations that
                 results (after discretising in space) is solved using the
                 double precision version of the stiff ordinary differential
                 equation solver DASSL. The numerical results clearly
                 demonstrate that the moving mesh methods are capable of
                 tracking the deflagration wave as it travels down the
                 detonator delay element more accurately and more efficiently
                 than a fixed mesh method.",
  keywords =    "moving mesh methods, detonator delay element, time dependent
                 partial differential equations, scaling, monitor functions."
}

@techreport{Cheng:2001:PIB,
  author = "Sheung Hun Cheng and Nicholas J. Higham",
  title = "Parallel Implementation of a Block Algorithm for Matrix 1-Norm
           Estimation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "No. 374",
  pages = "9",
  month = feb,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep374.pdf",
  abstract = "We describe a parallel Fortran~77 implementation, in ScaLAPACK
              style, of a block matrix 1-norm estimator of Higham and Tisseur.
              This estimator differs from that underlying the existing
              ScaLAPACK code, \texttt{PxLACON}, in that it iterates with a
              matrix with $t$ columns, where $t\ge 1$ is a parameter, rather
              than with a vector, and so the basic computational kernel is
              level 3 BLAS operations. Our experiments on an SGI Origin2000
              show that with $t=2$ or $4$ the new code offers better estimates
              than \texttt{PDLACON} with a similar execution time. Moreover,
              with $t>4$, estimates exact over $90\%$ of the time are achieved
              with execution time growing much slower than $t$.",
  keywords = ""
        }

@techreport{Higham:2001:DDH,
  author = "Nicholas J. Higham and Fran\c{c}oise Tisseur and
            Paul M. Van Dooren",
  title = "Detecting a Definite {Hermitian} Pair
           and a Hyperbolic or Elliptic Quadratic Eigenvalue Problem,
           and Associated Nearness Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "373",
  month = feb,
  pages = "18",
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep373.ps.gz",
  abstract = "An important class of generalized eigenvalue problems $Ax =
              \lambda Bx$ is those in which $A$ and $B$ are Hermitian and some
              real linear combination of them is definite. For the quadratic
              eigenvalue problem (QEP) $(\l^2 A + \l B + C)x = 0$ with
              Hermitian $A$, $B$ and $C$ and positive definite $A$, particular
              interest focuses on problems in which $(x^*Bx)^2 - 4
              (x^*Ax)(x^*Cx)$ is one-signed for all nonzero $x$---for the
              positive sign these problems are called hyperbolic and for the
              negative sign elliptic. The important class of overdamped
              problems arising in mechanics is a sub-class of the hyperbolic
              problems. For each of these classes of generalized and quadratic
              eigenvalue problems we show how to check that a putative member
              has the required properties and we derive the distance to the
              nearest problem outside the class. For definite pairs $(A,B)$ the
              distance is the Crawford number, and we derive bisection and
              level set algorithms both for testing its positivity and for
              computing it. Testing hyperbolicity of a QEP is shown to reduce
              to testing a related pair for definiteness. The distance to the
              nearest non-hyperbolic or non-elliptic $n\times n$ QEP is shown
              to be the solution of a global minimization problem with $n-1$
              degrees of freedom. Numerical results are given to illustrate the
              theory and algorithms.",
  keywords = "Generalized eigenvalue problem, definite pair, Crawford number,
              quadratic eigenvalue problem, hyperbolic system, elliptic system,
              overdamped system, gyroscopic system, bisection algorithm, level
              set algorithm, nearness problems"
        }

@techreport{Higham:2001:MPP,
  author = "Nicholas J. Higham and Fran\c{c}oise Tisseur",
  title = "More on Pseudospectra for Polynomial Eigenvalue Problems
           and Applications in Control Theory",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "372",
  month = jan,
  pages = "21",
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep372.ps.gz",
  abstract = "Definitions and characterizations of pseudospectra are given for
              rectangular matrix polynomials expressed in homogeneous form:
              $P(\a,\b) = (\a^dA_d + \a^{d-1}\b A_{d-1} + \cdots + \b^d A_0)$.
              It is shown that problems with infinite (pseudo)eigenvalues are
              elegantly treated in this framework. For such problems
              stereographic projection onto the Riemann sphere is shown to
              provide a convenient way to visualize pseudospectra. Lower bounds
              for the distance to the nearest non-regular polynomial and the
              nearest uncontrollable $d$th order system (with equality for
              standard state-space systems) are obtained in terms of
              pseudospectra, showing that pseudospectra are a fundamental tool
              for reasoning about matrix polynomials in areas such as control
              theory. How and why to incorporate linear structure into
              pseudospectra is also discussed by example.",
  keywords = "polynomial eigenvalue problem, $\lambda$-matrix, matrix
              polynomial, homogeneous form, pseudospectrum, stereographic
              projection, Riemann sphere, nearest non-regular polynomial,
              nearest uncontrollable system, structured perturbations"
        }

@techreport{Tisseur:2001:BEM,
  author = "Nicholas J. Higham and Fran\c{c}oise Tisseur",
  title = "Bounds for Eigenvalues of Matrix Polynomials",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = t-NAREP,
  number = "371",
  month = jan,
  pages = "16",
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep371.pdf",
  abstract = "Upper and lower bounds are derived for the absolute values of the
              eigenvalues of a matrix polynomial (or $\lambda$-matrix). The
              bounds are based on norms of the coefficient matrices and involve
              the inverses of the leading and trailing coefficient matrices.
              They generalize various existing bounds for scalar polynomials
              and single matrices. A variety of tools are used in the
              derivations, including block companion matrices, Gershgorin's
              theorem, the numerical radius, and associated scalar polynomials.
              Numerical experiments show that the bounds can be surprisingly
              sharp on practical problems.",
  keywords = "polynomial eigenvalue problem, $\lambda$-matrix, matrix
              polynomial, block companion matrix, Gershgorin's theorem,
              numerical radius"
        }

@techreport{Tisseur:2000:SQEP,
  author = "Fran\c{c}oise Tisseur and Karl Meerbergen",
  title = "The Quadratic Eigenvalue Problem",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 370",
  pages = "49",
  month = nov,
  year = "2000",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep370.ps.gz",
  abstract = "We survey the quadratic eigenvalue problem,
              treating its many applications, its
              mathematical properties, and a variety of numerical
              solution techniques.  Emphasis is given to exploiting
              both the structure of the matrices in the problem
              (dense, sparse, real, complex, Hermitian,
              skew-Hermitian) and the spectral properties of the
              problem. We classify the available choices of methods
              and catalogue available software.",
  keywords = "quadratic eigenvalue problem, eigenvalue, eigenvector,
             $\lambda$-matrix, matrix polynomial, second-order
             differential equation, overdamped system, gyroscopic
             system, linearization, backward error, pseudospectrum,
             condition number, Krylov methods, Arnoldi method, Lanczos
             method, Jacobi-Davidson method"
}

@techreport{Higham:2000:CNC,
  author = "Nicholas J. Higham",
  title = "Computing the Nearest Correlation Matrix---{A} Problem from Finance",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 369",
  month = oct,
  pages = 14,
  year = 2000,
  note = "Revised July 2001",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep369.pdf",
  abstract = "Given a symmetric matrix what is the nearest correlation matrix,
              that is, the nearest symmetric positive semidefinite matrix with
              unit diagonal? This problem arises in the finance industry, where
              the correlations are between stocks. For distance measured in the
              Frobenius norm, we characterize the solution using convex
              analysis. We show how the modified alternating projections method
              can be used to compute the solution. In the finance application
              the original matrix has many zero or negative eigenvalues; we
              show that the nearest correlation matrix has correspondingly many
              zero eigenvalues and that this fact can be exploited in the
              computation.",
  keywords = "correlation matrix, positive semidefinite matrix, nearness
              problem, convex analysis, alternating projections method"
}

@techreport{Paul:2000:DES,
  author = "Christopher Paul",
  title = "Designing Efficient Software for Solving Delay Differential
           Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 368",
  month = oct,
  pages = 13,
  year = 2000,
  note = "To appear in Journal of Computational Applied Mathematics",
  URL = "http://www.maths.man.ac.uk/~nareports/narep368.ps.gz",
  abstract = "In this paper the efficient implementation of numerical software
              for solving delay differential equations is addressed. Several
              strategies that have been developed over the past 25 years for
              improving the efficiency of delay differential equation solvers
              are described. Of particular interest is a new method of
              automatically constructing the network dependency graph used in
              tracking derivative discontinuities."
}

@techreport{Baker:2000:DAE,
  author = "Christopher Baker and Christopher Paul and Hongjiong Tian",
  title = "Differential Algebraic Equations with After-Effect",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 367",
  month = oct,
  pages = 20,
  year = 2000,
  note = "Submitted to Journal of Computational Applied Mathematics",
  URL = "http://www.maths.man.ac.uk/~nareports/narep367.ps.gz",
  abstract = "We consider the numerical solution of delay differential
             algebraic equations - they are differential algebraic equations
             with after-effect, or constrained delay differential equations.
             The general semi-explicit form of the problem consists of a set
             of delay differential equations along with a set of constraints
             that may involve retarded arguments. Even simply-stated problems
             of this type can give rise to difficult analytical and numerical
             problems. The more tractable examples can be shown to be
             equivalent to systems of delay or neutral delay differential
             equations. We shall explore, in particular through examples, some
             of the complexities and obstacles that can arise when solving
             these problems."
}

@techreport{MCCM:2001:AR,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 1999",
  type =         "Numerical Analysis Report",
  number =       "No. 364",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "24",
  month =        mar,
  year =         "2001",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep364.ps.gz",
}

@techreport{Dedieu:2001:PTH,
  author = "Jean-Pierre Dedieu and Fran\c{c}oise Tisseur",
  title = "Perturbation Theory for Homogeneous Polynomial Eigenvalue Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 363",
  pages = "21",
  month = mar,
  year = "2001",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep363.ps.gz",
  abstract = "We consider polynomial eigenvalue problems
              $P(A,\alpha,\beta)x=0$ in which the matrix polynomial is
              homogeneous in the eigenvalue $(\alpha,\beta)\in\mathbb{C}^2$.
              In this framework infinite eigenvalues are on the same
              footing as finite eigenvalues.  We view the problem in
              projective spaces to avoid normalization of the
              eigenpairs. We show that a polynomial eigenvalue problem
              is well-posed when its eigenvalues are simple. We define
              the condition numbers of a simple eigenvalue
              $(\a,\beta)$ and a corresponding eigenvector $x$ and
              show that the distance to the nearest ill-posed problem
              is equal to the reciprocal of the condition number of
              the eigenvector $x$.  We describe a bihomogeneous Newton
              method for the solution of the homogeneous polynomial
              eigenvalue problem.",
keywords = "polynomial eigenvalue problem, matrix polynomial,
            quadratic eigenvalue problem,condition number"
}

@techreport{Silvester:00:EPB,
  author = "David Silvester and Milan Mihajlovic",
  title = "A black-box multigrid preconditioner for the biharmonic equation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 362",
  month = aug,
  year = 2002,
  note = "Submitted to BIT",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep362.ps.gz",
  abstract = "We examine the convergence characteristics of a preconditioned
   Krylov subspace solver applied to the linear systems arising from low-order
   mixed finite element approximation of the biharmonic problem. The key
   feature of our approach is that the preconditioning can be realized using
   any ``black-box'' multigrid solver designed for the discrete Dirichlet
   Laplacian operator. This leads to preconditioned systems having an
   eigenvalue  distribution consisting of a tightly clustered set together
   with a small number of outliers. Numerical results show that the
   performance of the methodology is competitive with that of specialized
   fast iteration  methods that have been developed in the context of
   biharmonic problems."
}

@techreport{Baker:2000:CTM,
  author = "Christopher T. H. Baker and Evelyn Buckwar",
  title = "Continuous $\Theta$-methods for the Stochastic Pantograph
           Equation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 361",
  month = may,
  year = 2000,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep361.ps.gz",
  abstract = "We consider a stochastic version of the pantograph equation:
              \begin{eqnarray}\nonumber dX(t)  &= & \Big\{ a X(t)+ b X (qt)
              \Big\} dt + \Big\{\sigma_1 + \sigma_2\ X(t) + \sigma_3 X
              (qt)\Big\}dW(t), \\[-1mm]\nonumber && \label{spanto}
              \\[-1mm]\nonumber X(0)  &= & X_0, \nonumber \end{eqnarray} for $
              t \in [0,T]$, a given Wiener noise $W$ and $0 < q < 1$. This is
              an example of an It\^{o} stochastic delay differential equation
              with unbounded memory.
              We give the necessary analytical theory for existence and
              uniqueness of a strong solution of ($\ddagger$), and of strong
              approximations to the solution obtained by a continuous extension
              of the $\Theta$-Euler scheme ($\Theta \in [0,1]$). We establish
              ${\mathcal O}(h^\half)$ mean-square convergence of approximations
              obtained using a bounded mesh of uniform step $h$, rising in the
              case of additive noise only to ${\mathcal O}(h)$. Illustrative
              numerical examples are provided.",
  keywords = "Stochastic delay differential equation, continuous
              $\Theta$-method, mean-square convergence."
}

@techreport{Davies:2000:ACM,
  author = "Philip I. Davies and Nicholas J. Higham and Fran\c{c}oise Tisseur",
  title = "Analysis of the {Cholesky} Method with Iterative Refinement for
           Solving the Symmetric Definite Generalized Eigenproblem",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 360",
  month = jun,
  pages = 21,
  year = "2000",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep360.pdf",
  abstract = "A standard method for solving the symmetric definite generalized
              eigenvalue problem $Ax = \lambda Bx$, where $A$ is symmetric and
              $B$ is symmetric positive definite, is to compute a Cholesky
              factorization $B = LL^T$ (optionally with complete pivoting) and
              solve the equivalent standard symmetric eigenvalue problem $C y =
              \l y$ where $C = L^{-1} A L^{-T}$. Provided that a stable
              eigensolver is used, standard error analysis says that the
              computed eigenvalues are exact for $A+\dA$ and $B+\dB$ with
              $\max( \normt{\dA}/\normt{A}, \normt{\dB}/\normt{B} )$ bounded by
              a multiple of $\kappa_2(B)u$, where $u$ is the unit roundoff. We
              take the Jacobi method as the eigensolver and give a detailed
              error analysis that yields backward error bounds potentially much
              smaller than $\kappa_2(B)u$. To show the practical utility of our
              bounds we describe a vibration problem from structural
              engineering in which $B$ is ill conditioned yet the error bounds
              are small. We show how, in cases of instability, iterative
              refinement based on Newton's method can be used to produce
              eigenpairs with small backward errors. Our analysis and
              experiments also give insight into the popular Cholesky--QR
              method, in which the QR method is used as the eigensolver. We
              argue that it is desirable to augment current implementations of
              this method with pivoting in the Cholesky factorization.",
  keywords = "symmetric definite generalized eigenvalue problem, Cholesky
              method, Cholesky factorization with complete pivoting, Jacobi
              method, backward error analysis, rounding error analysis,
              iterative refinement, Newton's method, LAPACK, MATLAB."
   }

@techreport{Tisseur:2000:SPP,
  author = "Fran\c{c}oise Tisseur and Nicholas J. Higham",
  title = "Structured Pseudospectra for Polynomial Eigenvalue Problems,
           with Applications",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 359",
  month = apr,
  pages = 22,
  year = "2000",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep359.pdf",
  abstract = "Pseudospectra associated with the standard and generalized
              eigenvalue problems have been widely investigated in recent
              years. We extend the usual definitions in two respects, by
              treating the polynomial eigenvalue problem and by allowing
              structured perturbations of a type arising in control theory. We
              explore connections between structured pseudospectra, structured
              backward errors, and structured stability radii. Two main
              approaches for computing pseudospectra are described. One is
              based on a transfer function and employs a generalized Schur
              decomposition of the companion form pencil. The other, specific
              to quadratic polynomials, finds a solvent of the associated
              quadratic matrix equation and thereby factorizes the quadratic
              $\lambda$-matrix. Possible approaches for large, sparse problems
              are also outlined. A collection of examples from vibrating
              systems, control theory, acoustics and fluid mechanics is given
              to illustrate the techniques.",
  keywords = "polynomial eigenvalue problem, $\lambda$-matrix, matrix
              polynomial, pseudospectrum, stability radius, backward error,
              transfer function, quadratic matrix equation, solvent, structured
              perturbations, Orr--Sommerfeld equation"
   }

@techreport{Higham:2000:EPA,
  author = "Nicholas J. Higham",
  title = "Evaluating {Pad\'e} Approximants of the Matrix Logarithm",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 358",
  month = mar,
  pages = 10,
  year = 2000,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep358.pdf",
  abstract = "The inverse scaling and squaring method for evaluating the
              logarithm of a matrix takes repeated square roots to bring the
              matrix close to the identity, computes a Pad\'e approximant, and
              then scales back. We analyze several methods for evaluating the
              Pad\'e approximant, including Horner's method (used in some
              existing codes), suitably customized versions of the
              Paterson--Stockmeyer method and Van Loan's variant, and methods
              based on continued fraction and partial fraction expansions. The
              computational cost, storage, and numerical accuracy of the
              methods are compared. We find the partial fraction method to be
              the best method overall and illustrate the benefits it brings to
              a transformation-free form of the inverse scaling and squaring
              method recently proposed by Cheng, Higham, Kenney and Laub. We
              comment briefly on how the analysis carries over to the matrix
              exponential.",
  keywords = "matrix logarithm, Pad\'e approximation, inverse scaling and
              squaring method, Horner's method, Paterson--Stockmeyer method,
              continued fraction, partial fraction expansion"
}

@techreport{Tisseur:2000:SSHE,
  author = "Fran\c{c}oise Tisseur",
  title = "Stability of Structured {Hamiltonian} Eigensolvers",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 357",
  pages = "23",
  month = feb,
  year = "2000",
  URL =          "ftp://ftp.ma.man.ac.uk/pub/narep/narep357.ps.gz",
  abstract = "Various applications give rise to eigenvalue problems
              for which the matrices are Hamiltonian or
              skew-Hamiltonian and also symmetric or skew-symmetric.
              We define structured backward errors that are useful for
              testing the stability of numerical methods for the
              solution of these four classes of structured
              eigenproblems. We introduce the symplectic quasi-QR
              factorisation and show that for three of the classes it
              enables the structured backward error to be efficiently
              computed.  We also give a detailed rounding error
              analysis of some recently developed Jacobi-like
              algorithms of Fa{\ss}bender, Mackey \& Mackey for these
              eigenproblems.  Based on the direct solution of $4\times
              4$, and in one case $8\times 8$, structured subproblems
              these algorithms produce a complete basis of symplectic
              orthogonal eigenvectors for the two symmetric cases and
              a symplectic orthogonal basis for all the real invariant
              subspaces for the two skew-symmetric cases. We prove
              that, when the rotations are implemented using suitable
              formulae, the algorithms are strongly backward stable
              and we show that the QR algorithm does not have this
              desirable property",
keywords = "Backward error, Hamiltonian, skew-Hamiltonian,
            symmetric-Hamiltonian, skew-symmetric skew-Hamiltonian,
            structure-preserving, rounding error, symplectic,
            quaternion rotation"
}

@techreport{Sturgeon:2000:SSD,
  author = "J. A. Sturgeon and R. M. Thomas and I. Gladwell",
  title  = "Solving a Singular DAE Model of Unconfined Detonation",
  type   = "Numerical Analysis Report",
  number = "No. 356",
  institution = inst-mccm,
  address = inst-mccm:adr,
  month  = apr,
  year   = "2000",
  pages = "27",
  URL    = "http://www.maths.man.ac.uk/~nareports/narep356.ps.gz",
  abstract = "We consider a simplified model of an unconfined one-dimensional
              detonation problem, giving a brief survey of the history of the
              problem and of its numerical solution.  This problem with its
              mathematical features is typical of those solved commercially
              by ICI plc, and the specific values used for the chemical
              constants in the example are typical of those of interest.
              Unfortunately, not all obvious methods work well, because of the
              singular nature of the problem at the Chapman-Jouguet shock
              front.  We concentrate on shooting methods for the detonation
              problem based on backward differentiation formula integrators,
              and present a new analysis which explains how these methods
              work.  Finally, we outline some possibilities for further work,
              including discussing a more general detonation problem,
              previous solutions and potential future solution methods",
  keywords = "detonation problem, differential-algebraic equations, backward
              differentiation formulae, shooting methods, numerical analysis,
              Chapman-Jouguet shock front"
    }

@techreport{Langlois:1999:ALC,
  author = "Philippe Langlois",
  title = "Automatic Linear Correction of Rounding Errors",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 355",
  month = nov,
  pages = "24",
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep355.ps.gz",
  abstract = "A new automatic method to correct the first-order effect of
              floating point rounding errors on the result of a numerical
              algorithm is presented. A correcting term and a confidence
              threshold are computed using automatic differentiation,
              computation of elementary rounding error and running error
              analysis. Algorithms for which the accuracy of the result is not
              affected by higher order terms are identified. The correction is
              applied to the final result or to sensitive intermediate results.
              The properties and the efficiency of the method are illustrated
              with a sample numerical example.",
  keywords = "automatic error analysis, rounding error, floating point
              arithmetic"
        }

@techreport{Davies:1999:NSG,
  author = "Philip I. Davies and Nicholas J. Higham",
  title = "Numerically Stable Generation of
           Correlation Matrices and their Factors",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 354",
  month = nov,
  pages = "15",
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep354.pdf",
  abstract = "Correlation matrices---symmetric positive semidefinite matrices
              with unit diagonal---are important in statistics and in numerical
              linear algebra. For simulation and testing it is desirable to be
              able to generate random correlation matrices with specified
              eigenvalues (which must be nonnegative and sum to the dimension
              of the matrix). A popular algorithm of Bendel and Mickey takes a
              matrix having the specified eigenvalues and uses a finite
              sequence of Given rotations to introduce $1$s on the diagonal. We
              give improved formulae for computing the rotations and prove that
              the resulting algorithm is numerically stable. We show by example
              that the formulae originally proposed, which are used in certain
              existing Fortran implementations, can lead to serious
              instability. We also show how to modify the algorithm to generate
              a rectangular matrix with columns of unit 2-norm. Such a matrix
              represents a correlation matrix in factored form, which can be
              preferable to representing the matrix itself, for example when
              the correlation matrix is nearly singular to working precision.",
  keywords = "random correlation matrix, Bendel--Mickey algorithm, eigenvalues,
              singular value decomposition, test matrices, forward error
              bounds, relative error bounds, IMSL, NAG Library, Jacobi method"
        }

@techreport{Cheng:1999:ALM,
  author = "Sheung Hun Cheng and Nicholas J. Higham and Charles S. Kenney
            and Alan J. Laub",
  title = "Approximating the Logarithm of a Matrix to Specified Accuracy",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 353",
  pages = "14",
  month = nov,
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep353.ps.gz",
  abstract = "The standard inverse scaling and squaring algorithm for computing
              the matrix logarithm begins by transforming the matrix to Schur
              triangular form in order to facilitate subsequent matrix square
              root and Pad\'e approximation computations. A transformation-free
              form of this method is presented that exploits incomplete
              Denman--Beavers square root iterations and aims for a specified
              accuracy. The error introduced by using approximate square roots
              is accounted for by a novel splitting lemma for logarithms of
              matrix products. The number of square root stages and the degree
              of the final Pad\'e approximation are chosen to minimize the
              computational work. This new method is attractive for
              high-performance computation since it uses only the basic
              building blocks of matrix multiplication, LU factorizati on and
              matrix inversion.",
  keywords = "matrix logarithm, Pad\'e approximation, inverse scaling and
              squaring method, matrix square root, Denman--Beavers iteration"
        }

@techreport{Silvester:99:EPL,
  author = "David Silvester and Howard Elman and David Kay and Andrew Wathen",
  title = "Efficient Preconditioning of the Linearized {Navier-Stokes}
           Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 352",
  month = oct,
  year = "1999",
  note = "Solicited by the Journal of Computational and Applied Mathematics",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep352.ps.gz",
  abstract = "We outline a new class of robust and efficient methods for
              solving  subproblems that arise in the linearization and
              operator splitting of Navier-Stokes equations. We describe
              a very general strategy for preconditioning that has two
              basic building blocks; a multigrid V-cycle for the
              scalar convection-diffusion operator, and a multigrid V-cycle
              for a pressure Poisson operator. We present numerical
              experiments illustrating that  a simple implementation
              of our approach leads to an effective and robust solver
              strategy in that the convergence rate is independent of the grid,
              robust with respect to the time-step, and only deteriorates
              very slowly as the Reynolds number is increased."
}

@techreport{Ford:1999:HDN,
  author = "Neville J. Ford and Volker Wulf",
  title = "How Do Numerical Methods Perform for Delay Differential Equations
           Undergoing a {Hopf} Bifurcation?",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 351",
  pages = "10",
  month = sep,
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep351.ps.gz",
  abstract = "In this report we consider the numerical solution of delay
              differential equations that undergo a Hopf bifurcation. We show
              that, under suitable conditions on the choice of method, the Hopf
              bifurcation is maintained in the approximation, and that the
              class of bifurcation is also preserved.",
  keywords = "Delay differential equations, numerical methods, Hopf
              bifurcation"
  }

@techreport{Ford:1999:Charroots,
  author = "Neville J. Ford",
  title = "Numerical Approximation of the Characteristic Values for a Delay
           Differential Equation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 350",
  pages = "11",
  month = sep,
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep350.ps.gz",
  abstract = "In this report we consider the characteristic roots of numerical
              approximations of delay differential equations and their
              relationship to the characteristic values of the original
              problem. We are able to show that, as h -> 0, the characteristic
              roots are approximated to the order of the method.",
  keywords = "Delay differential equations, characteristic roots, numerical
              methods, stability coefficient"
        }

@techreport{Rihan:99:SAP,
  author = "Christopher T. H. Baker and Fathalla A. Rihan",
  title = "Sensitivity Analysis of Parameters in Modelling With
           Delay-Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 349",
  month = aug,
  pages = 22,
  year = 1999,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep349.ps.gz",
  abstract = "Many problems in bioscience for which
	observations are reported in the literature can be modelled by
	suitable functional differential equations incorporating a delay,
	parameterized by parameters $ p_1$,$p_2$,$\dots$,$p_L$.  Given such
	observations (which usually contain error or 'noise'), we may
	determine the parameters $\{p_i\}$ by fitting the model equations to
	the data and optimizing a measure of best fit.
	Our aim in this paper is to produce a new method to estimate (i) the
	sensitivity of the state variables to the parameter estimates $\{
	p_i\}$, (ii) the sensitivity of the parameter estimates to the
	observations and (iii) the nonlinearity effects for delay differential
	models.  The sensitivity of the parameter estimate to the observation
	is {\it low} if the sensitivity of the state variable to the parameter
	estimate is {\it high}.  Sensitivity coefficients are used to
	determine the covariance matrix of parameter estimates and hence to
	determine the standard deviations.  Numerical results are used to
	illustrate the results.",
  keywords = "Sensitivity analysis, parameter estimates, neutral
	      delay differential equation, time-lag, nonlinearity effect."
       }

@techreport{Norburn:99:FAS,
  author = "Sean Norburn and David Silvester ",
  title = "Fourier Analysis of Stabilised {Q1-Q1} Mixed
	   Finite Element Approximation",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 348",
  month = sep,
  year = 1999,
  note = "Submitted to the SIAM Journal on Numerical Analysis",
  URL = "http://www.maths.man.ac.uk/~nareports/narep348.ps.gz",
  abstract = "We use Fourier analysis  to investigate the instability
	      of an equal-order mixed finite element approximation
              method for elliptic incompressible flow equations.
              The lack of stability can be attributed to the fact that
	      the associated discrete LBB (Ladyzhenskaya-Babuska-Brezzi)
              constant tends to zero as the mesh size is reduced.
              We develop a stabilisation approach that is appropriate
	      to the periodic setting, and deduce optimal choices of
	      the associated stabilisation parameter."
}

@techreport{Higham:99:NAQ,
  author = "Nicholas J. Higham and Hyun-Min Kim",
  title = "Numerical Analysis of a Quadratic Matrix Equation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 347",
  month = aug,
  pages = 26,
  year = 1999,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep347.pdf",
  abstract = "The quadratic matrix equation $AX^2 + BX + C = 0$ in $n\times n$
              matrices arises in applications and is of intrinsic interest as
              one of the simplest nonlinear matrix equations. We give a
              complete characterization of solutions in terms of the
              generalized Schur decomposition and describe and compare various
              numerical solution techniques. In particular, we give a thorough
              treatment of functional iteration methods based on Bernoulli's
              method. Other methods considered include Newton's method with
              exact line searches, symbolic solution and continued fractions.
              We show that functional iteration applied to the quadratic matrix
              equation can provide an efficient way to solve the associated
              quadratic eigenvalue problem $(\lambda^2 A + \lambda B +
              C)x=0$.",
  keywords = "quadratic matrix equation, solvent, generalized Schur
              decomposition, scaling, functional iteration, Bernoulli's method,
              Newton's method, exact line searches, continued fractions,
              quadratic eigenvalue problem",
       }

@techreport{Tisseur:99:NMFPA,
  author = "Fran\c{c}oise Tisseur",
  title = "Newton's Method in Floating Point Arithmetic
           and Iterative Refinement of Generalized Eigenvalue Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 346",
  pages = "26",
  month = aug,
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep346.ps.gz",
  abstract = "We examine the behaviour of Newton's method in floating
             point arithmetic, allowing for extended precision in
             computation of the residual, inaccurate evaluation of the
             Jacobian and unstable solution of the linear systems. We
             bound the limiting accuracy and the smallest norm of the
             residual.  The application that motivates this work is
             iterative refinement for the generalized eigenvalue
             problem.  We show that iterative refinement by Newton's
             method can be used to improve the forward and backward
             errors of computed eigenpairs.",
keywords = " Newton's method, generalized eigenvalue problem, iterative
             refinement, backward error, forward error, rounding error
             analysis, limiting accuracy, limiting residual"
     }

@techreport{Baker:99:INA,
  author = "Christopher T. H. Baker and Evelyn Buckwar",
  title = "Introduction to the Numerical Analysis of Stochastic Delay
           Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 345",
  month = aug,
  year = 1999,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep345.ps.gz",
  abstract = "We consider the problem of strong approximations of
              the solution of stochastic delay differential equations
              of It\^{o} form
              \[
                \begin{array}{rcl}
                 dX(t)  &= &  f(t, X(t),X (t- \tau) ) \; dt
                        +  g(t, X(t),X (t- \tau) )\; dW(t),
                            \ t \in [0,T] \\
                  X(t)  &= & \Psi(t),  \qquad  t \in [-\tau,0],
                \end{array}
              \]
              for given $f, g$, Wiener noise $W$ and given $\tau > 0$,
              with a prescribed initial function $\Psi$. We indicate the
              nature of the equations of interest and the problems
              that arise in a discussion of their numerical treatment,
              provide necessary background material, and give a
              convergence proof for an Euler-Maruyama scheme. Some
              illustrative numerical examples are provided."
         }

@techreport{MCCM:2000:AR,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 1998",
  type =         "Numerical Analysis Report",
  number =       "No. 344",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "24",
  month =        mar,
  year =         "2000",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep344.ps.gz",
}

@techreport{Rihan:99:RUD,
  author = "Christopher T. H. Baker and Gennadii A. Bocharov and
            Fathalla A. Rihan",
  title = "A Report on the Use of Delay Differential Equations
           in Numerical Modelling in the Biosciences",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 343",
  month = aug,
  pages = 45,
  year = 1999,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep343.ps.gz",
  abstract = "We review the application of numerical techniques to investigate
	mathematical models of phenomena in the biosciences using delay
	differential equations. We show that there are {\em prima facie}
	reasons for using such models: $(i)$ they have a richer mathematical
	framework (compared with ordinary differential equations) for the
	analysis of biosystems dynamics, $(ii)$ they display better
	consistency with the nature of the underlying processes and predictive
	results. We now have suitable computational techniques to treat
	numerically the emerging models for the biosciences.",
  keywords = "Delay differential equations, dynamical systems,
	biological systems, deterministic and stochastic models,
	numerical modelling, parameter estimation, optimization."
       }

@techreport{Edwards:1999:EBS,
  author = "John T. Edwards and Jason A. Roberts",
  title = "On the Existence of Bounded Solutions of Discrete Analogues for a
           Nonlinear Integrodifferential Equation",
  institution = inst-MCCM,
  address = inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 342",
  pages = "14",
  month = sep,
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep342.ps.gz",
  abstract = "This work is concerned with the boundedness and  ultimately the
              asymptotic stability of solutions of discrete analogues of a
              nonlinear integrodifferential equation with convolution kernel.
              In previous work, the stability analysis of the nonlinear
              equation required the assumption that bounded solutions exist. We
              prove in this work that, for a certain choice of discretization
              of the equation, all solutions are                   bounded (and
              hence asymptotically stable) provided the equation's parameters
              fall within certain ranges.",
  keywords = "Volterra integrodifferential equations, numerical methods,
              difference equations, stability analysis"
       }

@techreport{Higham:99:BAM,
  author = "Nicholas J. Higham and Fran\c{c}oise Tisseur",
  title = "A Block Algorithm for Matrix 1-Norm Estimation,
           with an Application to 1-Norm Pseudospectra",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 341",
  month = apr,
  pages = 32,
  year = 1999,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep341.ps.gz",
  abstract = "The matrix 1-norm estimation algorithm used in LAPACK and various
              other software libraries and packages has proved to be a valuable
              tool. However, it has the limitations that it offers the user no
              control over the accuracy and reliability of the estimate and
              that it is based on level~2 BLAS operations. A block
              generalization of the 1-norm power method underlying the
              estimator is derived here and developed into a practical
              algorithm applicable to both real and complex matrices. The
              algorithm works with $n\times t$ matrices, where $t$ is a
              parameter. For $t=1$ the original algorithm is recovered, but
              with two improvements (one for real matrices and one for complex
              matrices). The accuracy and reliability of the estimates
              generally increases with $t$ and the computational kernels are
              level~3 BLAS operations for $t>1$. The last $t-1$ columns of the
              starting matrix are randomly chosen, giving the algorithm a
              statistical flavour. As a by-product of our investigations we
              identify a matrix for which the 1-norm power method takes the
              maximum number of iterations. As an application of the new
              estimator we show how it can be used to efficiently approximate
              1-norm pseudospectra.",
  keywords = "matrix 1-norm, matrix norm estimation, matrix condition number,
              condition number estimation, $p$-norm power method, 1-norm
              pseudospectrum, LAPACK, level 3 BLAS"
       }

@techreport{Ford:1999:NHBClass,
  author = "Neville J. Ford and Volker Wulf",
  title = "Numerical {Hopf} Bifurcation for a Class of
           Delay Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 340",
  pages = "17",
  month = feb,
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep340.ps.gz",
  abstract = "In this report we consider discretisation of parameter-dependent
              delay differential equations. We show that, if the delay
              differential equation undergoes a Hopf bifurcation, then the
              discrete scheme undergoes a Hopf bifurcation of the same type.
              The results of this report extend our previous results for the
              delay logistic equation to a wider class of problems.",
  keywords = "Delay differential equations, Hopf bifurcation, numerical methods,
              stability coefficient"
        }

@techreport{Higham:99:SQM,
  author = "Nicholas J. Higham and Hyun-Min Kim",
  title = "Solving a Quadratic Matrix Equation by
           {Newton's} Method with Exact Line Searches",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 339",
  month = jan,
  pages = 18,
  year = 1999,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep339.pdf",
  abstract = "We show how to incorporate exact line searches into Newton's
              method for solving the quadratic matrix equation $AX^2 + BX + C =
              0$, where $A$, $B$ and $C$ are square matrices. The line searches
              are relatively inexpensive and improve the global convergence
              properties of Newton's method in theory and in practice. We also
              derive a condition number for the problem and show how to compute
              the backward error of an approximate solution.",
  keywords = "quadratic matrix equation, solvent, Newton's method, generalized
              Sylvester equation, exact line searches, quadratic eigenvalue
              problem, condition number, backward error"
       }

@techreport{Davies:99:GTM,
  author = "Philip I. Davies and Nicholas J. Higham",
  title = "Generating Test Matrices for the One- and Two-Sided {Jacobi} Methods",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 338",
  month = jan,
  pages = "13",
  year = 1999,
  URL =          "http://www.maths.man.ac.uk/~nareports/narep338.ps.gz",
  abstract = "For a symmetric positive definite matrix the relative error in
              the eigenvalues computed by the two-sided Jacobi method is
              bounded in terms of the condition number of the matrix scaled to
              have unit diagonal. Similarly, for a general matrix the relative
              error in the singular values computed by the one-sided Jacobi
              method is bounded in terms of the condition number of the matrix
              scaled to have rows or columns of unit-2 norm. We show how to
              generate random matrices having these scalings and given
              eigenvalues and singular values, respectively. For the two-sided
              Jacobi method we apply an algorithm of Bendel and Mickey for
              generating random correlation matrices, with an improved formula
              for the rotations. We show how to modify the algorithm to
              generate matrices for the one-sided Jacobi method. Using these
              test matrices we show empirically that the forward error bounds
              for the one- and two-sided Jacobi methods are sharp.",
  keywords = "Jacobi method, correlation matrix, eigenvalues, singular value
              decomposition, test matrices, forward error bounds, relative
              error bounds"
       }

@techreport{Paul:99:TDDDE,
  author = "Christopher A.H. Paul",
  title = "The Treatment of Derivative Discontinuities in Differential
           Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 337",
  month = jan,
  pages = 11,
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep337.ps.gz",
  abstract = "The assumption of sufficiently smooth derivatives underlies much
              of the analysis of numerical methods for both ordinary and delay
              differential equations. However, derivative discontinuities can
              arise in ordinary differential equations and usually do arise in
              delay differential equations. In this paper, I review two of the
              approaches for treating derivative discontinuities, ``tracking''
              using the discontinuity tracking equations and ``detection and
              location'' using defect error control. I conclude that neither
              approach is ideal when it comes to the treatment of derivative
              discontinuities in delay differential algebraic equations.",
  keywords = "differential equations, discontinuities, numerical solution"
        }

@techreport{Higham:99:NSM,
  author = "Nicholas J. Higham",
  title = "A New \texttt{sqrtm} for \textsc{Matlab}",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 336",
  month = jan,
  pages = 11,
  year = "1999",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep336.ps.gz",
  abstract = "\Matlab's function \t{sqrtm} computes a square root of a matrix.
              We propose a replacement for the \t{sqrtm} in \Matlab~5.2 that is
              more accurate and returns useful information about the stability
              and conditioning of the problem.",
  keywords = "matrix square root, \Matlab, Schur decomposition, condition
              number, stability"
        }

@techreport{Tshelametse:98:ANT,
  author = "Ronald Tshelametse and Jack Williams",
  title = "A New Termination Criterion for Nonlinear Iterations in {ODE} Codes",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 334",
  month = oct,
  pages = 23,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep334.ps.gz",
  abstract = "Many codes that solve systems of stiff initial value problems for
              ordinary differential equations implement implicit linear
              multistep methods. This yields nonlinear equations which are
              commonly solved using some form of Newton iteration. For linear
              multistep methods we apply the theory of Dorsselaer and Spijker
              \cite{dorss94} to justify a new stopping criterion for the Newton
              iteration. The ideas are tested on the \textsc{Matlab} ODE suite
              of Shampine and Reichelt \cite{shamp97} and are shown to lead to
              a significant improvement in overall efficiency for some
              problems.",
  keywords = "stiff problems, IVPs, ODEs, BDFs, Newton method"
        }

@techreport{Higham:98:NAS,
  author = "Nicholas J. Higham",
  title = "Notes on Accuracy and Stability of Algorithms in
           Numerical Linear Algebra",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 333",
  month = aug,
  pages = 42,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep333.ps.gz",
  note = "To appear in proceedings of EPSRC Numerical Analysis Summer School,
          Leicester University, July 1998."
        }

@techreport{Tisseur:1998:BEC,
  author = "Fran\c{c}oise Tisseur",
  title = "Backward Error and Condition of Polynomial Eigenvalue Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 332",
  pages = "20",
  month = aug,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep332.ps.gz",
  abstract = "We develop normwise backward errors and condition
              numbers for the polynomial eigenvalue problem (PEP).
              The standard way of dealing
              with the PEP is to reformulate it as a generalized
	      eigenvalue problem (GEP). For
              the special case of the quadratic eigenvalue problem (QEP), we
              show that solving the QEP by applying the QZ algorithm to a
	      corresponding GEP can be backward unstable.  The QEP can
	      be reformulated as a GEP in many ways. We investigate
	      the sensitivity of a given eigenvalue to perturbations
	      in each ofthe GEP formulations and identify which
	      formulations are to be preferred for large and small
              eigenvalues respectively.",
  keywords = "Polynomial eigenvalue problem,
              quadratic eigenvalue problem, generalized eigenvalue
	      problem, backward error, condition number"
 }

@techreport{Barlow:1998:MAB,
  author = "Jesse L. Barlow",
  title = "More Accurate Bidiagonal Reduction for Computing the Singular Value
           Decomposition",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 331",
  month = aug,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep331.ps.gz",
  abstract = "Bidiagonal reduction is the preliminary stage for the fastest
              stable algorithms for computing the singular value decomposition.
              However, the best error bounds on bidiagonal reduction methods
              are of the form \[ A + \delta A = UBV^T , \] \[ \normtwo{\delta A
              } \leq \macheps f(n) \normtwo{A} \] where $B$ is bidiagonal, $U$
              and $V$ are orthogonal, $\macheps$ is machine precision, and
              $f(n)$ is a modestly growing function of the dimensions of $A$.
              A Givens-based bidiagonal reduction procedure is proposed that satisfies
              \[
                  A + \delta A = U (B + \delta B ) V^T,
              \]
              where $\delta B$ is bounded {\em componentwise} and $\delta A$
              satisfies a tighter columnwise bound. Thus the routine obtains
              more accurate singular values for matrices that have poor column
              scaling or those arising from rank revealing decompositions.",
  keywords = ""
        }

@techreport{Cheng:1998:NDP,
  author = "Sheung Hun Cheng and Nicholas J. Higham",
  title = "The Nearest Definite Pair for the
           {Hermitian} Generalized Eigenvalue Problem",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 325",
  pages = "15",
  month = may,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep325.ps.gz",
  abstract = "The generalized eigenvalue problem $Ax = \lambda Bx$ has special
              properties when $(A,B)$ is a Hermitian and definite pair. Given a
              general Hermitian pair $(A,B)$ it is of interest to find the
              nearest definite pair having a specified Crawford number $\delta
              > 0$. We solve the problem in terms of the inner numerical radius
              associated with the field of values of $A+iB$. We show that once
              the problem has been solved it is trivial to rotate the perturbed
              pair $(A+\dA,B+\dB)$ to a pair $(\widetilde{A},\widetilde{B})$
              for which $\lambda_{\min}(\widetilde{B})$ achieves its maximum
              value $\delta$, which is a numerically desirable property when
              solving the eigenvalue problem by methods that convert to a
              standard eigenvalue problem by ``inverting $B$''. Numerical
              examples are given to illustrate the analysis.",
  keywords = "Nearest definite pair, Crawford number, Hermitian pair,
              generalized eigenvalue problem, field of values, inner numerical
              radius, numerical radius"
   }

@techreport{Higham:1998:QRF,
  author = "Nicholas J. Higham",
  title = "{QR} Factorization with Complete Pivoting
           and Accurate Computation of the {SVD}",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 324",
  month = sep,
  pages = 26,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep324.ps.gz",
  abstract = "A new algorithm of Demmel et al.\ for computing the singular
              value decomposition (SVD) to high relative accuracy begins by
              computing a rank-revealing decomposition. Demmel et al.\ analyse
              the use of Gaussian elimination with complete pivoting (GECP) for
              computing the rank-revealing decomposition. We investigate the
              use of QR factorization with complete pivoting (that is, column
              pivoting together with row sorting or row pivoting) as an
              alternative to GECP, since this leads to a faster SVD algorithm.
              We derive a new componentwise backward error result for
              Householder QR factorization and combine it with the theory of
              Demmel et al.\ to show when high relative accuracy in the
              computed SVD can be expected. An a posteriori error bound is
              derived that gives useful estimates of the relative accuracy of
              the computed singular values. Numerical experiments confirm the
              theoretical predictions.",
  keywords = "QR factorization, Householder matrix, row pivoting, row sorting,
              column pivoting, complete pivoting, backward error analysis,
              singular value decomposition, relative accuracy"
     }

@techreport{Ford:1998:NHB,
  author = "Neville J. Ford and Volker Wulf",
  title = "Numerical {Hopf} Bifurcation for the Delay Logistic Equation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 323",
  pages = "20",
  month = apr,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep323.ps.gz",
  abstract = "This report is concerned with the approximate solution of the
              delay logistic equation close to the Hopf bifurcation point. We
              show that, for sufficiently small step sizes the bifurcation
              point of the original problem is approximated by a bifurcation
	      point in the discrete problem and that furthermore, the nature of
              the bifurcation is preserved.",
  keywords = "Delay differential equations, Hopf bifurcation, numerical
              methods, stability coefficient"
        }

@techreport{Ford:1998:BLP,
  author = "Neville J. Ford and Volker Wulf",
  title = "The Use of Boundary Locus Plots in the Identification of Bifurcation
           Points in Numerical Approximation of Delay Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 322",
  pages = "15",
  month = mar,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep322.ps.gz",
  abstract = "We use the boundary locus method as a tool to determine precise
             parameter values at which discrete analogues of delay
	     differential equations may have a Hopf bifurcation. We prove that,
	     for strictly stable multistep methods, the boundary of the
	     stability region for the discrete equation approximates the true
	     stability region to the order of the method and that Hopf
	     bifurcation points in the discrete scheme approximate the
	     bifurcation points for the original equation to the same order.",
  keywords = "Delay differential equations, Hopf bifurcations,
              Boundary locus method"
          }

@techreport{Cox:1998:BEB,
  author = "Anthony J. Cox and Nicholas J. Higham",
  title = "Backward Error Bounds for Constrained Least Squares Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 321",
  pages = "20",
  month = apr,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep321.ps.gz",
  abstract = "We derive an upper bound on the normwise backward error of an
              approximate solution to the equality constrained least squares
              problem $\min_{Bx=d}\|b-Ax\|_2$. Instead of minimizing over the
              four perturbations to $A$, $b$, $B$ and $d$, we fix those to $B$
              and $d$ and minimize over the remaining two; we obtain an
              explicit solution of this simplified minimization problem. Our
              experiments show that backward error bounds of practical use are
              obtained when $B$ and $d$ are chosen as the optimal normwise
              relative backward perturbations to the constraint system, and we
              find that when the bounds are weak they can be improved by direct
              search optimization. We also derive upper and lower backward
              error bounds for the problem of least squares minimization over a
              sphere: $\min_{\|x\|_2\le\alpha}\|b-Ax\|_2$.",
  keywords = "Equality constrained least squares problem, least squares
              minimization over a sphere, null space method, elimination
              method, method of weighting, backward error, backward stability"
    }

@techreport{MCCM:1998:AR,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 1997",
  type =         "Numerical Analysis Report",
  number =       "No. 320",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages =        "20",
  month =        mar,
  year =         "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep320.ps.gz",
}

@techreport{Cox:1998:RWB,
  author = "Anthony J. Cox and Nicholas J. Higham",
  title = "Row-Wise Backward Stable Elimination Methods
           for the Equality Constrained Least Squares Problem",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 319",
  pages = "18",
  month = mar,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep319.ps.gz",
  abstract = " It is well known that the solution of the equality constrained
               least squares (LSE) problem $\min_{Bx=d}\|b-Ax\|_2$ is the limit
               of the solution of the unconstrained weighted least squares
               problem $$ \min_x\left\| \bmatrix{ \mu d \cr b } - \bmatrix{\mu
               B \cr A } x \right\|_2 $$ as the weight $\mu$ tends to infinity,
               assuming that $\bmatrix{B^T & A^T \cr}^T$ has full rank. We
               derive a new method for the LSE problem by applying Householder
               QR factorization with column pivoting to this weighted problem
               and taking the limit analytically, with an appropriate rescaling
               of rows. The method obtained is a type of direct elimination
               method. We adapt existing error analysis for the unconstrained
               problem to obtain a row-wise backward error bound for the
               method. The bound shows that, provided row pivoting or row
               sorting is used, the method is well-suited to problems in which
               the rows of $A$ and $B$ vary widely in norm. As a by-product of
               our analysis, we show that the standard elimination method for
               solving the LSE problem satisfies a row-wise backward error
               bound of precisely the same form as the new method. We
               illustrate our results with numerical tests.",
  keywords = "Constrained least squares problem, weighted least squares
              problem, Householder QR factorization, Gaussian elimination,
              elimination method, rounding error analysis, backward stability,
              row pivoting, row sorting, column pivoting"
 }

@techreport{Tisseur:1998:PDC,
  author = "Fran\c{c}oise Tisseur and Jack J. Dongarra",
  title = "Parallelizing the Divide and Conquer Algorithm for
           the Symmetric Tridiagonal Eigenvalue Problem on
           Distributed Memory Architectures",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 317",
  pages = "28",
  month = mar,
  year = "1998",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep317.ps.gz",
  abstract = "We present a new parallel implementation of a divide and
              conquer algorithm for computing the spectral
              decomposition of a symmetric tridiagonal matrix on
              distributed memory architectures.  The implementation we
              develop differs from other implementations in that we use a two
              dimensional block cyclic distribution of the data, we
              use the L{\"o}wner theorem approach to compute orthogonal
              eigenvectors, and we introduce permutations before the
              back transformation of each rank-one update in order to
              make good use of deflation.  This algorithm yields the
              first scalable, portable and numerically stable parallel
              divide and conquer eigensolver.  Numerical results
              confirm the effectiveness of our algorithm.  We compare
              performance of the algorithm with that of the QR
              algorithm and of bisection followed by inverse iteration
              on an IBM SP2 and a cluster of Pentium PII's.",
  keywords = "Divide and conquer, symmetric eigenvalue problem,
              tridiagonal matrix, rank-one modification, parallel algorithm,
              ScaLAPACK, LAPACK, distributed memory architecture"
 }

@techreport{kay:97:APE,
  author = "David Kay and David Silvester",
  title = "A-Posteriori Error Estimation for Stabilised Mixed Approximations of
           the {Stokes} Equations.",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 316",
  month = dec,
  year = 1997,
  note = "Submitted to the SIAM Journal on Scientifc Computing",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep316.ps.gz",
  abstract = "The practical implementation  of the lowest-order \plpo\ (linear
              velocity, constant pressure) finite element method for the
              steady-state incompressible (Navier-)Stokes equations is
              addressed in this work. Three different types of a-posteriori
              error indicator are introduced and each is shown to be equivalent
              to the discretisation  error. Our numerical results show that
              these indicators can be used to drive an adaptive refinement
              process which is  specially tailored to create grids which
              conform to the requirements of the local stabilisation. It is
              also shown that the indicators provide an effective method for
              detecting local singularities in the flow solution."
   }

@techreport{Baker:1997:MMI2,
  author = "Christopher T. H. Baker and Gennadii A. Bocharov and
            Christopher A. H. Paul",
  title = "Mathematical Modelling Of The Interleukin-2 T-Cell System: A
           Comparative Study Of Approaches Based On Ordinary And Delay
           Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 314",
  pages = "17",
  month = nov,
  year = "1997",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep314.ps.gz",
  abstract = "Cell proliferation and differentiation phenomena are key issues
              in immunology, tumour growth and cell biology. We study the
              kinetics of cell growth in the immune system using mathematical
              models formulated in terms of ordinary and delay differential
              equations. We study how the suitability of the mathematical
              models depends on the nature of the cell growth data and the
              types of differential equations by minimizing an objective
              function to give a best-fit parameterized solution. We show that
              mathematical models that incorporate a time-lag in the cell
              division phase are more consistent with certain reported data.
              They also allow various cell proliferation characteristics to be
              estimated directly, such as the average cell-doubling time and
              the rate of commitment of cells to cell division. Specifically,
              we study the Interleukin-2-dependent cell division of
              phytohemagglutinin stimulated T-cells -- the model of which can
              be considered to be a general model of cell growth. We also
              review the numerical techniques available for solving delay
              differential equations and calculating the least-squares
              best-fit parameterized solution.",
  keywords = "Cell proliferation, Interleukin-2, mathematical modelling,
              parameter estimation, time-lag"
        }

@techreport{Baker:1997:MATL,
  author = "Christopher T. H. Baker and Gennadii A. Bocharova and
            Christopher A. H. Paul and Fathalla A. Rihan",
  title = "Modelling and Analysis of Time-Lags in Cell Proliferation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 313",
  pages = "16",
  month = nov,
  year = "1997",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep313.ps.gz",
  abstract = "In this paper, we present a systematic approach for obtaining
              qualitatively and quantitatively correct mathematical models of
              some biological phenomena with time-lags. Features of our
              approach are the development of a hierarchy of related models
              and the estimation of parameter values, along with their
              non-linear biases and standard deviations, for sets of
              experimental data.
              We demonstrate our method of solving parameter estimation
              problems for neutral delay differential equations by analyzing
              some models of cell growth that incorporate a time-lag in the
              cell division phase. We show that these models are more
              consistent with certain reported data than the classic
              exponential growth model. Although the exponential growth model
              provides estimates of some of the growth characteristics, such
              as the population-doubling time, the time-lag growth models can
              additionally provide estimates of: (i) the fraction of cells
              that are dividing, (ii) the rate of commitment of cells to cell
              division, (iii) the initial distribution of cells in the cell
              cycle, and (iv) the degree of synchronization of cells in the
              (initial) cell population.",
  keywords = "Cell proliferation, mathematical modelling, time-lag, neutral
              delay differential equation"
        }

@techreport{Ford:1997:SDA,
  author = "Neville J. Ford and John T. Edwards and Jason A. Roberts and
            Leonid E. Shaikhet",
  title = "Stability of a difference analogue for a nonlinear
           integro differential equation of convolution type",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 312",
  pages = "15",
  month = oct,
  year = "1997",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep312.ps.gz",
  abstract = "This paper is concerned with the stability of discrete
              analogues of nonlinear integro-differential equations.
              We apply the general method of Lyapunov functional
              construction to give sufficient conditions for
              the asymptotic stability of the zero solution",
  keywords = "Integro-differential equations, numerical stability,
              Lyapunov stability"
        }

@techreport{Ford:1997:NVI,
  author = "Neville J. Ford and Christopher T. H. Baker and Jason A. Roberts",
  title = "Nonlinear {Volterra} integro-differential equations -
	   stability and numerical stability of theta methods",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 311",
  pages = "15",
  month = oct,
  year = "1997",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep311.ps.gz",
  abstract = "This paper considers whether existing stability results
	      for a non-linear integro-differential equation are
	      preserved when a solution is found using numerical
	      approximations. We choose to analyse theta quadrature methods
	      and we apply the direct method of Lyapunov for the stability analysis",
  keywords = "Integro-differential equations, numerical stability, Lyapunov stability"
        }

@techreport{Ford:1997:QBA,
  author = "Neville J. Ford and John T. Edwards and Kurt Frischmuth",
  title = "Qualitative behaviour of approximate solutions of some {Volterra}
           integral equations with non-{Lipschitz} nonlinearity",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 310",
  pages = "17",
  month = oct,
  year = "1997",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep310.ps.gz",
  abstract = "In some recent work the  long-term behaviour of numerical
              solutions of the nonlinear convolution equation has been
              considered. In the present paper, we consider examples
              which do not satisfy classical  conditions guaranteeing
              existence and uniqueness  of  the  exact solutions and
              suitable behaviour  of  approximate solutions. We are able
              to give a strengthened  version of a theorem given
              originally by Corduneanu for certain kernels and
              nonlinearities. On that basis  we consider how
              certain numerical codes may fail. We explore and test methods
              known to preserve properties of the solutions in the linear
              case",
  keywords = "Volterra integral equations, non-Lipschitz kernel, numerical
              methods"
        }

@techreport{Ford:1997:ADM,
  author = "John T. Edwards and Jason A. Roberts and Neville J. Ford",
  title = "A comparison of {Adomian's} decomposition method and {Runge Kutta}
           methods for approximate solution of some predator prey
           model equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 309",
  pages = "17",
  month = oct,
  year = "1997",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep309.ps.gz",
  abstract = "Adomian's decomposition method (ADM) has been presented
             for example, by Bellman in 1986 as a significant powerful method
             for solving systems of nonlinear equations. A recent paper by
             Olek indicates that ADM is a powerful tool in the solution
             of some predator prey model equations and that it compares
             favourably with alternative approximate methods (both
             analytical and numerical). In this report we investigate
             whether this property holds more generally",
  keywords = "Adomian's Decomposition Method, Predator Prey equations"
        }

@techreport{Higham:1997:SBF,
  author = "Nicholas J. Higham",
  title = "Stability of Block $\mathrm{LDL^T}$ Factorization of a
           Symmetric Tridiagonal Matrix",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 308",
  month = sep,
  pages = "8",
  year = "1997",
  URL = "http://www.maths.man.ac.uk/~nareports/narep308.ps.gz",
  abstract = "For symmetric indefinite tridiagonal matrices, block
              $\mathrm{LDL^T}$ factorization without interchanges is shown to
              have excellent numerical stability when a pivoting strategy of
              Bunch is used to choose the dimension (1 or 2) of the pivots.",
  keywords = "tridiagonal matrix, symmetric indefinite matrix, diagonal
              pivoting method, $\mathrm{LDL^T}$ factorization, growth factor,
              numerical stability, rounding error analysis, LAPACK,
              LINPACK."
       }

@techreport{Norburn:1997:SSM,
  author = "Sean Norburn and David Silvester",
  title = "Stabilised Vs Stable Mixed Methods for Incompressible Flow",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 307",
  month = sep,
  year = "1997",
  note = "Submitted to Computer Methods in Applied Mechanics and Engineering",
  URL = "http://www.maths.man.ac.uk/~nareports/narep307.ps.gz",
  abstract = "The accuracy of low-order mixed finite element methods for the
              incompressible (Navier-)Stokes equations is investigated in this
              work. Some numerical experiments suggest that the lowest-order
              stabilised $P_1$--$P_0$ method (linear velocity, constant
              pressure) is more efficient than the alternative a-priori stable
              non-conforming Crouzeix-Raviart ($P_{-1}$--$P_0$) approach. The
              relative accuracy of stabilised $P_1$--$P_0$ and $P_1$--$P_1$
              (linear velocity, continuous linear pressure) is also assessed
              herein."
      }

@techreport{Cox:1997:ASN,
  author = "Anthony J. Cox and Nicholas J. Higham",
  title = "Accuracy and Stability of the Null Space Method for Solving the
           Equality Constrained Least Squares Problem",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 306",
  pages = "20",
  month = aug,
  year = "1997",
  URL = "http://www.maths.man.ac.uk/~nareports/narep306.ps.gz",
  abstract = "The null space method is a standard method for solving the linear
              least squares problem subject to equality constraints (the LSE
              problem). We show that three variants of the method, including
              one used in LAPACK that is based on the generalized QR
              factorization, are numerically stable. We derive two perturbation
              bounds for the LSE problem: one of standard form that is not
              attainable, and a bound that yields the condition number of the
              LSE problem to within a small constant factor. By combining the
              backward error analysis and perturbation bounds we derive an
              approximate forward error bound suitable for practical
              computation. Numerical experiments are given to illustrate the
              sharpness of this bound.",
  keywords = "Constrained least squares problem, null space method,
              rounding error analysis, condition number,
              generalized QR factorization, LAPACK"
        }

@techreport{Higham:1997:SIM,
  author = "Nicholas J. Higham",
  title = "Stable Iterations for the Matrix Square Root",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 305",
  month = apr,
  pages = "20",
  year = "1997",
  URL = "http://www.maths.man.ac.uk/~nareports/narep305.pdf",
  abstract = "Any matrix with no nonpositive real eigenvalues has a unique
              square root for which every eigenvalue lies in the open right
              half-plane. A link between the matrix sign function and this
              square root is exploited to derive both old and new iterations
              for the square root from iterations for the sign function. One
              new iteration is a quadratically convergent Schulz iteration
              based entirely on matrix multiplication; it converges only
              locally, but can be used to compute the square root of any
              nonsingular $M$-matrix. A new Pad\'e iteration well suited to
              parallel implementation is also derived and its properties
              explained. Iterative methods for the matrix square root are
              notorious for suffering from numerical instability. It is shown
              that apparently innocuous algorithmic modifications to the Pad\'e
              iteration can lead to instability, and a perturbation analysis is
              given to provide some explanation. Numerical experiments are
              included and some advice is offered on the choice of iterative
              method for computing the matrix square root.",
  keywords = "Matrix square root, matrix logarithm, matrix sign function,
              $M$-matrix, symmetric positive definite matrix, Pad\'e
              approximation, numerical stability, Newton's method, Schulz
              method"
      }

@techreport{MCCM:1997:AR,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 1996",
  type =         "Numerical Analysis Report",
  number =       "No. 304",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  pages = "16",
  year =         "1997",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep304.ps.gz",
}

@techreport{Braconnier:1997:CBE,
  author = "Thierry Braconnier and Fran\c{c}oise Chaitin-Chatelin",
  title = "Chaotic Behavior for Eigensolvers Applied on
           Highly Nonnormal Matrices in Finite Precision",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 303",
  pages = "15",
  month = may,
  year = "1997",
  URL = "http://www.maths.man.ac.uk/~nareports/narep303.ps.gz",
  abstract = "In this paper, we provide examples of nonnormal matrices where
              the power method is shown to be backward unstable. Investigating
              the behavior of the computed sequence of Rayleigh quotients, we
              show via two measures of chaos, i.e, the computation of the
              Hausdorff dimension (more precisely of the correlation dimension)
              and the computation of the Lyapunov exponent, that the computed
              dominant eigenvalue lies on an attractor which has a dimension
              larger than the dimension predicted by the theoretical analysis
              and which increases with the nonnormality of the matrix.",
  keywords = "Backward stability, Hausdorff dimension, correlation dimension,
              Lyapunov exponents, eigenvalues, departure from normality."
   }

@techreport{Williams:1997:LHC,
  author = "Jack Williams and Z. Kalogiratou",
  title = "The Local {Haar} Condition in Parameter
           Estimation for Second Order Ordinary Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 302",
  month = feb,
  pages = "17",
  year = "1997",
  URL = "http://www.maths.man.ac.uk/~nareports/narep302.ps.gz"
        }

@techreport{Cox:1997:SHQ,
  author = "Anthony J. Cox and Nicholas J. Higham",
  title = "Stability of {Householder QR} Factorization
           for Weighted Least Squares Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 301",
  pages = "17",
  month = feb,
  year = "1997",
  URL = "http://www.maths.man.ac.uk/~nareports/narep301.ps.gz",
  abstract = "For least squares problems in which the rows of the coefficient
              matrix vary widely in norm, Householder QR factorization (without
              pivoting) has unsatisfactory backward stability properties.
              Powell and Reid showed in 1969 that the use of both row and
              column pivoting leads to a desirable row-wise backward error
              result. We give a reworked backward error analysis in modern
              notation and prove two new results. First, sorting the rows by
              decreasing $\infty$-norm at the start of the factorization
              obviates the need for row pivoting. Second, row-wise backward
              stability is obtained for only one of the two possible choices of
              sign in the Householder vector.",
  keywords = "Weighted least squares problem, Householder matrix, QR
              factorization, row pivoting, column pivoting, backward error
              analysis"
        }

@techreport{Usman:1996:ESP,
  author = "Anila Usman and George Hall",
  title = "Equilibrium States for Predictor-Corrector Methods",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 300",
  month = dec,
  pages = "35",
  year = "1996",
  URL = "http://www.maths.man.ac.uk/~nareports/narep300.ps.gz"
        }

@techreport{Baker:1996:VCT,
  author = "Christopher T. H. Baker and Arsalang Tang",
  title = "Volterra Centennial Meetings---Invited talks at {Arlington and
           Tempe}",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 299",
  month = nov,
  pages = "28",
  year = "1996",
  URL = "http://www.maths.man.ac.uk/~nareports/narep299.ps.gz"
        }

@techreport{Higham:1996:FCS,
  author = "Nicholas J. Higham",
  title = "Factorizing Complex Symmetric Matrices with
           Positive Definite Real and Imaginary Parts",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 298",
  month = nov,
  pages = "12",
  year = "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep298.ps.gz",
  abstract = "Complex symmetric matrices whose real and imaginary parts are
              positive definite are shown to have a growth factor bounded by 2
              for LU factorization. This result adds to the classes of matrix
              for which it is known to be safe not to pivot in LU
              factorization. Block $\mathrm{LDL^T}$ \fact\ with the pivoting
              strategy of Bunch and Kaufman is also considered, and it is shown
              that for such matrices only $1\times 1$ pivots are used and the
              same growth factor bound of 2 holds, but that interchanges that
              destroy band structure may be made. The latter results hold
              whether the pivoting strategy uses the usual absolute value or
              the modification employed in LINPACK and LAPACK.",
  keywords = "Complex symmetric matrices, LU factorization, diagonal pivoting
              factorization, block $\mathrm{LDL^T}$ factorization,
              Bunch--Kaufman pivoting strategy, growth factor, band matrix,
              LINPACK, LAPACK"
        }

@techreport{Higham:1996:SBE,
  author = "Desmond J. Higham and Nicholas J. Higham",
  title = "Structured Backward Error and Condition of
           Generalized Eigenvalue Problems",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 297",
  month = nov,
  pages = "24",
  year = "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep297.pdf",
  abstract = "Backward errors and condition numbers are defined and evaluated
              for eigenvalues and eigenvectors of generalized eigenvalue
              problems. Both normwise and componentwise measures are used.
              Unstructured problems are considered first, and then the basic
              definitions are extended so that linear structure in the
              coefficient matrices (for example, Hermitian, Toeplitz or band
              structure) is preserved by the perturbations.",
  keywords = "Generalized eigenvalue problem, backward error, condition number,
              structured matrices"
       }

@techreport{Barlow:1996:GGE,
  author = "Jesse L. Barlow and Hongyuan Zha",
  title = "Growth in {Gaussian} Elimination, Orthogonal Matrices, and the
           {Euclidean} Norm",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 296",
  month = sep,
  year = "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep296.ps.gz",
  abstract = "It is shown that maximal growth for Gaussian elimination as
              measured in the Euclidean norm is achieved by orthogonal
              matrices. A new, sharp bound on that growth is given.",
  keywords = "LU factorization, orthogonal invariance, triangular matrix"
        }

@techreport{Higham:1996:MIM,
  author = "Nicholas J. Higham and Sheung Hun Cheng",
  title = "Modifying the Inertia of Matrices Arising in Optimization",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 295",
  month = sep,
  pages = "17",
  year = "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep295.ps.gz",
  abstract = "Applications in constrained optimization (and other areas)
              produce symmetric matrices with a natural block $2\times 2$
              structure. An optimality condition leads to the problem of
              perturbing the $(1,1)$ block of the matrix to achieve a specific
              inertia. We derive a perturbation of minimal norm, for any
              unitarily invariant norm, that increases the number of
              nonnegative eigenvalues by a given amount, and we show how it can
              be computed efficiently given a factorization of the original
              matrix. We also consider an alternative way to satisfy the
              optimality condition based on a projection approach. Theoretical
              tools developed here include an extension of Ostrowski's theorem
              on congruences and some lemmas on inertias of block $2\times 2$
              symmetric matrices.",
  keywords = "Inertia, optimization, nonlinear programming, unitarily
              invariant norm"
        }

@techreport{Higham:1996:TLA,
  author = "Nicholas J. Higham",
  title = "Testing Linear Algebra Software",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 294",
  month = aug,
  pages = "14",
  year = "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep294.ps.gz",
  abstract = "How can we test the correctness of a computer implementation of
              an algorithm such as Gaussian elimination, or the QR algorithm
              for the eigenproblem? This is an important question for program
              libraries such as LAPACK, that are designed to run on a wide
              range of systems. We discuss testing based on verifying known
              backward or forward error properties of the algorithms, with
              particular reference to the test software in LAPACK. Issues
              considered include the choice of bound to verify, computation of
              the backward error, and choice of test matrices. Some examples of
              bugs in widely used linear algebra software are described.",
  keywords = "Mathematical software, testing, linear algebra, backward error,
              LAPACK.",
  note = "To appear in
          {\em The Quality Of Numerical Software: Assessment and Enhancement},
          Ronald F. Boisvert, editor,
          Chapman and Hall, 1997, ISBN 0-412-80530-8"
        }

@techreport{Braconnier:1996:FAF,
  author = "Thierry Braconnier",
  title = "Fvpspack: {A Fortran} and {PVM} Package to Compute the Field of
           Values and Pseudospectra of Large Matrices",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 293",
  month = aug,
  year = "1996",
  pages = "34",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep293.ps.gz",
  abstract = "The field of values and pseudospectra are tools which yield
              insight into the spectral behavior of a matrix. For large sparse
              matrices, both sets can be efficiently computed using a Lanczos
              type method. Since both computations can be done in a natural
              parallel way, we have developed a package including Fortran and
              PVM routines. Experiments show that the PVM codes achieve
              excellent speed-ups and efficiency.",
   keywords = "Field of values, pseudospectra, PVM, parallel algorithm,
               speed-up."
        }

@techreport{Baker:1996:NAV,
  author =       "Christopher T. H. Baker",
  title =        "Numerical Analysis of {Volterra} Functional and Integral
                  Equations---{State} of the Art",
  type =         "Numerical Analysis Report",
  number =       "No. 292",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        sep,
  year =         "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep292.ps.gz",
  abstract =     "Volterra equations are, broadly speaking, evolutionary
                  equations incorporating memory. We indicate where they arise
                  in practice, discuss their significant features (notably
                  smoothness properties of the solutions), describe numerical
                  methods for their approximate solution, and indicate the
                  mathematical foundations of numerical procedures.",
}

@techreport{Smith:1996:IAL,
  author = "Antony Smith and David Silvester",
  title = "Implicit Algorithms and their Linearisation for the
           Transient Incompressible {Navier-Stokes} Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 290",
  month = oct,
  year = "1996",
  note = "Submitted to the IMA Journal of Numerical Analysis",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep290.ps.gz",
  abstract = "The issue of appropriate time discretisation methods for the
              unsteady incompressible Navier-Stokes equations is considered
              from a practical perspective here. Conventional implicit
              time-stepping algorithms are not feasible for long time
              simulations since they inherit the quadratic nonlinearity of the
              steady-state equations. As a result linearised versions of the
              ``pure'' algorithms are analysed herein. These have similar
              stability properties and comparable accuracy to the underlying
              nonlinear methods. "
}

@techreport{Cheng:1996:MCA,
  author = "Sheung Hun Cheng and Nicholas J. Higham",
  title = "A Modified {Cholesky} Algorithm Based on a Symmetric Indefinite
           Factorization",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 289",
  pages = "18",
  month = apr,
  year = "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep289.pdf",
  abstract = "Given a symmetric and not necessarily positive definite matrix
              $A$, a modified Cholesky algorithm computes a Cholesky
              factorization $P(A+E)P^T = R^T R$, where $P$ is a permutation
              matrix and $E$ is a perturbation chosen to make $A+E$ positive
              definite. The aims include producing a suitable small-normed $E$
              and making $A+E$ reasonably well conditioned. Modified Cholesky
              factorizations are widely used in optimization. We propose a new
              modified Cholesky algorithm based on a symmetric indefinite
              factorization computed using a new pivoting strategy of Ashcraft,
              Grimes and Lewis. We analyze the effectiveness of the algorithm,
              both in theory and practice, showing that the algorithm is
              competitive with the existing algorithms of Gill, Murray and
              Wright, and Schnabel and Eskow. Attractive features of the new
              algorithm include easy-to-interpret inequalities that explain the
              extent to which it satisfies its design goals, and the fact that
              it can be implemented in terms of existing software.",
  keywords = "Modified Cholesky factorization, optimization, Newton's method,
              symmetric indefinite factorization"
        }

@techreport{Higham:1996:RDN,
  author = "Nicholas J. Higham",
  title = "Recent Developments in Dense Numerical Linear Algebra",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 288",
  pages = "26",
  month = apr,
  year = "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep288.ps.gz",
  abstract = "We survey recent developments in dense numerical linear algebra,
              covering linear systems, least squares problems and
              eigenproblems. Topics considered include the design and analysis
              of block, partitioned and parallel algorithms, condition number
              estimation, componentwise error analysis, and the computation of
              practical error bounds. Frequent reference is made to LAPACK, the
              state of the art package of Fortran software designed to solve
              linear algebra problems efficiently and accurately on
              high-performance computers.",
  keywords = "Numerical linear algebra, linear system,
              least squares problem, eigenvalue problem, LAPACK",
  note = "Revised August 1996.
          To appear in
          {\em The State of the Art in Numerical Analysis},
          I. S. Duff and G. A. Watson, editors,
          Oxford University Press",
        }

@techreport{Blank:1996:NTD,
  author = "Luise Blank",
  title = "Numerical Treatment of Differential Equations of Fractional Order",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 287",
  pages = "18",
  month = mar,
  year = "1996",
  URL = "http://www.maths.man.ac.uk/~nareports/narep287.ps.gz",
  abstract = "A brief introduction to derivatives of fractional order is given
              and results on the qualitative behaviour of solutions of basic
              linear differential equations of fractional order are cited.
              The collocation approximation with polynomial splines is applied
              to differential equations of fractional order and the systems of
              equations characterizing the numerical solution are determined.
              In particular, the weight matrices resulting from the fractional
              derivative of the spline are deduced and decomposed for
              numerical implementation. The main result of this paper is the
              simplification of the numerical method under a specific
              smoothness condition on the chosen splines. Numerical results
              conclude the work and show that the expected qualitative
              behaviour is given by collocation approximation.",
  keywords = "Differential equations, fractional order, numerical treatment"
       }

@techreport{MCCM:1996:AR,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 1995",
  type =         "Numerical Analysis Report",
  number =       "No. 286",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month = mar,
  pages = "18",
  year =         "1996",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep286.ps.gz",
}

@techreport{Ford:1995:PTV,
  author = "Neville J. Ford and Christopher T. H. Baker",
  title = "Preserving transience in numerical solutions of {Volterra}
           integral equations of convolution type",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 284",
  pages = "13",
  month = nov,
  year = "1995",
  URL = "http://www.maths.man.ac.uk/~nareports/narep284.ps.gz",
  abstract = "We explore the property of transience in solutions of a nonlinear
              integral equation of convolution type under discretisation. We
              show that, under suitable conditions on the kernel function and
              with the choice of a numerical method based on a positive
              quadrature, transient solutions are approximated by transient
              solutions. We show how the analysis can be extended to a wider
              class of kernels",
  keywords = "Volterra integral equations, numerical methods, transient solution, stability
              analysis, positive quadrature"
        }

@techreport{Paul:1995:UGA,
  author = "Christopher A.H. Paul",
  title = "A User-Guide to {Archi}---An Explicit {Runge-Kutta} Code for Solving
           Delay and Neutral Differential Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 283",
  pages = "20",
  month = dec,
  year = "1995",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep283.ps.gz",
  abstract = "Archi is based on the fifth-order Dormand and Prince explicit
              Runge-Kutta method with a fifth-order Hermite interpolant due to
              Shampine. The choice of a Hermite interpolant was influenced by:
              1. The conclusion of Gladwell, Shampine, Baca and Brankin that
              for non-monotonic solutions, quintic Hermite interpolants are
              required in order to preserve the shape of a solution.
              2. In order to maintain the asymptotic correctness of a $p$-th
              order local error estimator, the local error in the approximation
              scheme used to evaluate delay terms must be a least order $p+1$.
              Archi has also been designed to track the propagation of
              discontinuities in a solution using the discontinuity tracking
              theory developed by Wille and Baker, and to solve vanishing-lag
              DDEs using either extrapolation or the predictor-corrector
              approach developed by Baker and Paul. In addition, Archi can
              solve a limited class of delay-integro-differential equations,
              specifically those equations that have a sufficiently smooth
              integrand, using the adaptive quadrature algorithm developed
              by Paul.",
  keywords = "Delay differential equations, numerical solution, user-guide"
        }

@techreport{Blank:1995:STR,
  author = "Luise Blank",
  title = "Stability-Type Results for Collocation Methods for {Volterra}
           Integral Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 282",
  pages = "19",
  month = dec,
  year = "1995",
  URL = "http://www.maths.man.ac.uk/~nareports/narep282.ps.gz",
  abstract = "Stability behaviour is investigated for polynomial spline
              collocation applied to Volterra integral equations with special
              emphasis on a weakly singular test equation. The characterization
              of the corresponding stability domain is related to stability
              results for solutions of finite recursion relations, as they
              occur for methods applied to ordinary differential equations.
              Estimates of the stability domains are established and, as a
              consequence, conditions for $A(\pi/2)$-stability are obtained.
              Exploiting these estimates, collocation approximations with
              constant, linear, and continuous splines are considered in more
              depth and several numerical illustrations are presented.",
  keywords = "Stability, Volterra integral equation, weakly singular kernel,
              collocation method"
       }

@techreport{Braconnier:1995:IOB,
  author = "Thierry Braconnier",
  title = "Influence of Orthogonality on the Backward Error
           and the Stopping Criterion for Krylov Methods",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 281",
  month = dec,
  year = "1995",
  pages = "30",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep281.ps.gz",
  abstract = "Many algorithms for solving linear systems, least squares
              problems or eigenproblems need to compute an orthonormal basis.
              The computation is commonly performed using a QR factorization
              computed using the classical or the modified Gram-Schmidt
              algorithm, the Householder algorithm, the Givens algorithm or the
              Gram-Schmidt algorithm with iterative reorthogonalization. For
              linear systems, it is well-known that the backward stability of
              the process depends crucially on the algorithm used for the QR
              factorization. For the eigenproblem, although textbooks warn
              users about the possible instability of eigensolvers due to loss
              of orthogonality, few theoretical results exist. In this paper we
              show that the loss of orthogonality of the computed basis can
              affect the reliability of the computed eigenpair. We also show
              that the classical residual norm $\norm{Ax-\lambda x}$ and the
              specific computed one using a Krylov method can differ because of
              the loss of orthogonality of the computed basis of the Krylov
              subspace. We give a bound which quantifies the difference between
              both residuals in terms of the loss of orthogonality.",
  keywords = "QR factorization, backward error, eigenvalues, Krylov methods"
        }

@techreport{Braconnier:1995:CFV,
  author = "Thierry Braconnier and Nicholas J. Higham",
  title = "Computing the Field of Values and Pseudospectra
           Using the Lanczos Method with Continuation",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 279",
  month = nov,
  year = "1995",
  pages = "20",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep279.ps.gz",
  abstract = "The field of values and pseudospectra are useful tools for
              understanding the behaviour of various matrix processes. To
              compute these subsets of the complex plane it is necessary to
              estimate one or two eigenvalues of a large number of parametrized
              Hermitian matrices; these computations are prohibitively
              expensive for large, possibly sparse, matrices, if done by use of
              the QR algorithm. We describe an approach based on the Lanczos
              method with selective reorthogonalization and Chebyshev
              acceleration that, when combined with continuation and a shift
              and invert technique, enables efficient and reliable computation
              of the field of values and pseudospectra for large matrices. The
              idea of using the Lanczos method with continuation to compute
              pseudospectra is not new, but in experiments reported here our
              algorithm is faster and more accurate than existing algorithms of
              this type.",
  keywords = "Field of values, pseudospectra, Lanczos method, Chebyshev
              acceleration, continuation."
        }

@techreport{Higham:1995:IRL,
  author = "Nicholas J. Higham",
  title = "Iterative Refinement and {LAPACK}",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 277",
  month = sep,
  year = "1995",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep277.ps.gz",
  abstract = "The technique of iterative refinement for improving the computed
              solution to a linear system was used on desk calculators and
              computers in the 1940s and has remained popular. In the 1990s
              iterative refinement is well supported in software libraries,
              notably in LAPACK. Although the behaviour of iterative refinement
              in floating point arithmetic is reasonably well understood, the
              existing theory is not sufficient to justify the use of fixed
              precision iterative refinement in all the LAPACK routines in
              which it is implemented. We present analysis that provides the
              theoretical support needed for LAPACK. The analysis covers both
              mixed and fixed precision iterative refinement with an arbitrary
              number of iterations, makes only a general assumption on the
              underlying solver, and is relatively short. We identify some
              remaining open problems.",
  keywords = "iterative refinement, rounding error analysis, backward error,
              condition number, LAPACK"
        }

@techreport{Higham:1995:TMT,
  author = "Nicholas J. Higham",
  title = "The {Test Matrix Toolbox} for \textsc{Matlab} (Version 3.0)",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 276",
  month = sep,
  year = "1995",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep276.ps.gz",
  abstract = "We describe version 3.0 of the Test Matrix Toolbox for
              Matlab~4.2. The toolbox contains a collection of test matrices,
              routines for visualizing matrices, routines for direct search
              optimization, and miscellaneous routines that provide useful
              additions to Matlab's existing set of functions. There are 58
              parametrized test matrices, which are mostly square, dense,
              nonrandom, and of arbitrary dimension. The test matrices include
              ones with known inverses or known eigenvalues; ill-conditioned or
              rank deficient matrices; and symmetric, positive definite,
              orthogonal, defective, involutary, and totally positive matrices.
              The visualization routines display surface plots of a matrix and
              its (pseudo-) inverse, the field of values, Gershgorin disks, and
              two- and three-dimensional views of pseudospectra. The direct
              search optimization routines implement the alternating directions
              method, the multidirectional search method and the Nelder--Mead
              simplex method. We explain the need for collections of test
              matrices and summarize the features of the collection in the
              toolbox. We give examples of the use of the toolbox and explain
              some of the interesting properties of the Frank matrix and magic
              square matrices. The leading comment lines from all the toolbox
              routines are listed.",
  keywords = "test matrix, Matlab, pseudospectrum, visualization, Frank matrix,
              magic square matrix, random matrix, direct search optimization"
        }

@techreport{MCCM:1995:AR,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report: {January--December} 1994",
  type =         "Numerical Analysis Report",
  number =       "No. 275",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  year =         "1995",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep275.ps.gz",
}

@techreport{Kalogiratou:1995:BCA,
  author = "Z. Kalogiratou and Jack Williams",
  title = "Best {Chebyshev} Approximation by Families of {ODEs} with a
           Fixed Initial Condition",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 274",
  month = jul,
  year = "1995",
  ftpaddress = "ftp.ma.man.ac.uk",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep274.ps.gz",
  abstract = "Best Chebyshev approximation of real-valued data is treated by
              approximations obtained from solutions of parameter dependent
              initial value problems in ordinary differential equations, in
              which the initial conditions are specified. The problem is a
              modified form of the problem treated by the authors in
              \cite{JH1}, \cite{JH2}, and \cite{JH3}. The fixing of the initial
              condition at $t=a$ requires that approximation be carried out on
              a set which excludes this point. Numerical examples illustrating
              the theory are presented."
}

@techreport{Silvester:1995:FRS,
  author = "David Silvester and Andrew Wathen ",
  title = "Fast and Robust Solvers for Time-Discretised
           Incompressible {Navier-Stokes} Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 273",
  month = jul,
  year = "1995",
  note = "To appear in proceedings of the 16th Dundee Bienniel
          Numerical Analysis Conference",
  ftpaddress = "ftp.ma.man.ac.uk",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep273.ps.gz",
  abstract = "We consider the design of robust and efficient methods for
              solving the Navier-Stokes equations governing laminar flow of a
              viscous incompressible fluid. Three fundamental issues will be
              assessed: the discretisation of the convective transport terms,
              the weak enforcement of the incompressibility constraint in a
              mixed finite element setting, and the solution of the indefinite
              (Stokes-like) systems arising at each time-level. Our aim is to
              prescribe a framework for adaptive error control. The essential
              ingredients are: an unconditionally stable time-discretisation;
              ``natural'' spatial discretisations which are (inf-sup) stable;
              and a fast iterative solution strategy generating iterates
              converging monotonically in an appropriate norm (which mimics the
              dissipation inherent in the continuous system). A distinguishing
              feature of our methodology is the use of multigrid
              preconditioning to accelerate the convergence of our Krylov
              subspace iterative solver. We motivate this with some analysis
              showing that the  contraction rate is bounded away from unity
              independently of the choice of the mixed finite element method
              and the subdivision parameter. Analysis and implementation of
              ``pure'' multigrid methods seems to be relatively complicated and
              (discretisation-) method dependent by comparison. "
     }

@techreport{Williams:1995:EUS,
  author = "Jack Williams",
  title = "Existence and Uniqueness of Solutions
           of the Algebraic Equations in the {BDF} Methods",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 272",
  month = jul,
  year = "1995",
  ftpaddress = "ftp.ma.man.ac.uk",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep272.ps.gz",
  abstract = "The paper presents a new result for the existence and uniqueness
              of solutions of the system of nonlinear algebraic equations
              obtained from the application of the BDF methods to stiff initial
              value problems in ordinary differential equations. Expressed in
              terms of the local truncation error, starting errors in the BDF
              and contractivity properties of the differential equation, the
              result is used to explain how with a certain form of local error
              control existence and uniqueness is guaranteed. Although the form
              of local error control is a theoretical procedure it relates well
              to the error control procedures used in practical codes and
              provides insight into how at a typical step the BDF yields
              solution values. A potential difficulty in trying to compute low
              accuracy numerical solutions is also identified."
}

@techreport{Baker:1995:BNS,
  author =       "C. T. H. Baker and C. A. H. Paul and D. R. Wille",
  title =        "A Bibliography on the Numerical Solution of Delay
                  Differential Equations",
  type =         "Numerical Analysis Report",
  number =       "No. 269",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jun,
  year =         "1995",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep269.ps.gz",
  abstract =     "The aim of this bibliography is to provide an introduction
                  to papers and technical reports in the field of delay
                  differential equations and related differential equations.
                  In addition to the title, authors and reference of an
                  article, we provide the abstract which, if the article has
                  previously appeared as a technical report, comes from the
                  published paper unless indicated by a `+'. The main interest
                  in this bibliography derives from the references to early
                  papers and technical reports in the field, as nowadays
                  on-line search facilities (such as BIDS in the U.K.) provide
                  access to the most recent publications. Although it is hoped
                  to keep the bibliography up-to-date, most immediate effort
                  will be invested in extending the collection of earlier
                  references -- as it is far from being exhaustive. The
                  up-to-date bibliography is only available by anonymous ftp,
                  since it is hoped to update it regularly.",
  keywords =     "Delay differential equations, numerical solution,
                  stability, applications"
}

@techreport{Hall:1995:ESM,
  author =       "George Hall",
  title =        "Equilibrium States for Multistep Methods",
  type =         "Numerical Analysis Report",
  number =       "No. 268",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jun,
  year =         "1995",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep268.ps.gz",
  abstract =     "When the stepsize in non-stiff ODE codes is restricted by
                  stability, an uneven pattern of stepsizes with many step
                  rejections is frequently observed. Results analysing this
                  behaviour have been obtained for Runge-Kutta methods, leading
                  to several papers attempting to improve stepsize control. It
                  is shown here that a similar analysis can be carried out for
                  multistep methods. The explicit Adams 2-step method is used
                  to illustrate the techniques required.",
  keywords =     "Multistep methods, stepsize selection, stability"
}

@techreport{Baker:1995:PPE,
  author =       "Christopher T. H. Baker and Christopher A. H. Paul",
  title =        "Pitfalls in Parameter Estimation for Delay Differential
                  Equations",
  type =         "Numerical Analysis Report",
  number =       "No. 267",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        may,
  year =         "1995",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep267.ps.gz",
  abstract =     "Given a set of data
                  $\{U(\gamma_{i}) \approx u(\gamma_{i};p^{*})\}$
                  corresponding to the delay differential equation (DDE)
                  $u'(t;p) = f(t,u(t;p),u(\alpha(t;p);p);p)$ for
                  $t \ge t_{0}(p)$, the basic problem addressed here is that
                  of calculating the vector $p^{*}$. (We also consider neutral
                  differential equations.) Most approaches to parameter
                  estimation calculate $p^{*}$ by minimizing a suitable
                  objective function that is assumed to be sufficiently smooth.
                  In this paper, we use the derivative discontinuity tracking
                  theory for DDEs to analyze how jumps can arise in the
                  derivatives of a natural objective function. These jumps can
                  occur when estimating parameters in lag functions and
                  estimating the position of the initial point, and as such are
                  not expected to occur in parameter estimation problems for
                  ordinary differential equations.",
  keywords =     "Parameter estimation, delay differential equations,
                  derivative discontinuities"
}

@techreport{Baker:1995:SND,
  author = "C. T. H. Baker and  G. A. Bocharov and C. A. H. Paul and
            A. A. Romanyukha",
  title = "A Short Note on Delay Effects in Cell Proliferation",
  type = "Numerical Analysis Report",
  number = "266",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month = may,
  year = "1995",
  URL = "http://www.maths.man.ac.uk/~nareports/narep266.ps.gz",
  abstract = "We show that a growth model for cell proliferation that
              incorporates a time-lag in the cell division phase is more
              consistent with certain reported data than the classical
              exponential growth model. Although both models provide
              estimates of some of the growth characteristics, such as
              the population doubling-time, the time-lag growth model
              additionally provides estimates of: $(i)$ the cell-doubling
              time, $(ii)$ the fraction of the cells that are dividing,
              $(iii)$ the rate of commitment to cell division and, $(iv)$
              the initial distribution of cells in the cell cycle.",
  keywords = "Parameter estimation, delay differential equations,
              mathematical modelling"
        }

@techreport{Higham:1995:SDP,
  author = "Nicholas J. Higham",
  title = "Stability of the Diagonal Pivoting Method with Partial Pivoting",
  institution = inst-MCCM,
  address =     inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 265",
  month = jul,
  year = "1995",
  URL = "http://www.maths.man.ac.uk/~nareports/narep265.ps.gz",
  abstract = "LAPACK and LINPACK both solve symmetric indefinite linear systems
              using the diagonal pivoting method with the partial pivoting
              strategy of Bunch and Kaufman (1977). No proof of the stability
              of this method has appeared in the literature. It is tempting to
              argue that the diagonal pivoting method is stable for a given
              pivoting strategy if the growth factor is small. We show that
              this argument is false in general, and give a sufficient
              condition for stability. This condition is not satisfied by the
              partial pivoting strategy, because the multipliers are unbounded.
              Nevertheless, using a more specific approach we are able to prove
              the stability of partial pivoting, thereby filling a gap in the
              body of theory supporting LAPACK and LINPACK.",
  keywords = "symmetric indefinite matrix, diagonal pivoting method,
              $\mathrm{LDL^T}$ factorization, partial pivoting, growth factor,
              numerical stability, rounding error analysis, LAPACK, LINPACK."
        }

@techreport{Elman:1995:FNI,
  author = "Howard Elman and David Silvester",
  title = "Fast Nonsymmetric Iterations and Preconditioning for
           {Navier-Stokes} Equations",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 263",
  month = jun,
  year = "1995",
  note = "To appear in SIAM Journal of Scientific Computing",
  keywords = "Preconditioning, Krylov subspace iterations,
              Navier-Stokes equations, convection-diffusion problems",
  URL =      "http://www.maths.man.ac.uk/~nareports/narep263.ps.gz",
  abstract = "Discretization and linearization of the steady-state
              Navier-Stokes equations gives rise to a nonsymmetric indefinite
              linear system of equations. In this paper, we introduce
              preconditioning techniques for such systems with the property
              that the eigenvalues of the preconditioned matrices are bounded
              independently of the mesh size used in the discretization. We
              confirm and supplement these analytic results with a series of
              numerical experiments indicating that Krylov subspace iterative
              methods for nonsymmetric systems display rates of convergence
              that are independent of the mesh parameter. In addition, we show
              that preconditioning costs can be kept small by using iterative
              methods for some intermediate steps performed by the
              preconditioner."
     }

@techreport{Silvester:95:SMF,
  author = "David Silvester",
  title = "Stabilised Mixed Finite Element Methods",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 262",
  month = aug,
  year = "1995",
  note = "To appear in Incompressible Flow and the Finite Element
          Method, by P. M. Gresho",
  URL = "http://www.maths.man.ac.uk/~nareports/narep262.ps.gz",
  abstract = "In this article we review a variety of techniques for stabilising
              mixed finite element methods for solving incompressible flow
              problems. Our emphasis is on low order finite element
              approximations, typically based on piecewise linear velocity with
              equal order or piecewise constant pressure. The optimal choice of
              the stabilisation/regularisation parameter is discussed in
              detail. "
  }

@techreport{Ford:1994:QBS,
  author = "Neville J. Ford and Christopher T. H. Baker",
  title = "Qualitative behaviour and stability of solutions of discretised
           non-linear {Volterra} integral equations of convolution type",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  type = "Numerical Analysis Report",
  number = "No. 261",
  pages = "11",
  month = aug,
  year = "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep261.ps.gz",
  abstract = "We provide a discrete analogue of a theorem given by Corduneanu
              for a nonlinear Volterra equation of convolution type. We are
              able to demonstrate the stability behaviour of some simple
              quadrature rules applied to an example equation and show that the
              qualitative behaviour of solutions is preserved under the
              discretisation for certain families of quadrature rule",
  keywords = "Volterra integral equations, numerical methods, stability analysis"
        }

@techreport{Paul:1994:PPC,
  author =       "Christopher A. H. Paul",
  title =        "Performance and Properties of a Class of Parallel Continuous
                  Explicit Runge-Kutta Methods for Ordinary and Delay
                  Differential Equations equations",
  type =         "Numerical Analysis Report",
  number =       "No. 260",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        oct,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep260.ps.gz",
  abstract =     "In this paper we examine the parallel implementation of the
                  iterated continuous extensions (ICEs) of Baker \& Paul. We
                  indicate how to construct arbitrarily high-order ICEs, and
                  discuss some of the strategies for choosing the `free'
                  parameters of the methods. The numerical results that we
                  present suggest that ICEs implemented in parallel provide an
                  easy and efficient way of writing dense-output ordinary
                  differential equation codes and delay differential equation
                  codes, even in the case that the lag vanishes and/or is
                  state-dependent.",
  keywords =     "Ordinary differential equations, delay differential
                  equations, parallel algorithms, numerical performance,
                  stability"
}

@techreport{Griffiths:1994:UME,
  author =       "D. F. Griffiths and D. J. Silvester",
  title =        "Unstable modes of the {$Q_1$--$P_0$} element",
  type =         "Numerical Analysis Report",
  number =       "No. 257",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        sep,
  year =         "1994",
  note =         "Submitted to SIAM J. Numer. Anal.",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep257.ps.gz",
  abstract =     "In this paper the unstable eigenmodes of $Q_1$--$P_0$
                 velocity/pressure finite element approximation for
                 incompressible flow problems are characterised. It is
                 shown that the inf-sup stability constant is $O(h)$ in
                 two dimensions and $O(h^2)$ in three dimensions. The
                 basic tool for our analysis is the method of modified
                 equations which is applied to finite difference
                 representations of the underlying finite element
                 equations. The asymptotic estimates are confirmed and
                 supplemented by numerical experiments.",
  keywords =     "mixed finite elements, inf-sup stability,
                 incompressible flow"
}

@techreport{Baxter:1994:NEI,
  author =       "B. J. C. Baxter and C. A. Micchelli",
  title =        "Norm Estimates for the {$l^2$}-Inverses of
                 Multivariate {Toeplitz} Matrices",
  type =         "Numerical Analysis Report",
  number =       "No. 255",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jul,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep255.ps.gz",
  abstract =     "We bound the $l^2$-norms of inverses of certain
                 Toeplitz matrices arising from P\'olya frequency
                 functions.",
  keywords =     "Radial basis functions, Toeplitz forms, Polya
                 frequency functions",
}

@techreport{Baxter:1994:RNR,
  author =       "B. J. C. Baxter and N. Sivakumar and J. D. Ward",
  title =        "Regarding the {$p$}-Norms of Radial Basis
                 Interpolation Matrices",
  type =         "Numerical Analysis Report",
  number =       "No. 254",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jul,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep254.ps.gz",
  abstract =     "Radial basis functions are linear combinations of
                 translates of a given radially symmetric function. In
                 this paper we provide some upper bounds on
                 $\|A^{-1}\|_p$ for $p \ge 1$, where $A = (\phi(x_j -
                 x_k))_{j,k=1}^n$ and $\phi$ is the radial basis
                 function. We also study some new applications of total
                 positivity in this context.",
  keywords =     "Radial basis functions, p-norms, total positivity",
}

@techreport{Wille:1994:ESC,
  author =       "David Richard Will\'e",
  title =        "Experiments in Stepsize Control for {Adams} Linear
  		  Multistep Methods",
  type =         "Numerical Analysis Report",
  number =       "No. 253",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        oct,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep253.ps.gz",
  abstract =     "In this paper we consider stepsize selection in one class
  		 of Adams linear multistep methods for ordinary differential
  		 equations. In particular, the exact form of the local error
  		 for a variable step method is considered and a class of new
  		 direct approximations proposed. The implications of this
  		 approach are then discussed and illustrations provided with
  		 numerical results.  Intended applications include differential
  		 equations with derivative discontinuities and initial
  		 stepsize strategies.",
  keywords =     "Adams methods, linear multistep methods, error estimators,
  		 stepsize control",
}

@techreport{Gladwell:1994:IPS,
  author =       "I. Gladwell and J. Wang and R. M. Thomas",
  title =        "Iterations and Predictors for Second Order Problems",
  type =         "Numerical Analysis Report",
  number =       "No. 252",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jul,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep252.ps.gz",
  abstract =     "This paper is concerned with the application of
                 implicit Runge-Kutta methods to initial value problems
                 consisting of second order ordinary differential
                 equations without first derivative terms. In
                 particular, we concentrate on implementing modified
                 Newton iterations and providing predictors for these
                 iterations.",
  keywords =     "Ordinary differential equations, initial value
                 problems, implicit Runge-Kutta methods, predictors,
                 iteration schemes",
}

@techreport{MCCM:1994:AR,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report 1993",
  type =         "Numerical Analysis Report",
  number =       "No. 251",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        may,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep251.ps.gz",
}

@techreport{Thomas:1994:FNT,
  author =       "R. M. Thomas and T. E. Simos and G. V. Mitsou",
  title =        "A Family of {Numerov}-Type Exponentially Fitted
                 Predictor-Corrector Methods for the Numerical
                 Integration of the Radial {Schr{\"o}dinger} Equation",
  type =         "Numerical Analysis Report",
  number =       "No. 249",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jul,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep249.ps.gz",
  abstract =     "A family of predictor-corrector exponential
                 Numerov-type methods is developed for the numerical
                 integration of the one-dimensional Schr{\"o}dinger
                 equation. The formula considered contains certain free
                 parameters which allow it to be fitted automatically to
                 exponential functions. The new methods are very simple
                 and integrate more exponential functions than both the
                 well known fourth order Numerov type exponentially
                 fitted methods and the sixth algebraic order
                 Runge-Kutta type methods. Numerical results also
                 indicate that the new methods are much more accurate
                 than the other exponentially fitted methods mentioned
                 above.",
  keywords =     "Schr{\"o}dinger equation, predictor-corrector methods,
                 Numerov-type methods, exponentially-fitted methods,
                 resonance problem",
}

@techreport{Baker:1994:INS,
  author =       "Christopher T. H. Baker and Christopher A. H. Paul and
                 David R. Wille",
  title =        "Issues in the Numerical Solution of Evolutionary Delay
                 Differential Equations",
  type =         "Numerical Analysis Report",
  number =       "No. 248",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        apr,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep248.ps.gz",
  abstract =     "Delay differential equations are of sufficient
                 importance in modelling real-life phenomena to merit
                 the attention of numerical analysts. In this paper, we
                 discuss key features of delay differential equations
                 and consider the main issues to be addressed when
                 constructing robust numerical codes for their solution.
                 We provide an introduction to the existing literature,
                 and in particular we indicate the approaches adopted by
                 the authors at Manchester.",
  keywords =     "Numerical solution, delay differential equations",
}

@techreport{Wille:1994:NSE,
  author =       "David R. Wille",
  title =        "New Stepsize Estimators for Linear Multistep Methods",
  type =         "Numerical Analysis Report",
  number =       "No. 247",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        mar,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep247.ps.gz",
  abstract =     "The selection of appropriate stepsizes for linear
                 multistep methods requires in general the solution, or
                 approximation, of a non-trivial polynomial root-finding
                 problem. Here we show how this can be reduced to the
                 integration of a simple ODE. Results are presented for
                 Adams-Bashforth-Moulton and backward differentiation
                 formulae using both `error-per-step' and
                 `error-per-unit-step' controls.",
  keywords =     "Linear multistep methods, error estimators, stepsize
                 control",
}

@techreport{Hanby:1994:CCS,
  author =       "R. F. Hanby and D. J. Silvester and J. W. Chew",
  title =        "A Comparison of Coupled and Segregated Iterative
                 Solution Techniques for Incompressible Swirling Flow",
  type =         "Numerical Analysis Report",
  number =       "No. 246",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        mar,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep246.ps.gz",
  abstract =     "In many popular solution algorithms for the
                 incompressible Navier-Stokes equations the coupling
                 between the momentum equations is neglected when the
                 linearized momentum equations are solved to update the
                 velocities. This is known to lead to poor convergence
                 in highly swirling flows where coupling between the
                 radial and tangential momentum equations is strong.
                 Here we propose a coupled solution algorithm in which
                 the linearized momentum and continuity equations are
                 solved simultaneously. Comparisons between the new
                 method and the well-known SIMPLEC method are
                 presented.",
  keywords =     "Navier-Stokes equations; finite differences;
                 unsymmetric linear systems; Krylov subspace methods",
}

@techreport{Hall:1994:NSS,
  author =       "George Hall",
  title =        "A New Stepsize Strategy for {Runge-Kutta} codes",
  type =         "Numerical Analysis Report",
  number =       "No. 245",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        mar,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep245.ps.gz",
  abstract =     "When the stepsize in Runge-Kutta codes is restricted
                 by stability, an uneven pattern of stepsizes with many
                 step rejections is frequently observed. A modified
                 strategy is proposed to smooth out this type of
                 behaviour. Several new estimates for the dominant
                 eigenvalue of the Jacobian are derived. It is shown
                 that such estimates can be used, in a strictly
                 controlled way, to improve the stepsize strategy. Some
                 numerical evidence is presented to show that the
                 modified strategy is effective on a set of widely used
                 test problems.",
  keywords =     "Runge-Kutta, stepsize selection, stability",
}

@techreport{Paul:1994:TSF,
  author =       "Christopher A. H. Paul",
  title =        "A Test Set of Functional Differential Equations",
  type =         "Numerical Analysis Report",
  number =       "No. 243",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        feb,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep243.ps.gz",
  abstract =     "This document is a collection of delay, neutral and
                 functional differential equations (mostly with
                 analytical solutions) taken from the mathematical
                 literature. Also included is the original source of the
                 equation, and other useful information. Further test
                 equations are welcomed, and will appear in the
                 anonymous ftp version of this document. The reference
                 solutions quoted were obtained using an explicit
                 fifth-order continuous Runge-Kutta method based on the
                 fifth-order Dormand and Prince method with a Hermite
                 interpolant.",
  keywords =     "Test set, functional differential equations",
}

@techreport{Papadimitriou:1993:KNA,
  author =       "Pythagoras Papadimitriou",
  title =        "The {KSR1}---{A} Numerical Analyst's Perspective",
  type =         "Numerical Analysis Report",
  number =       "No. 242",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        dec,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep242.ps.gz",
  abstract =     "The Kendall Square Research KSR1 is a virtual shared
                 memory MIMD computer. We give a description of the KSR1
                 aimed at a numerical analyst who wishes to use the KSR1
                 for research. The basic architecture of the KSR1 is
                 described. Parallel constructs and language extensions
                 provided in KSR Fortran are discussed, and numerical
                 software issues are also considered. We describe our
                 practical experiences in using the KSR1 and draw some
                 conclusions.",
  keywords =     "Kendall Square Research KSR1, Virtual shared memory,
                 ALLCACHE memory system, KSR Fortran, parallel
                 constructs, LAPACK, BLAS",
}

@techreport{Higham:1994:SCP,
  author =       "Nicholas J. Higham",
  title =        "A Survey of Componentwise Perturbation Theory in
                 Numerical Linear Algebra",
  type =         "Numerical Analysis Report",
  number =       "No. 241",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        feb,
  year =         "1994",
  note =         "In W.~Gautschi, editor, {\em {Mathematics of Computation}
                  1943--1993: {A} Half Century of Computational Mathematics},
                  volume~48 of {\em Proceedings of Symposia in Applied
                  Mathematics}, pages 49--77. American Mathematical Society,
                  Providence, RI, USA, 1994.",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep241.ps.gz",
  abstract =     "Perturbation bounds in numerical linear algebra are
                 traditionally derived and expressed using norms. Norm
                 bounds cannot reflect the scaling or sparsity of a
                 problem and its perturbation, and so can be unduly
                 weak. If the problem data and its perturbation are
                 measured componentwise, much smaller and more revealing
                 bounds can be obtained. A survey is given of
                 componentwise perturbation theory in numerical linear
                 algebra, covering linear systems, the matrix inverse,
                 matrix factorizations, the least squares problem, and
                 the eigenvalue and singular value problems. Most of the
                 results described have been published in the last five
                 years.",
  keywords =     "Componentwise perturbation bounds, componentwise
                 backward error, linear systems, matrix inverse, matrix
                 factorizations, least squares problem, underdetermined
                 system, eigenvalue problem, singular value problem",
}

@techreport{Higham:1993:PSV,
  author =       "Nicholas J. Higham and Pythagoras Papadimitriou",
  title =        "Parallel Singular Value Decomposition via the Polar
                 Decomposition",
  type =         "Numerical Analysis Report",
  number =       "No. 239",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        oct,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep239.ps.gz",
  abstract =     "A new method is described for computing the singular
                 value decomposition (SVD). It begins by computing the
                 polar decomposition and then computes the spectral
                 decomposition of the Hermitian polar factor. The method
                 is particularly attractive for shared memory parallel
                 computers with a relatively small number of processors,
                 because the polar decomposition can be computed
                 efficiently on such machines using an iterative method
                 developed recently by the authors. This iterative polar
                 decomposition method requires only matrix
                 multiplication and matrix inversion kernels for its
                 implementation and is designed for full rank matrices;
                 thus the proposed SVD method is intended for matrices
                 that are not too close to being rank-deficient. On the
                 Kendall Square KSR1 virtual shared memory computer the
                 new method is up to six times faster than a
                 parallelized version of the LAPACK SVD routine,
                 depending on the condition number of the matrix.",
  keywords =     "singular value decomposition, polar decomposition,
                 numerical stability, LAPACK, level 3 BLAS, Kendall
                 Square Research KSR1 computer",
}

@techreport{Wille:1994:SID,
  author =       "David R. Wille and Christopher T. H. Baker",
  title =        "Some Issues in the Detection and Location of
                 Derivative Discontinuities in Delay Differential
                 Equations",
  type =         "Numerical Analysis Report",
  number =       "No. 238",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        mar,
  year =         "1994",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep238.ps.gz",
  abstract =     "The detection and location of derivative
                 discontinuities is a central issue in the design of
                 robust solvers for delay differential problems. In this
                 paper we discuss a number of particular features
                 associated with one common discontinuity tracking
                 strategy before addressing a particular problem arising
                 for a class of state-dependent problems.",
  keywords =     "Delay differential equations, derivative
                 discontinuities",
}

@techreport{Higham:1993:TMT,
  author =       "Nicholas J. Higham",
  title =        "The {Test Matrix Toolbox} for {MATLAB}",
  type =         "Numerical Analysis Report",
  number =       "No. 237",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        dec,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep237.ps.gz",
  abstract-1 =   "We describe version 2.0 of the Test Matrix Toolbox for
                 Matlab 4. The toolbox contains a collection of test
                 matrices, routines for visualizing matrices, and
                 miscellaneous routines that provide useful additions to
                 Matlab's existing set of functions. There are 58
                 parametrized test matrices, which are mostly square,
                 dense, nonrandom, and of arbitrary dimension. The test
                 matrices include ones with known inverses or known
                 eigenvalues; ill-conditioned or rank deficient
                 matrices; and symmetric, positive definite, orthogonal,
                 defective, involutary, and totally positive matrices.",
  abstract-2 =   "We describe version 2.0 of the Test Matrix Toolbox for
                 Matlab The visualization routines display surface plots
                 of a matrix and its (pseudo-) inverse, the field of
                 values, Gershgorin disks, and two- and
                 three-dimensional views of pseudospectra. We explain
                 the need for collections of test matrices and summarize
                 the features of the collection in the toolbox. We give
                 examples of the use of the toolbox and explain some of
                 the interesting properties of the Frank and Pascal
                 matrices, and of random, magic square and companion
                 matrices. The leading comment lines from all the
                 toolbox routines are listed.",
  keywords =     "test matrix, Matlab, pseudospectrum, visualization,
                 Frank matrix, Pascal matrix, companion matrix, magic
                 square matrix, random matrix",
  note = "*** This report has been withdrawn and is replaced by
              Report No. 276. ***"
}

@techreport{Higham:1993:SPT,
  author =       "Nicholas J. Higham",
  title =        "Stability of Parallel Triangular System Solvers",
  type =         "Numerical Analysis Report",
  number =       "No. 236",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jul,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep236.ps.gz",
  abstract =     "Several parallel algorithms have been proposed for
                 solution of triangular systems. The stability of four
                 of them is analysed here: a fan-in algorithm, a block
                 elimination method, a method based on a factorized
                 power series expansion of the matrix inverse, and a
                 method based on a divide and conquer matrix inversion
                 technique. We derive new forward error and residual
                 bounds, including an improvement on the bounds of Sameh
                 and Brent for the fan-in algorithm. We identify a
                 forward error bound that holds not only for all the
                 methods described here, but for any triangular equation
                 solver that does not rely on algebraic cancellation;
                 among the implications of the bound is that any such
                 method is extremely accurate for certain special types
                 of triangular system.",
  keywords =     "triangular system, matrix inversion, parallel
                 algorithms, fan-in operation, numerical stability,
                 rounding error analysis",
}

@techreport{Higham:1993:MPFb,
  author =       "Nicholas J. Higham and Philip A. Knight",
  title =        "Matrix Powers in Finite Precision Arithmetic",
  type =         "Numerical Analysis Report",
  number =       "No. 235",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        aug,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep235.ps.gz",
  abstract =     "If $A$ is a square matrix with spectral radius less
                 than 1 then $A^k \to 0$ as $k \to \infty$, but the
                 powers computed in finite precision arithmetic may or
                 may not converge. We derive a sufficient condition for
                 $fl(A^k) \to 0$ as $k\to\infty$ and a bound on
                 $\norm{fl(A^k)}$, both expressed in terms of the Jordan
                 canonical form of $A$. Examples show that the results
                 can be sharp. We show that the sufficient condition can
                 be rephrased in terms of a pseudospectrum of $A$ when
                 $A$ is diagonalizable, under certain assumptions. Our
                 analysis leads to the rule of thumb that convergence or
                 divergence of the computed powers of $A$ can be
                 expected according as the spectral radius computed by
                 any backward stable algorithm is less than or greater
                 than 1.",
  keywords =     "matrix powers, rounding errors, Jordan canonical form,
                 nonnormal matrices, pseudospectrum",
}

@techreport{MCCM:1993:AR,
  author =       "Manchester Centre for Computational Mathematics",
  title =        "Annual Report 1992",
  type =         "Numerical Analysis Report",
  number =       "No. 233",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jun,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep233.pdf",
}

@techreport{Paul:1993:FCP,
  author =       "Christopher A. H. Paul",
  title =        "A {FORTRAN} 77 Code for Plotting Stability Regions",
  type =         "Numerical Analysis Report",
  number =       "No. 230",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        may,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep230.ps.gz",
  abstract =     "The standard technique for obtaining the stability
                 regions of numerical methods for ordinary differential
                 equations (ODEs) and delay differential equations
                 (DDEs) is the boundary-locus algorithm (BLA). However,
                 in the case of the DDE $y'(t) = \lambda y(t) + \mu
                 y(t-\tau)$ for $\lambda$ and $\mu$ real, the BLA often
                 fails to map out the stability region correctly. In
                 this paper we give a FORTRAN 77 listing of an
                 alternative stability boundary algorithm, which can be
                 used for both implicit and explicit numerical methods
                 applied to both ODEs and DDEs.",
  keywords =     "stability regions, boundary-locus algorithm,
                 Runge-Kutta methods",
}

@techreport{Baker:1993:GCT,
  author = "Christopher T. H. Baker and Christopher A. H. Paul",
  title = "A Global Convergence Theorem for a Class of Parallel Continuous
           Explicit Runge-Kutta Methods and Vanishing Lag Delay Differential
           Equations",
  institution = "University of Manchester",
  address = "England",
  type = "Numerical Analysis Report",
  number = "No. 229",
  month = may,
  year = "1993",
  keywords = "parallel continuous explicit Runge-Kutta methods, iterated
              continuous extensions, delay differential equations, vanishing
              lag",
  URL = "http://www.maths.man.ac.uk/~nareports/narep229.ps.gz",
  abstract = "Iterated continuous extensions (ICEs) are continuous explicit
              Runge-Kutta methods developed for the numerical solution of
              evolutionary problems in ordinary and delay differential
              equations (DDEs). ICEs have a particular r\^{o}le in the
              explicit solution of DDEs with vanishing lags. They may be
              regarded as parallel continuous explicit Runge-Kutta (PCERK)
              methods, as they allow advantage to be taken of parallel
              architectures. ICEs can also be related to a collocation
              method.
              The purpose of this paper is to provide a theorem giving the
              global order of convergence for variable-step implementations
              of ICEs applied to state-dependent DDEs with and without
              vanishing lags. Implications of the theoryfor the
              implementation of this class of methods are discussed and
              demonstrated. The results establish that our approach allows
              the construction of PCERK methodswhose order of convergence is
              restricted only by the continuity of the solution."
        }

@techreport{Higham:1993:PAC,
  author =       "Nicholas J. Higham and Pythagoras Papadimitriou",
  title =        "A Parallel Algorithm for Computing the Polar
                 Decomposition",
  type =         "Numerical Analysis Report",
  number =       "No. 226",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        jun,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep226.ps.gz",
  abstract =     "The polar decomposition $A=UH$ of a rectangular matrix
                 $A$, where $U$ is unitary and $H$ is Hermitian positive
                 semidefinite, is an important tool in various
                 applications, including aerospace computations, factor
                 analysis and signal processing. We consider a $p$th
                 order iteration for computing $U$ that involves $p$
                 independent matrix inversions per step and which is
                 hence very amenable to parallel computation. We show
                 that scaling the iterates speeds convergence of the
                 iteration but makes the iteration only conditionally
                 stable, with the backward error typically $\kappa_2(A)$
                 times bigger than the unit roundoff. In our
                 implementation of the iteration on the Kendall Square
                 Research KSR1 virtual shared memory MIMD computer we
                 take $p$ to be the number of processors ($p \le 16$ in
                 our experiments). Our code is found to be significantly
                 faster than two existing techniques for computing the
                 polar decomposition: one a Newton iteration, the other
                 based on the singular value decomposition.",
  keywords =     "polar decomposition, singular value decomposition,
                 numerical stability, LAPACK, level 3 BLAS, Kendall
                 Square Research KSR1 computer",
}

@techreport{Higham:1993:MSD,
  author =       "Nicholas J. Higham",
  title =        "The Matrix Sign Decomposition and Its Relation to the
                 Polar Decomposition",
  type =         "Numerical Analysis Report",
  number =       "No. 225",
  institution =  inst-MCCM,
  address =      inst-MCCM:adr,
  month =        apr,
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep225.ps.gz",
  abstract-1 =   "The sign function of a square matrix was introduced by
                 Roberts in 1971. We show that it is useful to regard
                 $S=\sign(A)$ as being part of a matrix sign
                 decomposition $A=SN$, where $N = (A^2)^{1/2}$. This
                 decomposition leads to the new representation $\sign(A)
                 = A(A^2)^{-1/2}$. Most results for the matrix sign
                 decomposition have a counterpart for the polar
                 decomposition $A=UH$, and vice versa.",
  abstract-2 =   "To illustrate this, we derive best approximation
                 properties of the factors $U$, $H$ and~$S$, determine
                 bounds for $\norm{A-S}$ and $\norm{A-U}$, and describe
                 integral formulas for $S$ and $U$. We also derive
                 explicit expressions for the condition numbers of the
                 factors $S$ and $N$. An important equation expresses
                 the sign of a block $2\times2$ matrix involving $A$ in
                 terms of the polar factor $U$ of $A$. We apply this
                 equation to a family of iterations for computing $S$ by
                 Pandey, Kenney and Laub, to obtain a new family of
                 iterations for computing $U$. The iterations have some
                 attractive properties, including suitability for
                 parallel computation.",
  keywords =     "matrix sign function, polar decomposition,
                 conditioning, best approximation, properties, parallel
                 iteration",
}

@techreport{Higham:1992:SPI,
  author =       "Nicholas J. Higham and Alex Pothen",
  title =        "Stability of the Partitioned Inverse Method for
                 Parallel Solution of Sparse Triangular Systems",
  type =         "Numerical Analysis Report",
  number =       "No. 222",
  institution =  inst-U-MANCHESTER,
  address =      inst-U-MANCHESTER:adr,
  month =        oct,
  year =         "1992",
  note =         "To appear in {\em SIAM J. Sci. Comput.}.",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep222.ps.gz",
  abstract =     "Several authors have recently considered a parallel
                 method for solving sparse triangular systems with many
                 right-hand sides. The method employs a partition into
                 sparse factors of the product form of the inverse of
                 the coefficient matrix. It is shown here that while the
                 method can be unstable, stability is guaranteed if a
                 certain scalar that depends on the matrix and the
                 partition is small, and that this scalar is small when
                 the matrix is well-conditioned. Moreover, when the
                 partition is chosen so that the factors have the same
                 sparsity structure as the coefficient matrix, the
                 backward error matrix can be taken to be sparse.",
  keywords =     "sparse matrix, triangular system, substitution
                 algorithm, parallel algorithm, rounding error analysis,
                 matrix inverse",
}

@Article{Higham:1993:FPB,
  author =       "Nicholas J. Higham and Philip A. Knight",
  title =        "Finite Precision Behavior of Stationary Iteration for
                 Solving Singular Systems",
  journal =      "Linear Algebra and Appl.",
  volume =       "192",
  pages =        "165--186",
  year =         "1993",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep215.ps.gz",
  abstract-1 =   "A stationary iterative method for solving a singular
                 system $Ax=b$ converges for any starting vector if
                 $\lim_{i\to\infty}G^i$ exists, where $G$ is the
                 iteration matrix, and the solution to which it
                 converges depends on the starting vector. We examine
                 the behaviour of stationary iteration in finite
                 precision arithmetic. A perturbation bound for singular
                 systems is derived and used to define forward stability
                 of a numerical method. A rounding error analysis
                 enables us to deduce conditions under which a
                 stationary iterative method is forward stable or
                 backward stable.",
  abstract-2 =   "A stationary iterative method for solving a singular
                 system The component of the forward error in the null
                 space of $A$ can grow linearly with the number of
                 iterations, but it is innocuous as long as the
                 iteration converges reasonably quickly. As special
                 cases, we show that when $A$ is symmetric positive
                 semi-definite the Richardson iteration with optimal
                 parameter is forward stable and, if $A$ also has unit
                 diagonal and property A, then the Gauss-Seidel method
                 is both forward and backward stable. Two numerical
                 examples are given to illustrate the analysis.",
  keywords =     "stationary iteration, singular system, Drazin inverse,
                 Jacobi method, Gauss-Seidel method, successive
                 over-relaxation, Richardson iteration, error analysis,
                 numerical stability, Neumann problem",
  mynote =       "Special issue for Proceedings of Workshop on
                 Computational Linear Algebra in Algebraic and Related
                 Problems",
}

@techreport{Paul:1992:FEV,
  author =       "Christopher A. H. Paul",
  title =        "A fast, efficient, very low storage, adaptive
                 quadrature scheme",
  type =         "Numerical Analysis Report",
  number =       "No. 213",
  institution =  inst-U-MANCHESTER,
  address =      inst-U-MANCHESTER:adr,
  month =        may,
  year =         "1992",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep213.ps.gz",
  abstract =     "In this report we review the types of quadrature
                 scheme that are available for the evaluation of
                 integrals. We present a new adaptive quadrature scheme
                 which avoids the storage requirements of normal
                 adaptive schemes. Finally we provide a comparison of
                 all the types of quadrature scheme considered for a set
                 of test integrals.",
  keywords =     "adaptive quadrature, minimal storage",
}

@techreport{Paul:1992:ERK,
  author =       "Christopher A. H. Paul and Christopher T. H. Baker",
  title =        "Explicit {Runge-Kutta} Methods for the Numerical
                 Solution of Singular Delay Differential Equations",
  type =         "Numerical Analysis Report",
  number =       "No. 212",
  institution =  inst-U-MANCHESTER,
  address =      inst-U-MANCHESTER:adr,
  month =        apr,
  year =         "1992",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep212.ps.gz",
  abstract =     "In this paper we are concerned with the development of
                 an explicit Runge-Kutta scheme for the numerical
                 solution of delay differential equations (DDEs) where
                 one or more delay lies in the current Runge-Kutta
                 interval. The scheme presented is also applicable to
                 the numerical solution of {\em Volterra functional
                 equations} (VFEs), although the theory is not covered
                 in this paper. We also derive the stability equations
                 of the scheme for the ODE $y'(t) = \lambda y(t)$, and
                 the DDE $y'(t) = \lambda y(t) + \mu y(t-\tau)$, where
                 the delay $\tau$ and the Runge-Kutta stepsize $H$ are
                 both constant. In the case of the DDE, we consider the
                 two distinct cases: $(i)$ $\tau \ge H$, and $(ii)$
                 $\tau < H$. We evaluate the performance of the scheme
                 by solving several types of {\em singular} DDE and a
                 VFE.",
  keywords =     "explicit Runge-Kutta, singular delay differential
                 equations",
}

@Article{Higham:1993:PTB,
  author = "Nicholas J. Higham",
  title = "Perturbation Theory and Backward Error for {$AX-XB=C$}",
  journal = "BIT",
  volume = 33,
  pages = "124-136",
  year = "1993",
  URL = "http://www.maths.man.ac.uk/~nareports/narep211.ps.gz",
  abstract = "Because of the special structure of the equations $AX-XB=C$ the
              usual relation for linear equations ``$\mbox{backward error} =
              \mbox{relative residual}$'' does not hold, and application of the
              standard perturbation result for $Ax=b$ yields a perturbation
              bound involving ${\rm sep}(A,B)^{-1}$ that is not always
              attainable.  An expression is derived for the backward error of
              an approximate solution $Y$; it shows that the backward error can
              exceed the relative residual by an arbitrary factor.  A sharp
              perturbation bound is derived and it is shown that the condition
              number it defines can be arbitrarily smaller than the ${\rm
              sep}(A,B)^{-1}$-based quantity that is usually used to measure
              sensitivity.  For practical error estimation using the residual
              of a computed solution an ``LAPACK-style'' bound is shown to be
              efficiently computable and potentially much smaller than a
              sep-based bound. A Fortran~77 code has been written that solves
              the Sylvester equation and computes this bound, making use of
              LAPACK routines.",
  keywords = "Sylvester equation, Lyapunov equation, backward error,
              perturbation bound, condition number, error estimate, LAPACK."
}

@techreport{Baker:1992:ESA,
  author =       "Christopher T. H. Baker and John C. Butcher and
                 Christopher A. H. Paul",
  title =        "Experience of {STRIDE} applied to delay differential
                 equations",
  type =         "Numerical Analysis Report",
  number =       "No. 208",
  institution =  inst-U-MANCHESTER,
  address =      inst-U-MANCHESTER:adr,
  month =        feb,
  year =         "1992",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep208.ps.gz",
  abstract =     "STRIDE is intended as a robust adaptive code for
                 solving {\em initial value problems} for {\em ordinary
                 differential equations} (ODEs). The acronym STRIDE
                 stands for {\bf ST}able {\bf R}unge-Kutta {\bf
                 I}ntegrator for {\bf D}ifferential {\bf E}quations. Our
                 purpose here is to report on its adaptation for the
                 numerical solution of a set of delay and neutral
                 differential equations.",
  keywords =     "STRIDE, delay and neutral differential equations",
}

@InProceedings{Higham:1993:MPFa,
  author =       "Nicholas J. Higham and Philip A. Knight",
  editor =       "Carl D. Meyer and Robert J. Plemmons",
  booktitle =    "Linear Algebra, Markov Chains, and Queueing Models",
  title =        "Componentwise Error Analysis for Stationary Iterative
                 Methods",
  volume =       "48",
  publisher =    pub-SV,
  address =      pub-SV:adr,
  pages =        "29--46",
  year =         "1993",
  series =       "IMA Volumes in Mathematics and its Applications",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep206.ps.gz",
  asbtract =     "How small can a stationary iterative method for
                 solving a linear system $Ax=b$ make the error and the
                 residual in the presence of rounding errors? We give a
                 componentwise error analysis that provides an answer to
                 this question and we examine the implications for
                 numerical stability. The Jacobi, Gauss-Seidel and
                 successive over-relaxation methods are all found to be
                 forward stable in a componentwise sense and backward
                 stable in a normwise sense, provided certain conditions
                 are satisfied that involve the matrix, its splitting,
                 and the computed iterates. We show that the stronger
                 property of componentwise backward stability can be
                 achieved using one step of iterative refinement in
                 fixed precision, under suitable assumptions.",
  keywords =     "stationary iteration, Jacobi method, Gauss-Seidel
                 method, successive over-relaxation, error analysis,
                 numerical stability",
}

@techReport{Paul:1991:SBR,
  author =       "Christopher A. H. Paul and Christopher T. H. Baker",
  title =        "Stability boundaries revisited --- {Runge-Kutta}
                 methods for delay differential equations",
  type =         "Numerical Analysis Report",
  number =       "No. 205",
  institution =  inst-U-MANCHESTER,
  address =      inst-U-MANCHESTER:adr,
  month =        dec,
  year =         "1991",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep205.ps.gz",
  abstract =     "We discuss the practical determination of stability
                 regions when various Runge-Kutta methods, combined with
                 continuous extensions, are applied with a fixed
                 stepsize $h$ to the linear delay differential test
                 equation $ y'(t) = \lambda y(t) +\mu y(t-\tau) \quad (t
                 \ge 0)$, with fixed delay $\tau$. The delay $\tau$ is
                 not limited to an integer multiple of the stepsize
                 $h$.",
  keywords =     "delay differential equation, continuous extension,
                 stability",
}

@techReport{Paul:1991:DDD,
  author =       "Christopher A. H. Paul",
  title =        "Developing a Delay Differential Equation Solver",
  type =         "Numerical Analysis Report",
  number =       "No. 204",
  institution =  inst-U-MANCHESTER,
  address =      inst-U-MANCHESTER:adr,
  month =        dec,
  year =         "1991",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep204.ps.gz",
  abstract =     "We discuss briefly various phenomena to be tested when
                 designing a robust code for delay differential
                 equations: choice of interpolant; tracking of
                 discontinuities; vanishing delays; and problems arising
                 from floating-point arithmetic.",
  keywords =     "delay differential equations, interpolation,
                 discontinuities, vanishing delay, real arithmetic",
}

@InProceedings{Higham:1990:HAG,
  author =       "Nicholas J. Higham",
  editor =       "D. F. Griffiths and G. A. Watson",
  booktitle =    "Numerical Analysis 1989, Proceedings of the 13th
                 Dundee Conference",
  title =        "How Accurate is {Gaussian} Elimination?",
  volume =       "228",
  publisher =    pub-LONGMAN,
  address =      pub-LONGMAN:adr,
  pages =        "137--154",
  year =         "1990",
  series =       "Pitman Research Notes in Mathematics",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep173.ps.gz",
  abstract =     "J. H. Wilkinson put Gaussian elimination (GE) on a
                 sound numerical footing in the 1960s when he showed
                 that with partial pivoting the method is stable in the
                 sense of yielding a small backward error. He also
                 derived bounds proportional to the condition number
                 $\kappa(A)$ for the forward error $\|x-\widehat x\|$,
                 where $\widehat x$ is the computed solution to $Ax=b$.
                 More recent work has furthered our understanding of GE,
                 largely through the use of componentwise rather than
                 normwise analysis. We survey what is known about the
                 accuracy of GE in both the forward and backward error
                 senses. Particular topics include: classes of matrix
                 for which it is advantageous not to pivot; how to
                 estimate or compute the backward error; iterative
                 refinement in single precision; and how to compute
                 efficiently a bound on the forward error.",
  keywords =     "Gaussian elimination, partial pivoting, rounding error
                 analysis, backward error, forward error, condition
                 number, iterative refinement in single precision,
                 growth factor, componentwise bounds, condition
                 estimator",
}

@InProceedings{Higham:1990:ACD,
  author =       "Nicholas J. Higham",
  editor =       "M. G. Cox and S. J. Hammarling",
  booktitle =    "Reliable Numerical Computation",
  title =        "Analysis of the {Cholesky} Decomposition of a
                 Semi-definite Matrix",
  publisher =    pub-OUP,
  address =      pub-OUP:adr,
  pages =        "161--185",
  year =         "1990",
  URL =          "http://www.maths.man.ac.uk/~nareports/narep128.ps.gz",
  abstract =     "Perturbation theory is developed for the \Chol\ of an
                 $n \times n$ symmetric \psd\ matrix $A$ of rank~$r$.
                 The matrix $W=\All^{-1}\A{12}$ is found to play a key
                 role in the \pert\ bounds, where $\All$ and $\A{12}$
                 are $r \times r$ and $r \times (n-r)$ submatrices of
                 $A$ respectively. A backward error analysis is given;
                 it shows that the computed Cholesky factors are the
                 exact ones of a matrix whose distance from $A$ is
                 bounded by
                 $4r(r+1)\bigl(\norm{W}+1\bigr)^2u\norm{A}+O(u^2)$,
                 where $u$ is the unit roundoff. For the \cpiv\ strategy
                 it is shown that $\norm{W}^2 \le {1 \over 3}(n-r)(4^r-
                 1)$, and empirical evidence that $\norm{W}$ is usually
                 small is presented. The overall conclusion is that the
                 \Cholalg\ with \cpiv\ is stable for semi-definite
                 matrices. Similar \pert\ results are derived for the
                 \QR\ with \colpiv\ and for the LU decomposition with
                 \cpiv. The results give new insight into the
                 reliability of these decompositions in rank
                 estimation.",
  keywords =     "Cholesky decomposition, positive semi-definite matrix,
                 perturbation theory, backward error analysis, QR
                 decomposition, rank estimation, LINPACK",
}