This is a major rewrite of the packege with following features:

- works with Julia 1.0+
- uses multi-threading
- memory usage is improved

N.B. The latest master is obtained by
`pkg> add Arrowhead#master`

.

Usage:

```
using Arrowhead
# Set the dimension
n=10
# Generate SymArrow
A=GenSymArrow(n,n)
# Compute the eigenvalue decomposition
E,Info=eigen(A)
# Extract eigenvalues
Λ=E.values
# Extract eigenvectors
U=E.vectors
```

For furthaer examples see the file runtests.jl.

The package contains routines for **forward stable** algorithms which compute:

- all eigenvalues and eigenvectors of a real symmetric arrowhead matrices,
- all eigenvalues and eigenvectors of rank-one modifications of diagonal matrices (DPR1), and
- all singular values and singular vectors of half-arrowhead matrices.

The last class of matrices typically appears in SVD updating problems. The algorithms and their analysis are given in the references.

Eigen/singular values are computed in a forward stable manner. Eigen/singular vectors are computed entrywise to almost full accuracy, so they are automatically mutually orthogonal. The algorithms are based on a shift-and-invert approach. Only a single element of the inverse of the shifted matrix eventually needs to be computed with double the working precision.

The package also contains routines for applications:

- divide-and-conquer routine for symmetric tridiagonal eigenvalue problem
- roots of real polynomials with real distinct roots.

The file `arrowhead1.jl`

contains definitions of types
`SymArrow`

(arrowhead) and `SymDPR1`

. Full matrices are accessible
with `Matrix(A)`

.

The file `arrowhead3.jl`

contains routines to generate random symmetric
arrowhead and DPR1 matrices, `GenSymArrow`

and `GenSymDPR1`

, respectively,
three functions `inv()`

which compute various inverses of `SymArrow`

matrices, two functions `bisect()`

which compute outer eigenvalues of
`SymArrow`

and `SymDPR1`

matrices, the main computational function `eigen()`

which
computes the k-th eigenpair of an ordered unreduced `SymArrow`

,
and the driver function `eigen()`

which computes all eigenvalues and
eigenvectors of a `SymArrow`

.

The file `arrowhead4.jl`

contains three functions `inv()`

which compute
various inverses of `SymDPR1`

matrices, the main computational function `eigen()`

which computes the k-th eigenpair of an ordered unreduced `SymDPR1`

,
and the driver function `eigen()`

which computes all eigenvalues and
eigenvectors of a `SymDPR1`

.

The file `arrowhead5.jl`

contains definition of type `HalfArrow`

. The type is
of the form `[Diagonal(A.D) A.z]`

where either
`length(A.z)=length(A.D)`

or `length(A.z)=length(A.D)+1`

, thus giving two possible
forms of the SVD rank one update. The file `arrowhead6.jl`

contains
the function `doubledot()`

, three functions `inv()`

which compute
various inverses of `HalfArrow`

matrices, the main computational function `svd()`

which computes the k-th singular value triplet `u, σ, v`

of an ordered
unreduced `HalfArrow`

, and the driver function `svd()`

which computes all
singular values and vectors of a `HalfArrow`

.

The file `arrowhead7.jl`

contains a simple function `tdc()`

which implements
divide-and-conquer method for `SymTridiagonal`

matrices by spliting the matrix
in two parts and connecting the parts via eigenvalue decomposition of
arrowhead matrix.

The file `arrowhead7.jl`

conatains the function `rootsah()`

which computes the
roots of `Int32`

, `Int64`

, `Float32`

and `Float64`

polynomials with all distinct
real roots. The computation is
forward stable. The program uses `SymArrow`

form of companion matrix in
barycentric coordinates and
the corresponding `eigen()`

function specially designed for this case.
The file also contains three functions `inv()`

. Similarly, the file
`arrowhead8.jl`

conatains the function `rootsah()`

which computes the
roots of `BigInt`

and `BigFloat`

polynomials with all distinct real roots.
The file also contains function `rootsWDK()`

, an implementation of the
Weierstrass-Durand-Kerner polynomial root finding algorithm.

The functions for arrowhead and half-arrowhead matrices were developed and analysed by Jakovcevic Stor, Barlow and Slapnicar (2013) (see also arXiv:1302.7203). The routines for DPR1 matrices are described and analysed in Jakovcevic Stor, Barlow and Slapnicar (2015) (the paper is freely downloadable until Nov 15, 2015, see also arXiv:1405.7537). The polynomial root finder is described and analyzed in Jakovcevic Stor and Slapnicar (2015).

Double the working precision is implemented by using routines by
T. J. Dekker (1971) from the package DoubleDouble
by Simon Byrne. This package is not maintained any more and is deprecated in favor of
DoubleFloats. However,
we are using lightweight and faster implementation in the file `DoubleDouble.jl`

which is now ported to Julia 1.0+.

Highly appreciated help and advice came from Jiahao Chen, Andreas Noack, Jake Bolewski and Simon Byrne.

09/12/2014

5 months ago

118 commits