Skip to content

Commit

Permalink
Merge pull request #8467 from axsk/ordschur
Browse files Browse the repository at this point in the history 
RFC: add ordschur function to base/linalg
  • Loading branch information
@andreasnoack
andreasnoack committed Sep 30, 2014
2 parents 43fcd82 + 7ee06f1 commit 8b88617
Show file tree
Hide file tree
Showing 6 changed files with 136 additions and 1 deletion.
2 changes: 2 additions & 0 deletions base/exports.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -666,6 +666,8 @@ export
lyap,
norm,
null,
ordschur!,
ordschur,
peakflops,
pinv,
qr,
Expand Down
2 changes: 2 additions & 0 deletions base/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,6 +85,8 @@ export
lyap,
norm,
null,
ordschur!,
ordschur,
peakflops,
pinv,
qr,
Expand Down
5 changes: 5 additions & 0 deletions base/linalg/factorization.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -659,6 +659,11 @@ function schur(A::AbstractMatrix)
SchurF[:T], SchurF[:Z], SchurF[:values]
end

ordschur!{Ty<:BlasFloat}(Q::StridedMatrix{Ty}, T::StridedMatrix{Ty}, select::Array{Int}) = Schur(LinAlg.LAPACK.trsen!(select, T , Q)...)
ordschur{Ty<:BlasFloat}(Q::StridedMatrix{Ty}, T::StridedMatrix{Ty}, select::Array{Int}) = ordschur!(copy(Q), copy(T), select)
ordschur!{Ty<:BlasFloat}(schur::Schur{Ty}, select::Array{Int}) = (res=ordschur!(schur.Z, schur.T, select); schur[:values][:]=res[:values]; res)
ordschur{Ty<:BlasFloat}(schur::Schur{Ty}, select::Array{Int}) = ordschur(schur.Z, schur.T, select)

