From 7ee06f1a76c6d90252cea760725163ef3acfabc1 Mon Sep 17 00:00:00 2001 From: Alex Date: Wed, 24 Sep 2014 17:47:33 +0200 Subject: [PATCH] added ordschur to base/linalg --- base/exports.jl | 2 + base/linalg.jl | 2 + base/linalg/factorization.jl | 5 ++ base/linalg/lapack.jl | 98 ++++++++++++++++++++++++++++++++++++ doc/stdlib/linalg.rst | 18 ++++++- test/linalg1.jl | 12 +++++ 6 files changed, 136 insertions(+), 1 deletion(-) diff --git a/base/exports.jl b/base/exports.jl index d46ceecaba98c..80c62050b6116 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -666,6 +666,8 @@ export lyap, norm, null, + ordschur!, + ordschur, peakflops, pinv, qr, diff --git a/base/linalg.jl b/base/linalg.jl index 560e4325acca6..bf29a3915fe47 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -85,6 +85,8 @@ export lyap, norm, null, + ordschur!, + ordschur, peakflops, pinv, qr, diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index 450349b5631bb..3366f08169cce 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -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} diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index 15be3b058c5d9..dad24731dcfd9 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -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 diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 20a50f1f7b339..3a08ea92672fd 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -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]``. diff --git a/test/linalg1.jl b/test/linalg1.jl index d75e903c63116..9280b1f4ff86f 100644 --- a/test/linalg1.jl +++ b/test/linalg1.jl @@ -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])