add subspace methods

note that testing of subspace stuff is cursory for now.

src/GenericSchur.jl

@@ -7,7 +7,7 @@ import LinearAlgebra: lmul!, mul!, checksquare
 # This is the public interface of the package.
 # Wrappers like `schur` and `eigvals` should just work.
 import LinearAlgebra: schur!, eigvals!, eigvecs
-export triangularize
+export triangularize, eigvalscond, subspacesep
 
 schur!(A::StridedMatrix{T}; kwargs...) where {T} = gschur!(A; kwargs...)
 
@@ -662,4 +662,8 @@ end
 include("vectors.jl")
 include("triang.jl")
 
+include("norm1est.jl")
+include("sylvester.jl")
+include("ordschur.jl")
+
 end # module

src/norm1est.jl

@@ -0,0 +1,100 @@
+# This file is part of GenericSchur.jl, released under the MIT "Expat" license
+
+# The method in this file is derived from LAPACK's zlacon.
+# LAPACK is released under a BSD license, and is
+# Copyright:
+# Univ. of Tennessee
+# Univ. of California Berkeley
+# Univ. of Colorado Denver
+# NAG Ltd.
+
+
+# Hager's one-norm estimator, with modifications by N.J. Higham
+"""
+    norm1est!(applyA!,applyAH!,y::Vector) => γ
+
+Estimate the 1-norm of a linear operator `A` expressed as functions which
+apply `A` and `adjoint(A)` to a vector such as `y`.
+
+cf. N.J. Higham, SIAM J. Sci. Stat. Comp. 11, 804 (1990)
+"""
+function norm1est!(applyA!, applyAH!, y::AbstractVector{Ty}; maxiter=5) where Ty
+    n = length(y)
+    x = fill(one(Ty)/n,n)
+    y .= zero(Ty)
+    applyA!(x)
+    (n == 1) && return abs(x[1])
+    γ = norm(x,1)
+    tiny = safemin(real(Ty))
+    for i=1:n
+        absxi = abs(x[i])
+        if absxi > tiny
+            x[i] /= absxi
+        else
+            x[i] = one(Ty)
+        end
+    end
+    applyAH!(x)
+    ax0,j0 = _findamax(x)
+    for iter = 1:maxiter
+        x .= zero(Ty)
+        x[j0] = one(Ty)
+        applyA!(x)
+        y .= x
+        oldγ = γ
+        γ = norm(y,1)
+        if γ <= oldγ
+            break
+        end
+        for i=1:n
+            absxi = abs(x[i])
+            if absxi > tiny
+                x[i] /= absxi
+            else
+                x[i] = one(Ty)
+            end
+        end
+        applyAH!(x)
+        jlast = j0
+        ax0, j0 = _findamax(x)
+        if abs(x[jlast]) == ax0
+            break
+        end
+    end
+    # alternative estimate for tricky cases (see Higham 1990)
+    v = x # reuse workspace
+    isign = 1
+    for i in 1:n
+        v[i] = isign * (1+(i-1)/(n-1))
+        isign = -isign
+    end
+    applyA!(v)
+    t = 2*norm(v,1) / (3*n)
+    return max(t,γ)
+end
+
+function _findamax(x::AbstractVector{T}) where T
+    ax0 = abs(x[1])
+    i0 = 1
+    for i=2:length(x)
+        ax = abs(x[i])
+        if ax > ax0
+            ax0 = ax
+            i0 = i
+        end
+    end
+    return ax0,i0
+end
+
+function _findamax(x::AbstractVector{T}) where T <: Complex
+    ax0 = abs2(x[1])
+    i0 = 1
+    for i=2:length(x)
+        ax = abs2(x[i])
+        if ax > ax0
+            ax0 = ax
+            i0 = i
+        end
+    end
+    return sqrt(ax0),i0
+end

src/ordschur.jl

