Commit 9deee46
authored
Bunch-Kaufman factorization support for generic number types and inertia computations (#51487)
### Introduction
This PR adds a generic implementation of the Bunch-Kaufman factorization
in native Julia code, and a generic implementation of a inertia
calculation function. Right now Julia only support the Bunch-Kaufman
factorization for `Float32`, `Float64` and their complex variants. This
is because the factorization is handled by LAPACK, which only supports
these types. To extend support to generic number types, I translated the
LAPACK implementation to native Julia code, and the code performs the
factorization in-place. I also included the function `inertia` that
computes the number of positive, negative, and zero eigenvalues of an $n
\times n$ Bunch-Kaufman factorized matrix in $\mathcal{O}(n)$ time.
### Changes
- `bunchkaufman` and `bunchkaufman!` now work for any `AbstractFloat`,
`Rational` and their complex variants. Behavior for previously supported
types is not changed (LAPACK is used when possible). `bunchakaufman!`
does not support `Integer` types, as in general the factorization lies
in the arithmetic closure of the number type (the rationals for the
integers). On the other hand, `bunchakaufman` supports `Integer` types,
by making an internal conversion to `Rational{BigInt}`.
- `ldiv!` for a `BunchKaufman` factorization has extended support for
generic number types with type stability.
- Previously, Julia extracted the diagonal factor of an $n \times n$
`BunchKaufman` object by making a copy of the matrix and then calling a
LAPACK function (`dsyconvf`, `csyconvf`, etc., depending on the number
type). This function also computes the triangular factor, so it runs in
$\mathcal{O}(n^2)$ time. Now Julia uses a native implementation of the
LAPACK function with small modifications, so it computes the diagonal
factor in $\mathcal{O}(n)$ time, without making a new copy of the
matrix.
- Added the function `inertia` that computes the number of positive,
negative and zero eigenvalues of the diagonal factor of an $n \times n$
`BunchKaufman` object, in case that the matrix is real symmetric or
Hermitian. For complex symmetric matrices, `inertia` only computes the
number of zero eigenvalues of the diagonal factor. `inertia` runs in
$\mathcal{O}(n)$ time and only uses arithmetic and real absolute value
operations. Therefore, `inertia` can be used for matrices of any generic
number type, including `Rational`. In particular, for rational matrices
the output of `inertia` is exact (unless a positive tolerance is
specified).
- Unit tests of the `BunchKaufman` library has been adapted to handle
low precision number types (approximate comparisons with tolerance
`sqrt(eps(Float64))` do not make sense when the floating point type is
`Float16`, for example). The test-set now also runs on the following
types: `Float16, Complex{Float16}, BigFloat, Complex{BigFloat},
Complex{Int}, BigInt, Complex{BigInt}, Rational{BigInt},
Complex{Rational{BigInt}}`. Unit tests for the `inertia` function have
been added too.1 parent 3f4cfc6 commit 9deee46
File tree
3 files changed
+1263
-27
lines changed- stdlib/LinearAlgebra
- src
- test
3 files changed
+1263
-27
lines changed| Original file line number | Diff line number | Diff line change | |
|---|---|---|---|
| |||
115 | 115 | | |
116 | 116 | | |
117 | 117 | | |
| 118 | + | |
118 | 119 | | |
119 | 120 | | |
120 | 121 | | |
| |||
0 commit comments