The Schur decomposition is the workhorse for eigensystem analysis of dense matrices. The diagonal eigen-decomposition of normal (especially Hermitian) matrices is an important special case, but for non-normal matrices the Schur form is often more useful.

The purpose of this package is to extend the `schur!`

and related
functions of the standard library to number types not handled by
LAPACK, such as `Complex{BigFloat}, Complex{Float128}`

(from
Quadmath.jl), etc. For these, the `schur!`

, `eigvals!`

, and `eigen!`

functions in the `LinearAlgebra`

standard library are overloaded here,
and may be accessed through the usual `schur`

, `eigvals`

, and `eigen`

wrappers:

```
A = your_matrix_generator()
S = schur(A)
```

The result `S`

is a `LinearAlgebra.Schur`

object, with the properties `T`

,
`Z=vectors`

, and `values`

.

The unexported `gschur`

and `gschur!`

functions are available for types
normally handled by the LAPACK wrappers in `LinearAlgebra`

.

For square matrices of complex element type, this package provides a full Schur decomposition:

```
A::StridedMatrix{C} where {C <: Complex} == Z * T * adjoint(Z)
```

where `T`

is upper-triangular and `Z`

is unitary, both with the same element
type as `A`

. (See below for real matrices.)

The algorithm is essentially the unblocked, serial, single-shift Francis (QR) scheme used in the complex LAPACK routines. Balancing is also available.

Right and left eigenvectors are available from complex Schur factorizations, using

```
S = schur(A)
VR = eigvecs(S)
VL = eigvecs(S,left=true)
```

As of v0.4, eigenvectors as returned from our `eigen`

and `eigvecs`

methods
for the standard problem have unit Euclidean norm. This accords with the
current (undocumented) behavior of `LinearAlgebra`

methods. (Previously a
convention based on low-level LAPACK routines was used here.)

A quasi-triangular "real Schur" decomposition of real matrices is also provided:

```
A::StridedMatrix{R} where {R <: Real} == Z * T * transpose(Z)
```

where `T`

is quasi-upper-triangular and `Z`

is orthogonal, both with the
same element type as `A`

. This is what you get by invoking the above-mentioned
functions with matrix arguments whose element type `T <: Real`

.
The result is in standard form, so
pair-blocks (and therefore rank-2 invariant subspaces) should be fully resolved.

Eigenvectors are not currently available for the "real Schur" forms.
But don't despair; one can convert a standard quasi-triangular real `Schur`

into a complex `Schur`

with the `triangularize`

function provided here.

The accuracy of eigenvalues and eigenvectors may be improved for some
matrices by use of a similarity transform which reduces the matrix
norm. This is done by default in the `eigen!`

method, and may also be
handled explicitly via the `balance!`

function provided here:

```
Ab, B = balance!(copy(A))
S = schur(Ab)
v = eigvecs(S)
lmul!(B, v) # to get the eigenvectors of A
```

More details are in the function docstring. Although the balancing function also does permutations to isolate trivial subspaces, the Schur routines do not yet exploit this opportunity for reduced workload.

Methods for reordering a Schur decomposition (`ordschur`

) and computing
condition numbers (`eigvalscond`

) and subspace separation (`subspacesep`

)
are provided.
Tests to date suggest that behavior is analogous
to the LAPACK routines on which the implementation is based.

Methods for the generalized eigenvalue problem (matrix pencils) with
`Complex`

element types are available as of release 0.3.0;
in particular, extension of `schur(A,B)`

from LinearAlgebra.
The algorithms are translated from LAPACK, but this implementation has
had limited testing. (Note that it is easy to check the decomposition
of a particular case ex post facto.)

Corresponding functions for reordering and condition estimation are included. Tests to date suggest that behavior is analogous to LAPACK.

Right eigenvectors of generalized problems are available with
`V = eigvecs(S::GeneralizedSchur{<:Complex})`

. Column `j`

of `V`

satisfies
`S.beta[j] * A * v ≈ S.alpha[j] * B * v`

.
These currently have a peculiar norm intended to be compatible with LAPACK
conventions.

This package incorporates or elaborates several methods from Andreas Noack's GenericLinearAlgebra.jl package, and includes translations from LAPACK code.

11/03/2018

4 days ago

63 commits