JuliaLang / JuliaLang/LinearAlgebra.jl
`schurfact!` fails on `StridedMatrix`
- Dominant language
- Julia
- Stars
- 77
- Forks
- 65
- Avg merge
- 3d 23h
- Merged PRs (30d)
- 10
Description
I noticed that the following fails:
``` julia
B = randn(100,100);
A = view(B,1:50,1:50)
schurfact!(A)
```
The reason is that the `Schur` immutable assumes that `T` and `Z` are of the same type:
``` julia
immutable Schur{Ty<:BlasFloat, S<:AbstractMatrix} <: Factorization{Ty}
T::S
Z::S
values::Vector
Schur(T::AbstractMatrix{Ty}, Z::AbstractMatrix{Ty}, values::Vector) = new(T, Z, values)
end
```
which is not the case, because `T` will be stored in the input matrix `A` (in this case a `SubArray`) whereas `Z` will be newly allocated in `LAPACK.gees!` using `similar` on `A`, which produces a regular `Array`.
I could simply fix this by widening the definition of the `Schur` type (similar to how `Eigen` is currently defined), but wanted to consult the linear algebra/lapack experts (my apologies for pinging @andreasnoack , @kshyatt , @tkelman , @jiahao ) to see if that makes sense and whether any other plans or changes are in the pipeline that I could contribute to.
For example, the type of `values` is abstract (`Vector`) so I could also change this (like it is in `Eigen`). However, the output of `schurfact` will still be type unstable because of this line in the output of `LAPACK.gees!`
``` julia
all(wi .== 0) ? wr : complex(wr, wi)
```
Would it be that bad to just always return complex values, with the advantages of having type stability? I noticed the current behaviour also appears in the various function related to `eigfact` etc.
Another question/remark: I noticed that throughout `LAPACK`, new arrays are sometimes allocated with `Array{$elty}` and sometimes with `similar`. Given that the input array is always a `StridedMatrix`, I don't think `similar` can produce anything else than an `Array` (I checked for `SubArray` and for `SharedArray`).
With some input, I would be happy to prepare a PR.
Contributor guide
No contributing guide indexed for this repository
Research direction
Start by reproducing the StridedMatrix example and reading the Schur definition alongside the Eigen definition. Trace LAPACK.gees! to understand its similar-based allocation and the wr/wi value construction, then check the related eigfact behavior. Done requires an agreed, type-correct behavior for schurfact! on views and a clear decision on returned value types.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- julia
- Domain
- hpc
- Issue type
- Bug
- Difficulty
- 5/5
- Estimated time
- Over a week
- Activity status
- Stale
- Clarity
- Needs clarification
- Newbie friendliness
- 32/100