@@ -0,0 +1,122 @@
+# This file is part of GenericSchur.jl, released under the MIT "Expat" license
+# Portions derived from LAPACK, see below.
+
+# invariant subspace computations using complex Schur decompositions
+
+import LinearAlgebra: _ordschur!, _ordschur
+
+function _ordschur(T::StridedMatrix{Ty},
+                    Z::StridedMatrix{Ty},
+                   select::Union{Vector{Bool},BitVector}) where Ty <: Complex
+    _ordschur!(copy(T), copy(Z), select)
+end
+
+function _ordschur!(T::StridedMatrix{Ty},
+                    Z::StridedMatrix{Ty},
+                    select::Union{Vector{Bool},BitVector}) where Ty <: Complex
+    # suppress most checks since this is an internal function expecting
+    # components of a Schur object
+    n = size(T,1)
+    ks = 0
+    for k=1:n
+        if select[k]
+            ks += 1
+            if k != ks
+                _trexchange!(T,Z,k,ks)
+            end
+        end
+    end
+    triu!(T)
+    T,Z,diag(T)
+end
+
+"""
+reorder `T` by a unitary similarity transformation so that `T[iold,iold]`
+is moved to `T[inew,inew]`. Also right-multiply `Z` by the same transformation.
+"""
+function _trexchange!(T,Z,iold,inew)
+    # this is Algorithm 7.6.1 in Golub & van Loan 4th ed.
+    n = size(T,1)
+    if (n < 1) || (iold == inew)
+        return
+    end
+    if iold < inew
+        krange = iold:inew-1
+    else
+        krange = iold-1:-1:inew
+    end
+    for k in krange
+        Tkk = T[k,k]
+        Tpp = T[k+1,k+1]
+        G,_ = givens(T[k,k+1], Tpp - Tkk, k, k+1)
+        if k+1 <= n
+            lmul!(G,T)
+        end
+        rmul!(T,adjoint(G))
+        Z === nothing || rmul!(Z,adjoint(G))
+    end
+end
+
+# eigvalscond and subspacesep are derived from LAPACK's ztrsen.
+# LAPACK is released under a BSD license, and is
+# Copyright:
+# Univ. of Tennessee
+# Univ. of California Berkeley
+# Univ. of Colorado Denver
+# NAG Ltd.
+
+"""
+    eigvalscond(S::Schur,nsub::Integer)
+
+Estimate the reciprocal of the condition number of the `nsub` leading
+eigenvalues of `S`. (Use `ordschur` to move a subspace of interest
+to the front of `S`.)
+
+See the LAPACK User's Guide for details of interpretation.
+"""
+function eigvalscond(S::Schur{Ty},nsub) where Ty
+    n = size(S.T,1)
+    # solve T₁₁ R - R T₂₂ = σ T₁₂
+    R = S.T[1:nsub,nsub+1:n] # copy
+    R, scale = trsylvester!(view(S.T,1:nsub,1:nsub),view(S.T,nsub+1:n,nsub+1:n),
+                    R)
+    rnorm = norm(R) # Frobenius norm, as desired
+    s = (rnorm > 0) ?
+        scale / (sqrt(scale*scale / rnorm + rnorm ) * sqrt(rnorm)) :
+        one(real(Ty))
+    s
+end
+
+"""
+    subspacesep(S::Schur,nsub::Integer)
+
+Estimate the reciprocal condition of the separation angle for the
+invariant subspace corresponding
+to the leading block of size `nsub` of a Schur decomposition.
+(Use `ordschur` to move a subspace of interest
+to the front of `S`.)
+
+See the LAPACK User's Guide for details of interpretation.
+"""
+function subspacesep(S::Schur{Ty}, nsub) where Ty
+    n = size(S.T,1)
+    scale = one(real(Ty))
+    function f(X)
+        # solve T₁₁ R - R T₂₂ = σ X
+        R, s = trsylvester!(view(S.T,1:nsub,1:nsub),
+                            view(S.T,nsub+1:n,nsub+1:n),
+                            reshape(X,nsub,n-nsub))
+        scale = s
+        R
+    end
+    function fH(X)
+        # solve T₁₁ᴴ R - R T₂₂ᴴ = σ X
+        R, s = adjtrsylvester!(view(S.T,1:nsub,1:nsub),
+                               view(S.T,nsub+1:n,nsub+1:n),
+                               reshape(X,nsub,n-nsub))
+        scale = s
+        R
+    end
+    est = norm1est!(f,fH,zeros(Ty,nsub*(n-nsub)))
+    return scale / est
+end

src/sylvester.jl

