From ed8d0f420d79fb7d9b9093adf744b10f2207a956 Mon Sep 17 00:00:00 2001 From: RalphAS Date: Sun, 6 Jan 2019 22:17:09 -0500 Subject: [PATCH] add subspace methods note that testing of subspace stuff is cursory for now. --- src/GenericSchur.jl | 6 ++- src/norm1est.jl | 100 ++++++++++++++++++++++++++++++++++++ src/ordschur.jl | 122 ++++++++++++++++++++++++++++++++++++++++++++ src/sylvester.jl | 116 +++++++++++++++++++++++++++++++++++++++++ test/ordschur.jl | 71 ++++++++++++++++++++++++++ test/runtests.jl | 3 +- 6 files changed, 416 insertions(+), 2 deletions(-) create mode 100644 src/norm1est.jl create mode 100644 src/ordschur.jl create mode 100644 src/sylvester.jl create mode 100644 test/ordschur.jl diff --git a/src/GenericSchur.jl b/src/GenericSchur.jl index f5f7227..a0e2fdd 100644 --- a/src/GenericSchur.jl +++ b/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 diff --git a/src/norm1est.jl b/src/norm1est.jl new file mode 100644 index 0000000..bbe2963 --- /dev/null +++ b/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 diff --git a/src/ordschur.jl b/src/ordschur.jl new file mode 100644 index 0000000..e64eb05 --- /dev/null +++ b/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 diff --git a/src/sylvester.jl b/src/sylvester.jl new file mode 100644 index 0000000..c58cb79 --- /dev/null +++ b/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 diff --git a/test/ordschur.jl b/test/ordschur.jl new file mode 100644 index 0000000..25d61fe --- /dev/null +++ b/test/ordschur.jl @@ -0,0 +1,71 @@ + +function checkord(A::Matrix{Ty}, tol=10) where {Ty<:Complex} + n = size(A,1) + ulp = eps(real(Ty)) + + S = schur(A) + T2 = copy(S.T) + Z2 = copy(S.Z) + select = fill(true,n) + for i=1:(n>>1) + select[i] = false + end + T2,Z2,v = invoke(LinearAlgebra._ordschur!, + Tuple{StridedMatrix{T}, StridedMatrix{T}, + Union{Vector{Bool},BitVector}} where T <: Complex, + T2,Z2,select) + + # usual tests for Schur + @test all(tril(T2,-1) .== 0) + @test norm(Z2*T2*Z2' - A) / (n * norm(A) * ulp) < tol + @test norm(I - Z2 * Z2') / (n * ulp) < tol + + # make sure we got the ones we asked for + nwanted = count(select) + wwanted = [S.values[j] for j=1:n if select[j]] + errs = zeros(nwanted) + for i=1:nwanted + w = T2[i,i] + errs[i] = minimum(abs.(wwanted .- w)) / (ulp + abs(w)) + end + @test all(errs .< tol) +end + +using LinearAlgebra.LAPACK: trsen!, BlasInt + +function checkecond(A::Matrix{T}, nsub=3) where T + n = size(A,1) + select = fill(false,n) + inds = rand(1:n,nsub) + for i in inds + select[i] = true + end + S = schur(A) + + @assert real(T) <: BlasFloat + T1 = copy(S.T) + Z1 = copy(S.Z) + fselect = BlasInt.([x ? 1 : 0 for x in select]) + T1,Z1,w1,s1,sep1 = trsen!('B','V',fselect,T1,Z1) + + S2 = ordschur(S, select) + s2 = eigvalscond(S2, count(select)) + sep2 = subspacesep(S2, count(select)) + + @test s1 ≈ s2 + @test sep1 ≈ sep2 + +end + +for Ty in [Complex{Float64}] + @testset "ordschur $Ty" begin + n = 32 + A = rand(Ty,n,n) + checkord(A) + end + @testset "subspace condition $Ty" begin + n = 32 + A = rand(Ty,n,n) + checkecond(A) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 6a15966..fed7662 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -31,8 +31,9 @@ function godunov(T) A,vals end - include("real.jl") include("complex.jl") +include("ordschur.jl") + end # module