immutable GeneralizedSchur{Ty<:BlasFloat} <: Factorization{Ty}
S::Matrix{Ty}
T::Matrix{Ty}
Expand Down
98 changes: 98 additions & 0 deletions base/linalg/lapack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3508,6 +3508,104 @@ for (gees, gges, elty, relty) in
end
end
end
# Reorder Schur forms
for (trsen, elty) in
((:dtrsen_,:Float64),
(:strsen_,:Float32))
@eval begin
function trsen!(select::Array{Int}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
# * .. Scalar Arguments ..
# CHARACTER COMPQ, JOB
# INTEGER INFO, LDQ, LDT, LIWORK, LWORK, M, N
# DOUBLE PRECISION S, SEP
# * ..
# * .. Array Arguments ..
# LOGICAL SELECT( * )
# INTEGER IWORK( * )
# DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ), WR( * )
chkstride1(T, Q)
n = chksquare(T)
ld = max(1, n)
wr = similar(T, $elty, n)
wi = similar(T, $elty, n)
m = sum(select)
work = Array($elty, 1)
lwork = blas_int(-1)
iwork = Array(BlasInt, 1)
liwork = blas_int(-1)
info = Array(BlasInt, 1)
select = convert(Array{BlasInt}, select)

for i = 1:2
ccall(($(string(trsen)), liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{Void}, Ptr{Void},
Ptr{$elty}, Ptr {BlasInt}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{BlasInt}),
&'N', &'V', select, &n,
T, &ld, Q, &ld,
wr, wi, &m, C_NULL, C_NULL,
work, &lwork, iwork, &liwork,
info)
@lapackerror
if i == 1 # only estimated optimal lwork, liwork
lwork = blas_int(real(work[1]))
liwork = blas_int(real(iwork[1]))
work = Array($elty, lwork)
iwork = Array(BlasInt, liwork)
end
end
T, Q, all(wi .== 0) ? wr : complex(wr, wi)
end
end
end

for (trsen, elty) in
((:ztrsen_,:Complex128),
(:ctrsen_,:Complex64))
@eval begin
function trsen!(select::Array{Int}, T::StridedMatrix{$elty}, Q::StridedMatrix{$elty})
# * .. Scalar Arguments ..
# CHARACTER COMPQ, JOB
# INTEGER INFO, LDQ, LDT, LWORK, M, N
# DOUBLE PRECISION S, SEP
# * ..
# * .. Array Arguments ..
# LOGICAL SELECT( * )
# COMPLEX Q( LDQ, * ), T( LDT, * ), W( * ), WORK( * )
chkstride1(T, Q)
n = chksquare(T)
ld = max(1, n)
w = similar(T, $elty, n)
m = sum(select)
work = Array($elty, 1)
lwork = blas_int(-1)
info = Array(BlasInt, 1)
select = convert(Array{BlasInt}, select)

for i = 1:2
ccall(($(string(trsen)), liblapack), Void,
(Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{Void}, Ptr{Void},
Ptr{$elty}, Ptr {BlasInt},
Ptr{BlasInt}),
&'N', &'V', select, &n,
T, &ld, Q, &ld,
w, &m, C_NULL, C_NULL,
work, &lwork,
info)
@lapackerror
if i == 1 # only estimated optimal lwork, liwork
lwork = blas_int(real(work[1]))
work = Array($elty, lwork)
end
end
T, Q, w
end
end
end

### Rectangular full packed format

Expand Down
18 changes: 17 additions & 1 deletion doc/stdlib/linalg.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -288,12 +288,28 @@ Linear algebra functions in Julia are largely implemented by calling functions f

.. function:: schurfact!(A)

Computer the Schur factorization of ``A``, overwriting ``A`` in the process. See :func:`schurfact`
Computes the Schur factorization of ``A``, overwriting ``A`` in the process. See :func:`schurfact`

.. function:: schur(A) -> Schur[:T], Schur[:Z], Schur[:values]

See :func:`schurfact`

.. function:: ordschur(Q, T, select) -> Schur

Reorders the Schur factorization of a real matrix ``A=Q*T*Q'`` according to the logical array ``select`` returning a Schur object ``F``. The selected eigenvalues appear in the leading diagonal of ``F[:Schur]`` and the the corresponding leading columns of ``F[:vectors]`` form an orthonormal basis of the corresponding right invariant subspace. A complex conjugate pair of eigenvalues must be either both included or excluded via ``select``.

.. function:: ordschur!(Q, T, select) -> Schur

Reorders the Schur factorization of a real matrix ``A=Q*T*Q'``, overwriting ``Q`` and ``T`` in the process. See :func:`ordschur`

.. function:: ordschur(S, select) -> Schur

Reorders the Schur factorization ``S`` of type ``Schur``.

.. function:: ordschur!(S, select) -> Schur

Reorders the Schur factorization ``S`` of type ``Schur``, overwriting ``S`` in the process. See :func:`ordschur`

.. function:: schurfact(A, B) -> GeneralizedSchur

Computes the Generalized Schur (or QZ) factorization of the matrices ``A`` and ``B``. The (quasi) triangular Schur factors can be obtained from the ``Schur`` object ``F`` with ``F[:S]`` and ``F[:T]``, the left unitary/orthogonal Schur vectors can be obtained with ``F[:left]`` or ``F[:Q]`` and the right unitary/orthogonal Schur vectors can be obtained with ``F[:right]`` or ``F[:Z]`` such that ``A=F[:left]*F[:S]*F[:right]'`` and ``B=F[:left]*F[:T]*F[:right]'``. The generalized eigenvalues of ``A`` and ``B`` can be obtained with ``F[:alpha]./F[:beta]``.
Expand Down
12 changes: 12 additions & 0 deletions test/linalg1.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,18 @@ debug && println("Schur")
@test istriu(f[:Schur]) || iseltype(a,Real)
end

debug && println("Reorder Schur")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
# use asym for real schur to enforce tridiag structure
# avoiding partly selection of conj. eigenvalues
ordschura = eltya <: Complex ? a : asym
S = schurfact(ordschura)
select = rand(range(0,2), n)
O = ordschur(S, select)
bool(sum(select)) && @test_approx_eq S[:values][find(select)] O[:values][1:sum(select)]
@test_approx_eq O[:vectors]*O[:Schur]*O[:vectors]' ordschura
end

debug && println("Generalized Schur")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
f = schurfact(a[1:5,1:5], a[6:10,6:10])
Expand Down

0 comments on commit 8b88617

Please sign in to comment.