@@ -0,0 +1,116 @@
+# This file is part of GenericSchur.jl, released under the MIT "Expat" license
+
+# The methods in this file are derived from LAPACK's ztrsyl.
+# LAPACK is released under a BSD license, and is
+# Copyright:
+# Univ. of Tennessee
+# Univ. of California Berkeley
+# Univ. of Colorado Denver
+# NAG Ltd.
+
+"""
+    trsylvester!(A,B,C) => X, σ
+
+solve the Sylvester equation ``A X - X B = σ C``
+for upper triangular and square `A` and `B`, overwriting `C`
+and setting `σ` to avoid overflow.
+"""
+function trsylvester!(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedVecOrMat{T}; possign=false) where {T}
+    m = checksquare(A)
+    n = checksquare(B)
+    ((size(C,1) == m) && (size(C,2) == n)) || throw(DimensionMismatch(
+        "dimensions of C $(size(C)) must match A, ($m,$m), and B, ($n,$n)"))
+
+    scale = one(real(T))
+
+    tiny = eps(real(T))
+    small = safemin(real(T)) * m * n / tiny
+    bignum = one(real(T)) / small
+    smin = max(small, tiny * norm(A,Inf) * norm(B,Inf))
+
+    isgn = possign ? one(real(T)) : -one(real(T))
+    ierr = 0
+    for l=1:n
+        for k=m:-1:1
+            # FIXME: these should be dotu(), but I can't find anything usable in stdlib::LA
+            # WARNING: this relies on isempty() pass in checkindex() for view()
+            suml = sum(view(A,k,k+1:m) .* view(C,k+1:m,l))
+            sumr = sum(view(C,k,1:l-1) .* view(B,1:l-1,l))
+            v = C[k,l] - (suml + isgn * sumr)
+            scaloc = one(real(T))
+            a11 = A[k,k] + isgn * B[l,l]
+            da11 = abs1(a11)
+            if da11 <= smin
+                a11 = smin
+                da11 = smin
+                ier = 1
+            end
+            db = abs1(v)
+            if (da11 < 1) && (db > 1)
+                if (db > bignum * da11)
+                    # println("scaling by $db")
+                    scaloc = 1 / db
+                end
+            end
+            x11 = (v * scaloc) / a11
+            if scaloc != 1
+                lmul!(scaloc, C)
+                scale *= scaloc
+            end
+            C[k,l] = x11
+        end
+    end
+    return C, scale
+end
+
+"""
+    adjtrsylvesterh!(A,B,C) => X, σ
+
+solve the Sylvester equation ``Aᴴ X - X Bᴴ = σ C``,
+for upper triangular `A` and `B`, overwriting `C`
+and setting `σ` to avoid overflow.
+"""
+function adjtrsylvester!(A::StridedMatrix{T},B::StridedMatrix{T},C::StridedVecOrMat{T}; possign=false) where {T}
+    m = checksquare(A)
+    n = checksquare(B)
+    ((size(C,1) == m) && (size(C,2) == n)) || throw(DimensionMismatch(
+        "dimensions of C $(size(C)) must match A, ($m,$m), and B, ($n,$n)"))
+
+    scale = one(real(T))
+
+    tiny = eps(real(T))
+    small = safemin(real(T)) * m * n / tiny
+    bignum = one(real(T)) / small
+    smin = max(small, tiny * norm(A,Inf) * norm(B,Inf))
+
+    isgn = possign ? one(real(T)) : -one(real(T))
+    for l=n:-1:1
+        for k=1:m
+            suml = dot(view(A,1:k-1,k),view(C,1:k-1,l))
+            # WARNING: this relies on isempty() pass in checkindex() for view()
+            sumr = dot(view(C,k,l+1:n),view(B,l,l+1:n))
+            v = C[k,l] - (suml + isgn * conj(sumr))
+            scaloc = one(real(T))
+            a11 = conj(A[k,k] + isgn * B[l,l])
+            da11 = abs1(a11)
+            if da11 < smin
+                a11 = smin
+                da11 = smin
+                ier = 1
+            end
+            db = abs1(v)
+            if (da11 < 1) && (db > 1)
+                if (db > bignum * da11)
+                    scaloc = 1 / db
+                end
+            end
+            x11 = (v * scaloc) / a11
+            if scaloc != 1
+                lmul!(scaloc, C)
+                scale *= scaloc
+            end
+            C[k,l] = x11
+        end
+    end
+    return C, scale
+end