diff --git a/docs/make.jl b/docs/make.jl index e1768b1..1f5b0a0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,7 +7,7 @@ makedocs(warnonly = true, authors = "Andreas Varga", format = Documenter.HTML(prettyurls = false), pages = [ - "Home" => "index.md", + "Overview" => "index.md", "Library" => [ "lyapunov.md", "plyapunov.md", diff --git a/docs/src/meoperators.md b/docs/src/meoperators.md index aab0dcf..89ba430 100644 --- a/docs/src/meoperators.md +++ b/docs/src/meoperators.md @@ -1,17 +1,16 @@ # Linear Operators Related to Matrix Equation Solvers -## Lyapunov Operators +## Lyapunov and Lyapunov-like Operators ```@docs lyapop invlyapop -tlyapop -hlyapop -tulyapop -hulyapop +lyaplikeop +tulyaplikeop +hulyaplikeop ``` -## Sylvester Operators +## Sylvester and Sylvester-like Operators ```@docs sylvop diff --git a/src/MatrixEquations.jl b/src/MatrixEquations.jl index fb7c524..b5c659e 100644 --- a/src/MatrixEquations.jl +++ b/src/MatrixEquations.jl @@ -28,7 +28,8 @@ export sylvckr, sylvdkr, gsylvkr, sylvsyskr, dsylvsyskr, tsylvckr, hsylvckr, csylvckr, tsylvdkr, hsylvdkr, csylvdkr, tlyapckr, hlyapckr export opnorm1, opnorm1est, oprcondest, opsepest export lyapop, invlyapop, sylvop, invsylvop, sylvsysop, invsylvsysop, trmatop, eliminationop, duplicationop -export tulyapop, hulyapop, tlyapop, hlyapop, gsylvop +export tulyaplikeop, hulyaplikeop, lyaplikeop, gsylvop +export tulyapop, hulyapop include("meutil.jl") include("sylvester.jl") diff --git a/src/iterative_methods.jl b/src/iterative_methods.jl index 84a75f3..27a1c06 100644 --- a/src/iterative_methods.jl +++ b/src/iterative_methods.jl @@ -267,53 +267,72 @@ function ghsylvi(A::Vector{TA}, B::Vector{TB}, C::Vector{TC}, D::Vector{TD}, E:: return reshape(xt,LT.mx,LT.nx), info end """ - tlyapci(A, C; adj = false, abstol, reltol, maxiter) -> (X,info) + tlyapci(A, C, isig = +1; adj = false, abstol, reltol, maxiter) -> (X,info) -Compute for a rectangular `A` and a symmetric `C` a solution `X` of the continuous T-Lyapunov matrix equation +Compute a solution `X` of the continuous T-Lyapunov matrix equation - A*X +transpose(X)*transpose(A) = C if adj = false, + A*X +isig*transpose(X)*transpose(A) = C if adj = false, or - A*transpose(X) + X*transpose(A) = C if adj = true. + A*transpose(X) + isig*X*transpose(A) = C if adj = true, + +where for `isig = 1`, `C` is a symmetric matrix and for `isig = -1`, `C` is a skew-symmetric matrix. For a matrix `A`, a least-squares solution `X` is determined using a conjugate gradient based iterative method applied to a suitably defined T-Lyapunov linear operator `L:X -> Y` such that `L(X) = C` or `norm(L(X) - C)` is minimized. The keyword arguments `abstol` (default: `abstol = 0`) and `reltol` (default: `reltol = sqrt(eps())`) can be used to provide the desired tolerance for the accuracy of the computed solution and the keyword argument `maxiter` can be used to set the maximum number of iterations (default: maxiter = 1000). """ -function tlyapci(A::AbstractMatrix{T}, C::AbstractMatrix{T}; adj = false, abstol = zero(float(real(T))), reltol = sqrt(eps(float(real(T)))), maxiter = 1000) where {T} +function tlyapci(A::AbstractMatrix{T}, C::AbstractMatrix{T}, isig::Int = 1; adj = false, abstol = zero(float(real(T))), reltol = sqrt(eps(float(real(T)))), maxiter = 1000) where {T} m = LinearAlgebra.checksquare(C) ma, n = size(A) ma == m || throw(DimensionMismatch("A and C have incompatible dimensions")) - issymmetric(C) || throw(ArgumentError("C must be symmetric")) - LT = tlyapop(A; adj) - xt, info = cgls(LT,vec(C); abstol, reltol, maxiter) + abs(isig) == 1 || error(" isig must be either 1 or -1") + if isig == 1 + issymmetric(C) || error("C must be symmetric for isig = 1") + # temporary fix to avoid false results for DoubleFloats + # C == transpose(C) || error("C must be symmetric for isig = 1") + else + iszero(C+transpose(C)) || error("C must be skew-symmetric for isig = -1") + end + LT = lyaplikeop(A; adj, isig, htype = false) + xt, info = cgls(LT, vec(C); abstol, reltol, maxiter) info.flag == 1 || @warn "convergence issues: info = $info" return adj ? reshape(xt,m,n) : reshape(xt,n,m), info end """ - hlyapci(A, C; adj = false, abstol, reltol, maxiter) -> (X,info) + hlyapci(A, C, isig = +1; adj = false, abstol, reltol, maxiter) -> (X,info) + -Compute for a rectangular `A` and a hermitian `C` a solution `X` of the continuous H-Lyapunov matrix equation +Compute a solution `X` of the continuous H-Lyapunov matrix equation - A*X +X'*A' = C if adj = false, + A*X + isig*X'*A' = C if adj = false, or - A*X' + X*A' = C if adj = true. + A*X' + isig*X*A' = C if adj = true, + +where for `isig = 1`, `C` is a hermitian matrix and for `isig = -1`, `C` is a skew-hermitian matrix. For a matrix `A`, a least-squares solution `X` is determined using a conjugate gradient based iterative method applied to a suitably defined T-Lyapunov linear operator `L:X -> Y` such that `L(X) = C` or `norm(L(X) - C)` is minimized. The keyword arguments `abstol` (default: `abstol = 0`) and `reltol` (default: `reltol = sqrt(eps())`) can be used to provide the desired tolerance for the accuracy of the computed solution. The keyword argument `maxiter` can be used to set the maximum number of iterations (default: maxiter = 1000). """ -function hlyapci(A::AbstractMatrix{T}, C::AbstractMatrix{T}; adj = false, abstol = zero(float(real(T))), reltol = sqrt(eps(float(real(T)))), maxiter = 1000) where {T} +function hlyapci(A::AbstractMatrix{T}, C::AbstractMatrix{T}, isig::Int = 1; adj = false, abstol = zero(float(real(T))), reltol = sqrt(eps(float(real(T)))), maxiter = 1000) where {T} m = LinearAlgebra.checksquare(C) ma, n = size(A) ma == m || throw(DimensionMismatch("A and C have incompatible dimensions")) - ishermitian(C) || throw(ArgumentError("C must be hermitian")) - LT = hlyapop(A; adj) + abs(isig) == 1 || error(" isig must be either 1 or -1") + if isig == 1 + ishermitian(C) || error("C must be hermitian for isig = 1") + # temporary fix to avoid false results for DoubleFloats + # C == C' || error("C must be hermitian for isig = 1") + else + iszero(C+C') || error("C must be skew-hermitian for isig = -1") + end + LT = lyaplikeop(A; adj, isig, htype = true) xt, info = cgls(LT,vec(C); abstol, reltol, maxiter) info.flag == 1 || @warn "convergence issues: info = $info" return adj ? reshape(xt,m,n) : reshape(xt,n,m), info @@ -324,11 +343,12 @@ end Compute for an upper triangular `U` and a symmetric `Q` an upper triangular solution `X` of the continuous T-Lyapunov matrix equation - U*transpose(X) + X*transpose(U) = Q if adj = false, + transpose(U)*X + transpose(X)*U = Q if adj = false, -or +or + + U*transpose(X) + X*transpose(U) = Q if adj = true. - transpose(U)*X + transpose(X)*U = Q if adj = true. For a `n×n` upper triangular matrix `U`, a least-squares upper-triangular solution `X` is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator `L:X -> Y`, which maps upper triangular matrices `X` @@ -341,7 +361,7 @@ function tulyapci(U::AbstractMatrix{T}, Q::AbstractMatrix{T}; adj = false, abst n == LinearAlgebra.checksquare(Q) || throw(DimensionMismatch("U and Q have incompatible dimensions")) istriu(U) || throw(ArgumentError("U must be upper triangular")) issymmetric(Q) || throw(ArgumentError("Q must be symmetric")) - LT = tulyapop(adj ? U : transpose(U)) + LT = tulyaplikeop(U; adj) xt, info = cgls(LT,triu2vec(Q); abstol, reltol, maxiter) info.flag == 1 || @warn "convergence issues: info = $info" return vec2triu(xt), info @@ -351,11 +371,11 @@ end Compute for an upper triangular `U` and a hermitian `Q` an upper triangular solution `X` of the continuous H-Lyapunov matrix equation - U*X' + X*U' = Q if adj = false, + U'*X + X'*U = Q if adj = false, or - U'*X + X'*U = Q if adj = true. + U*X' + X*U' = Q if adj = true. For a `n×n` upper triangular matrix `U`, a least-squares upper-triangular solution `X` is determined using a conjugate-gradient based iterative method applied to a suitably defined T-Lyapunov linear operator `L:X -> Y`, which maps upper triangular matrices `X` @@ -367,8 +387,8 @@ function hulyapci(U::AbstractMatrix{T}, Q::AbstractMatrix{T}; adj = false, abst n = LinearAlgebra.checksquare(U) n == LinearAlgebra.checksquare(Q) || throw(DimensionMismatch("U and Q have incompatible dimensions")) istriu(U) || throw(ArgumentError("U must be upper triangular")) - ishermitian(Q) || throw(ArgumentError("Q must be symmetric")) - LT = hulyapop(adj ? U : U') + ishermitian(Q) || throw(ArgumentError("Q must be hermitian")) + LT = hulyaplikeop(U; adj) xt, info = cgls(LT,triu2vec(Q); abstol, reltol, maxiter) info.flag == 1 || @warn "convergence issues: info = $info" return vec2triu(xt), info diff --git a/src/lyapunov.jl b/src/lyapunov.jl index 26709ce..b114ce0 100644 --- a/src/lyapunov.jl +++ b/src/lyapunov.jl @@ -2236,13 +2236,13 @@ end """ X = tulyapc!(U, Q; adj = false) -Compute for a nonsingular upper triangular `U` and a symmetric `Q` an upper triangular solution `X` of the continuous T-Lyapunov matrix equation +Compute for an upper triangular `U` and a symmetric `Q` an upper triangular solution `X` of the continuous T-Lyapunov matrix equation - U*transpose(X) + X*transpose(U) = Q if adj = false, + transpose(U)*X + transpose(X)*U = Q if adj = false, or - transpose(U)*X + transpose(X)*U = Q if adj = true. + U*transpose(X) + X*transpose(U) = Q if adj = true. The solution `X` overwrites the matrix `Q`, while `U` is unchanged. @@ -2261,6 +2261,21 @@ function tulyapc!(U::AbstractMatrix{T}, Q::AbstractMatrix{T}; adj = false, absto issymmetric(Q) || throw(ArgumentError("Q must be symmetric")) if !any(iszero.(diag(U))) if adj + # U*transpose(X) + X*transpose(U) = Q + for i in n:-1:1 + for j in n:-1:1 + if i < j + Q[i, j] -= transpose(view(U, i, j:n))*view(Q, j, j:n) + elseif i == j + Q[i, j] /= 2 + else + Q[i, j] = 0 + end + end + Ut = transpose(LinearAlgebra.UpperTriangular(view(U, i:n, i:n))) + LinearAlgebra.rdiv!(view(Q, i:i, i:n), Ut) + end + else # transpose(U)*X + transpose(X)*U = Q for j in 1:n for i in 1:n @@ -2275,24 +2290,9 @@ function tulyapc!(U::AbstractMatrix{T}, Q::AbstractMatrix{T}; adj = false, absto Ut = transpose(LinearAlgebra.UpperTriangular(view(U, 1:j, 1:j))) LinearAlgebra.ldiv!(Ut, view(Q, 1:j, j)) end - else - # U*transpose(X) + X*transpose(U) = Q - for i in n:-1:1 - for j in n:-1:1 - if i < j - Q[i, j] -= transpose(view(U, i, j:n))*view(Q, j, j:n) - elseif i == j - Q[i, j] /= 2 - else - Q[i, j] = 0 - end - end - Ut = transpose(LinearAlgebra.UpperTriangular(view(U, i:n, i:n))) - LinearAlgebra.rdiv!(view(Q, i:i, i:n), Ut) - end end else - LT = tulyapop(adj ? U : transpose(U)) + LT = tulyaplikeop(U; adj) xt, info = MatrixEquations.cgls(LT,triu2vec(Q); abstol, reltol, maxiter) info.flag == 1 || @warn "convergence issues: info = $info" Q[:] = vec2triu(xt) @@ -2302,13 +2302,13 @@ end """ X = hulyapc!(U, Q; adj = false) -Compute for a nonsingular upper triangular `U` and a hermitian `Q` an upper triangular solution `X` of the continuous H-Lyapunov matrix equation +Compute for an upper triangular `U` and a hermitian `Q` an upper triangular solution `X` of the continuous H-Lyapunov matrix equation - U*X' + X*U' = Q if adj = false, + U'*X + X'*U = Q if adj = false, or - U'*X + X'*U = Q if adj = true. + U*X' + X*U' = Q if adj = true. The solution `X` overwrites the matrix `Q`, while `U` is unchanged. @@ -2327,38 +2327,38 @@ function hulyapc!(U::AbstractMatrix{T}, Q::AbstractMatrix{T}; adj = false, absto ishermitian(Q) || throw(ArgumentError("Q must be hermitian")) if !any(iszero.(diag(U))) if adj - # U'*X + X'*U = Q - for j in 1:n - for i in 1:n - if i < j - Q[i, j] -= view(Q, 1:i, i)'*view(U, 1:i, j) - elseif i == j - Q[i, j] /= 2 - else - Q[i, j] = 0 - end - end - Ut = LinearAlgebra.UpperTriangular(view(U, 1:j, 1:j))' - LinearAlgebra.ldiv!(Ut, view(Q, 1:j, j)) - end - else # U*X' + X*U' = Q for i in n:-1:1 - for j in n:-1:1 - if i < j - Q[i, j] -= view(Q, j, j:n)'*view(U, i, j:n) - elseif i == j - Q[i, j] /= 2 - else - Q[i, j] = 0 - end - end - Ut = LinearAlgebra.UpperTriangular(view(U, i:n, i:n))' - LinearAlgebra.rdiv!(view(Q, i:i, i:n), Ut) - end - end + for j in n:-1:1 + if i < j + Q[i, j] -= view(Q, j, j:n)'*view(U, i, j:n) + elseif i == j + Q[i, j] /= 2 + else + Q[i, j] = 0 + end + end + Ut = LinearAlgebra.UpperTriangular(view(U, i:n, i:n))' + LinearAlgebra.rdiv!(view(Q, i:i, i:n), Ut) + end + else + # U'*X + X'*U = Q + for j in 1:n + for i in 1:n + if i < j + Q[i, j] -= view(Q, 1:i, i)'*view(U, 1:i, j) + elseif i == j + Q[i, j] /= 2 + else + Q[i, j] = 0 + end + end + Ut = LinearAlgebra.UpperTriangular(view(U, 1:j, 1:j))' + LinearAlgebra.ldiv!(Ut, view(Q, 1:j, j)) + end + end else - LT = hulyapop(adj ? U : U') + LT = hulyaplikeop(U; adj) xt, info = cgls(LT,triu2vec(Q); abstol, reltol, maxiter) info.flag == 1 || @warn "convergence issues: info = $info" Q[:] = vec2triu(xt) diff --git a/src/meoperators.jl b/src/meoperators.jl index 389abe2..b8627b4 100644 --- a/src/meoperators.jl +++ b/src/meoperators.jl @@ -104,25 +104,6 @@ function LinearMaps._unsafe_mul!(x::AbstractVector, AT::LinearMaps.TransposeMap{ end return x end -# function LinearMaps._unsafe_mul!(x::AbstractVector, AT::LinearMaps.AdjointMap{T,<:eliminationop{T}}, y::AbstractVector) where {T} -# # X = vec2triu(y) -# n = AT.lmap.dim -# ZERO = zero(eltype(y)) -# @inbounds begin -# k = 1 -# for j = 1:n -# for i = 1:j -# x[(i-1)*n + j] = y[k] -# k += 1 -# end -# for i = j+1:n -# x[(i-1)*n + j] = ZERO -# end -# end -# end -# return x -# end - eliminationop(A) = eliminationop{eltype(A)}(LinearAlgebra.checksquare(A)) """ @@ -184,229 +165,231 @@ function LinearMaps._unsafe_mul!(x::AbstractVector, AT::LinearMaps.TransposeMap{ end return x end -# function LinearMaps._unsafe_mul!(x::AbstractVector, AT::LinearMaps.AdjointMap{T,<:duplicationop{T}}, y::AbstractVector) where {T} -# n = AT.lmap.dim -# # Y = reshape(y, n, n) -# # x[:] = triu2vec(Y+transpose(Y)-Diagonal(Y)) -# @inbounds begin -# k = 1 -# for j = 1:n -# for i = 1:j -# #x[k] = i == j ? Y[j,j] : Y[i,j] + Y[j,i] -# x[k] = i == j ? y[(j-1)*n + j] : y[(i-1)*n + j] + y[(j-1)*n + i] -# k += 1 -# end -# end -# end -# return x -# end - duplicationop(A) = duplicationop{eltype(A)}(LinearAlgebra.checksquare(A)) -struct TLyapMap{T,TA <: AbstractMatrix,CD <: ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T} +struct LyapLikeMap{T,TA <: AbstractMatrix,CD <: ContinuousOrDiscrete,TH<:TtypeOrHtype} <: LyapunovMatrixEquationsMaps{T} A::TA adj::Bool - function TLyapMap{T,TA,CD}(A::TA, adj=false) where {T,TA<:AbstractMatrix{T},CD} - return new{T,TA,CD}(A, adj) + isig::Int + function LyapLikeMap{T,TA,CD,TH}(A::TA; isig = 1, adj=false) where {T,TA<:AbstractMatrix{T},CD,TH} + abs(isig) == 1 || throw(ArgumentError("only 1 or -1 values allowed for isig; got isig = $isig")) + return new{T,TA,CD,TH}(A, adj, isig) end end -TLyapMap(A::TA, ::CD = Continuous(), adj::Bool = false) where {T,TA<:AbstractMatrix{T},CD<:ContinuousOrDiscrete} = - TLyapMap{T,TA,CD}(A, adj) -#LinearAlgebra.transpose(L::TLyapunovMap{<:Any,<:Any,CD}) where {CD} = TLyapunovMap(L.U', CD(), !L.adj) +LyapLikeMap(A::TA, ::CD = Continuous(), ::TH = Ttype(); adj::Bool = false, isig::Int = 1) where {T,TA<:AbstractMatrix{T},CD<:ContinuousOrDiscrete,TH<:TtypeOrHtype} = + LyapLikeMap{T,TA,CD,TH}(A; adj, isig) -function Base.size(L::TLyapMap) +function Base.size(L::LyapLikeMap) m, n = size(L.A) return (m*m, m*n) end """ - L = tlyapop(A; adj = false) + L = lyaplikeop(A; isig = 1, adj = false, htype = false) -Define, for an `m x n` matrix `A`, the continuous T-Lyapunov operator `L:X -> A*X+transpose(X)*transpose(A)` -if `adj = false` and `L:X -> A*transpose(X)+X*transpose(A)` if `adj = true`. +For a matrix `A`, define for `adj = false` the continuous T-Lyapunov operator `L:X -> A*X+isig*transpose(X)*transpose(A)` if `htype = false` +or the continuous H-Lyapunov operator `L:X -> A*X+isig*X'*A'` if `htype = true`, or +define for `adj = true` the continuous T-Lyapunov operator `L:X -> A*transpose(X)+X*transpose(A)` if `htype = false`, +or the continuous H-Lyapunov operator `L:X -> A*X'+isig*X*A'` if `htype = true`. """ -tlyapop(A::AbstractMatrix; disc=false, adj=false) = TLyapMap(A, ifelse(disc, Discrete(), Continuous()), adj) - -function LinearMaps._unsafe_mul!(y::AbstractVector, L::TLyapMap{T,TA,Continuous}, x::AbstractVector) where {T,TA} +function lyaplikeop(A::AbstractMatrix; disc=false, adj=false, isig = 1, htype = false) + LyapLikeMap(A, ifelse(disc, Discrete(), Continuous()), ifelse(htype, Htype(), Ttype()); adj, isig) +end +function LinearMaps._unsafe_mul!(y::AbstractVector, L::LyapLikeMap{T,TA,Continuous,Ttype}, x::AbstractVector) where {T,TA} m, n = size(L.A) T1 = promote_type(T, eltype(x)) X = L.adj ? reshape(convert(AbstractVector{T1}, x), m, n) : reshape(convert(AbstractVector{T1}, x), n, m) Y = L.adj ? L.A*transpose(X) : L.A*X - y[:] .= (Y+transpose(Y))[:] + y[:] .= L.isig == 1 ? (Y+transpose(Y))[:] : (Y-transpose(Y))[:] return y end -function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.TransposeMap{T,<:TLyapMap{T}}, y::AbstractVector) where {T} +function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.TransposeMap{T,<:LyapLikeMap{T,<:Any,Continuous,Ttype}}, y::AbstractVector) where {T} m = size(LT.lmap.A,1) T1 = promote_type(T, eltype(y)) Y = reshape(convert(AbstractVector{T1}, y), m, m) - x[:] = LT.lmap.adj ? ((Y+transpose(Y))*LT.lmap.A)[:] : (transpose(LT.lmap.A)*(Y+transpose(Y)))[:] + temp = LT.lmap.isig ==1 ? Y+transpose(Y) : Y-transpose(Y) + x[:] = LT.lmap.adj ? (temp*LT.lmap.A)[:] : (transpose(LT.lmap.A)*temp)[:] return x end -function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.AdjointMap{T,<:TLyapMap{T}}, y::AbstractVector) where {T} +function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.AdjointMap{T,<:LyapLikeMap{T,<:Any,Continuous,Ttype}}, y::AbstractVector) where {T} m = size(LT.lmap.A,1) T1 = promote_type(T, eltype(y)) Y = reshape(convert(AbstractVector{T1}, y), m, m) - x[:] = LT.lmap.adj ? ((Y+transpose(Y))*conj(LT.lmap.A))[:] : (LT.lmap.A'*(Y+transpose(Y)))[:] + temp = LT.lmap.isig ==1 ? Y+transpose(Y) : Y-transpose(Y) + x[:] = LT.lmap.adj ? (temp*conj(LT.lmap.A))[:] : (LT.lmap.A'*temp)[:] return x end -# -struct HLyapMap{T,TA <: AbstractMatrix,CD <: ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T} - A::TA - adj::Bool - function HLyapMap{T,TA,CD}(A::TA, adj=false) where {T,TA<:AbstractMatrix{T},CD} - return new{T,TA,CD}(A, adj) - end -end -HLyapMap(A::TA, ::CD = Continuous(), adj::Bool = false) where {T,TA<:AbstractMatrix{T},CD<:ContinuousOrDiscrete} = - HLyapMap{T,TA,CD}(A, adj) -#LinearAlgebra.transpose(L::HLyapunovMap{<:Any,<:Any,CD}) where {CD} = HLyapunovMap(L.U', CD(), !L.adj) - -function Base.size(L::HLyapMap) - m, n = size(L.A) - return (m*m, m*n) -end -""" - L = hlyapop(A; adj = false) - -Define, for an `m x n` matrix `A`, the continuous T-Lyapunov operator `L:X -> A*X+X'*A'` -if `adj = false` and `L:X -> A*X'+X*A'` if `adj = true`. -""" -function hlyapop(A::AbstractMatrix; disc=false, adj=false) - HLyapMap(A, ifelse(disc, Discrete(), Continuous()), adj) -end -function LinearMaps._unsafe_mul!(y::AbstractVector, L::HLyapMap{T,TA,Continuous}, x::AbstractVector) where {T,TA} +function LinearMaps._unsafe_mul!(y::AbstractVector, L::LyapLikeMap{T,TA,Continuous,Htype}, x::AbstractVector) where {T,TA} m, n = size(L.A) T1 = promote_type(T, eltype(x)) X = L.adj ? reshape(convert(AbstractVector{T1}, x), m, n) : reshape(convert(AbstractVector{T1}, x), n, m) Y = L.adj ? L.A*X' : L.A*X - y[:] .= (Y+Y')[:] + y[:] .= L.isig == 1 ? (Y+Y')[:] : (Y-Y')[:] return y end -function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.TransposeMap{T,<:HLyapMap{T}}, y::AbstractVector) where {T} +function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.TransposeMap{T,<:LyapLikeMap{T,<:Any,Continuous,Htype}}, y::AbstractVector) where {T} m = size(LT.lmap.A,1) T1 = promote_type(T, eltype(y)) Y = reshape(convert(AbstractVector{T1}, y), m, m) - x[:] = LT.lmap.adj ? ((Y+transpose(Y))*LT.lmap.A)[:] : (transpose(LT.lmap.A)*(Y+transpose(Y)))[:] + temp = LT.lmap.isig ==1 ? Y+transpose(Y) : Y-transpose(Y) + x[:] = LT.lmap.adj ? (temp*LT.lmap.A)[:] : (transpose(LT.lmap.A)*temp)[:] return x end -function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.AdjointMap{T,<:HLyapMap{T}}, y::AbstractVector) where {T} +function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.AdjointMap{T,<:LyapLikeMap{T,<:Any,Continuous,Htype}}, y::AbstractVector) where {T} m = size(LT.lmap.A,1) T1 = promote_type(T, eltype(y)) Y = reshape(convert(AbstractVector{T1}, y), m, m) - x[:] = LT.lmap.adj ? ((Y+Y')*LT.lmap.A)[:] : (LT.lmap.A'*(Y+Y'))[:] + temp = LT.lmap.isig ==1 ? Y+Y' : Y-Y' + x[:] = LT.lmap.adj ? (temp*LT.lmap.A)[:] : (LT.lmap.A'*temp)[:] return x end -struct UTTLyapunovMap{T,TU <: AbstractMatrix,CD <: ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T} +struct UTLyapLikeMap{T,TU <: AbstractMatrix,CD <: ContinuousOrDiscrete,TH<:TtypeOrHtype} <: LyapunovMatrixEquationsMaps{T} U::TU adj::Bool - function UTTLyapunovMap{T,TU,CD}(U::TU, adj=false) where {T,TU<:AbstractMatrix{T},CD} + function UTLyapLikeMap{T,TU,CD,TH}(U::TU, adj=false) where {T,TU<:AbstractMatrix{T},CD,TH} LinearAlgebra.checksquare(U) - (!adj && istriu(U.parent)) || (adj && istriu(U)) || error("U must be upper triangular") - return new{T,TU,CD}(U, adj) + istriu(U) || error("U must be upper triangular") + #(!adj && istriu(U.parent)) || (adj && istriu(U)) || error("U must be upper triangular") + return new{T,TU,CD,TH}(U, adj) end end -UTTLyapunovMap(U::TU, ::CD = Continuous(), adj::Bool = false) where {T,TU<:AbstractMatrix{T},CD<:ContinuousOrDiscrete} = - UTTLyapunovMap{T,TU,CD}(U, adj) +UTLyapLikeMap(U::TU, ::CD = Continuous(), ::TH = Ttype(); adj::Bool = false) where {T,TU<:AbstractMatrix{T},CD<:ContinuousOrDiscrete,TH<:TtypeOrHtype} = + UTLyapLikeMap{T,TU,CD,TH}(U, adj) -function Base.size(L::UTTLyapunovMap) +function Base.size(L::UTLyapLikeMap) n = size(L.U, 1) N = Int(n * (n + 1) / 2) return (N, N) end """ - L = tulyapop(U) - -Define, for an upper triangular matrix `U`, the continuous T-Lyapunov operator `L:X -> U*X+transpose(X)*transpose(U)`. + L = tulyaplikeop(U; adj = false) - L = tulyapop(transpose(U)) +Define, for an upper triangular matrix `U`, the continuous T-Lyapunov operator `L:X -> transpose(U)*X+transpose(X)*U`, if `adj = false`, +or `L:X -> U*transpose(X)+X*transpose(U)` if `adj = true`. +""" +function tulyaplikeop(U::AbstractMatrix; disc=false, adj = false) + UTLyapLikeMap(U, ifelse(disc, Discrete(), Continuous()), Ttype(); adj ) +end +""" + L = hulyaplikeop(U; adj = false) -Define, for an upper triangular matrix `U`, the continuous T-Lyapunov operator `L:X -> U*transpose(X)+X*transpose(U)`. +Define, for an upper triangular matrix `U`, the continuous H-Lyapunov operator `L:X -> U'*X+X'*U`, if `adj = false`, +`L:X -> U*X'+X*U'` if `adj = true`. """ -function tulyapop(U::AbstractMatrix; disc=false) - UTTLyapunovMap(U, ifelse(disc, Discrete(), Continuous()), !(typeof(U) <: Transpose)) +function hulyaplikeop(U::AbstractMatrix; disc=false, adj = false) + UTLyapLikeMap(U, ifelse(disc, Discrete(), Continuous()), Htype(); adj) end -function LinearMaps._unsafe_mul!(y::AbstractVector, L::UTTLyapunovMap{T,TU,Continuous}, x::AbstractVector) where {T,TU} + +function LinearMaps._unsafe_mul!(y::AbstractVector, L::UTLyapLikeMap{T,TU,Continuous,Ttype}, x::AbstractVector) where {T,TU} T1 = promote_type(T, eltype(x)) X = vec2triu(convert(AbstractVector{T1}, x), her=false) - Y = L.adj ? L.U*transpose(X) : L.U*X + Y = L.adj ? L.U*transpose(X) : transpose(L.U)*X y[:] .= triu2vec(Y+transpose(Y)) return y end -function LinearMaps._unsafe_mul!(y::AbstractVector, LT::LinearMaps.TransposeMap{T,<:UTTLyapunovMap{T}}, x::AbstractVector) where {T} +#LinearAlgebra.adjoint(L::UTLyapLikeMap{T,TU,Continuous,Ttype}) where {T,TU}= LinearAlgebra.transpose(L) +#LinearAlgebra.transpose(L::UTLyapLikeMap{T,TU,Continuous,Ttype}) where {T,TU}= LinearAlgebra.adjoint(L) + +function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.TransposeMap{T,<:UTLyapLikeMap{T,<:Any,Continuous,Ttype}}, y::AbstractVector) where {T} n = size(LT.lmap.U,2) - T1 = promote_type(T, eltype(x)) - X = vec2triu(convert(AbstractVector{T1}, x), her=false) - y[:] = LT.lmap.adj ? triu2vec((X+transpose(X))*LT.lmap.U) : triu2vec(transpose(LT.lmap.U)*(X+transpose(X))) - return y + T1 = promote_type(T, eltype(y)) + Y = vec2triu(convert(AbstractVector{T1}, y), her=false) + x[:] = LT.lmap.adj ? triu2vec((Y+transpose(Y))*LT.lmap.U) : triu2vec(LT.lmap.U*(Y+transpose(Y))) + #x[:] = LT.lmap.adj ? triu2vec(Y*LT.lmap.U+Y'*conj(LT.lmap.U)) : triu2vec(LT.lmap.U*Y+conj(LT.lmap.U)*Y') + return x end -function LinearMaps._unsafe_mul!(y::AbstractVector, LT::LinearMaps.AdjointMap{T,<:UTTLyapunovMap{T}}, x::AbstractVector) where {T} +function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.AdjointMap{T,<:UTLyapLikeMap{T,<:Any,Continuous,Ttype}}, y::AbstractVector) where {T} n = size(LT.lmap.U,2) + T1 = promote_type(T, eltype(y)) + Y = vec2triu(convert(AbstractVector{T1}, y), her=false) + #x[:] = LT.lmap.adj ? triu2vec((Y+Y')*conj(LT.lmap.U)) : triu2vec(conj(LT.lmap.U)*(Y+Y')) + x[:] = LT.lmap.adj ? triu2vec((Y+transpose(Y))*conj(LT.lmap.U)) : triu2vec(conj(LT.lmap.U)*(Y+transpose(Y))) + return x +end + +function LinearMaps._unsafe_mul!(y::AbstractVector, L::UTLyapLikeMap{T,TU,Continuous,Htype}, x::AbstractVector) where {T,TU} T1 = promote_type(T, eltype(x)) X = vec2triu(convert(AbstractVector{T1}, x), her=false) - y[:] = LT.lmap.adj ? triu2vec((X+transpose(X))*+conj(LT.lmap.U)) : triu2vec(LT.lmap.U'*(X+transpose(X))) + Y = L.adj ? L.U*X' : L.U'*X + y[:] .= triu2vec(Y+Y') return y end -struct UTHLyapunovMap{T,TU <: AbstractMatrix,CD <: ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T} - U::TU - adj::Bool - function UTHLyapunovMap{T,TU,CD}(U::TU, adj=false) where {T,TU<:AbstractMatrix{T},CD} - LinearAlgebra.checksquare(U) - (!adj && istriu(U.parent)) || (adj && istriu(U)) || error("U must be upper triangular") - return new{T,TU,CD}(U, adj) - end +# function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.TransposeMap{T,<:UTLyapLikeMap{T,TU,Continuous,Htype}}, y::AbstractVector) where {T,TU} +# n = size(LT.lmap.U,2) +# T1 = promote_type(T, eltype(y)) +# Y = vec2triu(convert(AbstractVector{T1}, y), her=false) +# x[:] = LT.lmap.adj ? triu2vec(Y*LT.lmap.U+Y'*conj(LT.lmap.U)) : triu2vec(LT.lmap.U*Y+conj(LT.lmap.U)*Y') +# return x +# end + +function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.AdjointMap{T,<:UTLyapLikeMap{T,TU,Continuous,Htype}}, y::AbstractVector) where {T,TU} + n = size(LT.lmap.U,2) + T1 = promote_type(T, eltype(y)) + Y = vec2triu(convert(AbstractVector{T1}, y), her=false) + # the following ensures Matrix(L') == Matrix(L)' + #x[:] = LT.lmap.adj ? triu2vec(Y*LT.lmap.U+Y'*conj(LT.lmap.U)) : triu2vec(LT.lmap.U*Y+conj(LT.lmap.U)*Y') + x[:] = LT.lmap.adj ? triu2vec((Y+Y')*LT.lmap.U) : triu2vec(LT.lmap.U*(Y+Y')) + return x end +# function LinearMaps._unsafe_mul!(y::AbstractVector, LT::LinearMaps.AdjointMap{T,<:UTHLyapunovMap{T}}, x::AbstractVector) where {T} +# n = size(LT.lmap.U,2) +# T1 = promote_type(T, eltype(x)) +# X = vec2triu(convert(AbstractVector{T1}, x), her=false) +# y[:] = LT.lmap.adj ? triu2vec((X+X')*LT.lmap.U) : triu2vec(LT.lmap.U'*(X+X')) +# return y +# end -UTHLyapunovMap(U::TU, ::CD = Continuous(), adj::Bool = false) where {T,TU<:AbstractMatrix{T},CD<:ContinuousOrDiscrete} = - UTHLyapunovMap{T,TU,CD}(U, adj) -#LinearAlgebra.transpose(L::UTHLyapunovMap{<:Any,<:Any,CD}) where {CD} = LinearAlgebra.adjoint(L) -function Base.size(L::UTHLyapunovMap) - n = size(L.U, 1) - N = Int(n * (n + 1) / 2) - return (N, N) -end +# struct UTLyapLikeMap{T,TU <: AbstractMatrix,CD <: ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T} +# U::TU +# adj::Bool +# function UTLyapLikeMap{T,TU,CD}(U::TU, adj=false) where {T,TU<:AbstractMatrix{T},CD} +# LinearAlgebra.checksquare(U) +# (!adj && istriu(U.parent)) || (adj && istriu(U)) || error("U must be upper triangular") +# return new{T,TU,CD}(U, adj) +# end +# end -""" - L = hulyapop(U) +# UTLyapLikeMap(U::TU, ::CD = Continuous(), adj::Bool = false) where {T,TU<:AbstractMatrix{T},CD<:ContinuousOrDiscrete} = +# UTLyapLikeMap{T,TU,CD}(U, adj) -Define, for an upper triangular matrix `U`, the continuous T-Lyapunov operator `L:X -> U*X+X'*U'`. +# #LinearAlgebra.transpose(L::UTLyapLikeMap{<:Any,<:Any,CD}) where {CD} = LinearAlgebra.adjoint(L) - L = hulyapop(U') +# function Base.size(L::UTLyapLikeMap) +# n = size(L.U, 1) +# N = Int(n * (n + 1) / 2) +# return (N, N) +# end -Define, for an upper triangular matrix `U`, the continuous T-Lyapunov operator `L:X -> U*X'+X*U'`. -""" -function hulyapop(U::AbstractMatrix; disc=false) - UTHLyapunovMap(U, ifelse(disc, Discrete(), Continuous()), !(typeof(U) <: Adjoint)) -end -function LinearMaps._unsafe_mul!(y::AbstractVector, L::UTHLyapunovMap{T,TU,Continuous}, x::AbstractVector) where {T,TU} - T1 = promote_type(T, eltype(x)) - X = vec2triu(convert(AbstractVector{T1}, x), her=false) - Y = L.adj ? L.U*X' : L.U*X - y[:] .= triu2vec(Y+Y') - return y -end +# function LinearMaps._unsafe_mul!(y::AbstractVector, L::UTLyapLikeMap{T,<:Any,Continuous,Htype}, x::AbstractVector) where {T,TU} +# T1 = promote_type(T, eltype(x)) +# X = vec2triu(convert(AbstractVector{T1}, x), her=false) +# Y = L.adj ? L.U*X' : L.U*X +# y[:] .= triu2vec(Y+Y') +# return y +# end -function LinearMaps._unsafe_mul!(y::AbstractVector, LT::LinearMaps.AdjointMap{T,<:UTHLyapunovMap{T}}, x::AbstractVector) where {T} - n = size(LT.lmap.U,2) - T1 = promote_type(T, eltype(x)) - X = vec2triu(convert(AbstractVector{T1}, x), her=false) - y[:] = LT.lmap.adj ? triu2vec((X+X')*LT.lmap.U) : triu2vec(LT.lmap.U'*(X+X')) - return y -end +# function LinearMaps._unsafe_mul!(y::AbstractVector, LT::LinearMaps.AdjointMap{T,<:UTLyapLikeMap{T}}, x::AbstractVector) where {T} +# n = size(LT.lmap.U,2) +# T1 = promote_type(T, eltype(x)) +# X = vec2triu(convert(AbstractVector{T1}, x), her=false) +# y[:] = LT.lmap.adj ? triu2vec((X+X')*LT.lmap.U) : triu2vec(LT.lmap.U'*(X+X')) +# return y +# end struct LyapunovMap{T,TA <: AbstractMatrix,CD <: ContinuousOrDiscrete} <: LyapunovMatrixEquationsMaps{T} A::TA @@ -1360,11 +1343,6 @@ struct GeneralizedSylvesterMap{T,TA <: AbstractMatrix,TB <: AbstractMatrix,TC <: end end -# LinearAlgebra.adjoint(L::GeneralizedSylvesterMap{T}) where T = -# GeneralizedSylvesterMap(L.A', L.B', L.C', L.D') -# LinearAlgebra.transpose(L::GeneralizedSylvesterMap{T}) where T = -# GeneralizedSylvesterMap(L.A', L.B', L.C', L.D') - """ M = sylvop(A, B, C, D) diff --git a/test/test_clyap.jl b/test/test_clyap.jl index aae8543..3cc7256 100644 --- a/test/test_clyap.jl +++ b/test/test_clyap.jl @@ -390,36 +390,36 @@ end @testset "Continuous positive Lyapunov-like equations" begin n = 5 - + Ty = Float64 # nonsingular U for Ty in (Float64, Float32, BigFloat, Double64) Ty == Float64 ? reltol = eps(float(100*n)) : reltol = eps(100*n*one(Ty)) U = triu(rand(Ty,n,n)); X0 = triu(rand(Ty,n,n)); Q = Matrix(Symmetric(transpose(U)*X0 + transpose(X0)*U)) - @time X = tulyapc!(U, copy(Q); adj = true) + @time X = tulyapc!(U, copy(Q); adj = false) @test norm(transpose(U)*X + transpose(X)*U - Q)/norm(X) < reltol Q = Matrix(Symmetric(U*transpose(X0) + X0*transpose(U))) - @time X = tulyapc!(U, copy(Q); adj = false); + @time X = tulyapc!(U, copy(Q); adj = true); @test norm(U*transpose(X) + X*transpose(U)- Q)/norm(X) < reltol U = triu(rand(Ty,n,n)+im*rand(Ty,n,n)); X0 = triu(rand(Ty,n,n)+im*rand(Ty,n,n)); Q = Matrix(Symmetric(transpose(U)*X0 + transpose(X0)*U)) - @time X = tulyapc!(U, copy(Q); adj = true) + @time X = tulyapc!(U, copy(Q); adj = false) @test norm(transpose(U)*X + transpose(X)*U - Q)/norm(X) < reltol Q = Matrix(Symmetric(U*transpose(X0) + X0*transpose(U))) - @time X = tulyapc!(U, copy(Q); adj = false); + @time X = tulyapc!(U, copy(Q); adj = true); @test norm(U*transpose(X) + X*transpose(U)- Q)/norm(X) < reltol Q = Matrix(Hermitian(U'*X0 + X0'*U)) - @time X = hulyapc!(U, copy(Q); adj = true) + @time X = hulyapc!(U, copy(Q); adj = false) @test norm(U'*X + X'*U - Q)/norm(X) < reltol Q = Matrix(Hermitian(U*X0' + X0*U')) - @time X = hulyapc!(U, copy(Q); adj = false); + @time X = hulyapc!(U, copy(Q); adj = true); @test norm(U*X' + X*U'- Q)/norm(X) < reltol end @@ -435,7 +435,7 @@ end @time X = tulyapc!(U, copy(Q); adj = false) @test norm(transpose(U)*X + transpose(X)*U - Q)/norm(X) < reltol - LT = MatrixEquations.tulyapop(transpose(U)) + LT = tulyaplikeop(U; adj = false) x2,info=MatrixEquations.cgls(LT,triu2vec(Q),reltol=1.e-14,maxiter=1000); X2 = vec2triu(x2); @test norm(transpose(U)*X2 + transpose(X2)*U - Q)/norm(X2) < reltol @@ -443,7 +443,7 @@ end @time X = tulyapc!(U, copy(Q); adj = true) @test norm(U*transpose(X) + X*transpose(U)- Q)/norm(X) < reltol - LT = MatrixEquations.tulyapop(U) + LT = tulyaplikeop(U; adj = true) x2,info=MatrixEquations.cgls(LT,triu2vec(Q),reltol=1.e-14,maxiter=1000); X2 = vec2triu(x2); @test norm(U*transpose(X2) + X2*transpose(U)- Q)/norm(X2) < reltol @@ -454,7 +454,7 @@ end @time X = tulyapc!(U, copy(Q); adj = false) @test norm(transpose(U)*X + transpose(X)*U - Q)/norm(X) < reltol - LT = MatrixEquations.tulyapop(transpose(U)) + LT = tulyaplikeop(U; adj = false) x2,info=MatrixEquations.cgls(LT,triu2vec(Q),reltol=1.e-14,maxiter=1000); X2 = vec2triu(x2); @test norm(transpose(U)*X2 + transpose(X2)*U - Q)/norm(X2) < reltol @@ -463,7 +463,7 @@ end @time X = tulyapc!(U, copy(Q); adj = true); @test norm(U*transpose(X) + X*transpose(U)- Q)/norm(X) < reltol - LT = MatrixEquations.tulyapop(U) + LT = tulyaplikeop(U; adj = true) x2,info=MatrixEquations.cgls(LT,triu2vec(Q),reltol=1.e-14,maxiter=1000); X2 = vec2triu(x2); @test norm(U*transpose(X2) + X2*transpose(U)- Q)/norm(X) < reltol @@ -474,7 +474,7 @@ end @time X = hulyapc!(U, copy(Q); adj = false) @test norm(U'*X + X'*U - Q)/norm(X) < reltol - LT = MatrixEquations.hulyapop(U') + LT = hulyaplikeop(U; adj = false) x2,info=MatrixEquations.cgls(LT,triu2vec(Q),reltol=1.e-14,maxiter=1000); X2 = vec2triu(x2); @test norm(U'*X2 + X2'*U - Q)/norm(X2) < reltol @@ -482,7 +482,7 @@ end @time X = hulyapc!(U, copy(Q); adj = true); @test norm(U*X' + X*U'- Q)/norm(X) < reltol - LT = MatrixEquations.hulyapop(U) + LT = hulyaplikeop(U; adj = true) x2,info=MatrixEquations.cgls(LT,triu2vec(Q),reltol=1.e-14,maxiter=1000); X2 = vec2triu(x2); @test norm(U*X2' + X2*U'- Q)/norm(X) < reltol diff --git a/test/test_iterative.jl b/test/test_iterative.jl index 27a154a..d936dfc 100644 --- a/test/test_iterative.jl +++ b/test/test_iterative.jl @@ -44,6 +44,8 @@ using IterativeSolvers X, info = tulyapci(U, Q, adj = false, reltol=1.e-14, maxiter=1000); @test norm(transpose(U)*X + transpose(X)*U - Q)/norm(X) < reltol + + Q = Matrix(Symmetric(U*transpose(X0) + X0*transpose(U))) @@ -132,7 +134,7 @@ using IterativeSolvers @testset "Generalized T/H-Sylvester" begin a = [rand(2,3)]; b = [I]; c = [I]; d = [a[1]']; L = gsylvop(a,b,c,d;nx=2) - L1 = tlyapop(a[1]) + L1 = lyaplikeop(a[1]) @test Matrix(L) == Matrix(L1) @test Matrix(L') == Matrix(L1') diff --git a/test/test_mecondest.jl b/test/test_mecondest.jl index b687962..f0b9f3f 100644 --- a/test/test_mecondest.jl +++ b/test/test_mecondest.jl @@ -5,23 +5,13 @@ using LinearMaps using MatrixEquations using Test -# function check_ctranspose(op::LinearMaps.LinearMap{T}) where {T <: Union{AbstractFloat, Complex}} -# (m, n) = size(op) -# x = rand(T,n) -# y = rand(T,m) -# yAx = dot(y, op * x) -# xAty = dot(x, op' * y) -# ε = eps(real(eltype(op))) -# return abs(yAx - conj(xAty)) < (abs(yAx) + ε) * ε^(1 / 3) -# end -function check_ctranspose(op::LinearMaps.LinearMap) +function check_ctranspose(op::LinearMaps.LinearMap{T}) where {T <: Union{AbstractFloat, Complex}} (m, n) = size(op) - T1 = promote_type(Float64,eltype(op)) - x = rand(T1,n) - y = rand(T1,m) + x = rand(T,n) + y = rand(T,m) yAx = dot(y, op * x) xAty = dot(x, op' * y) - ε = eps(real(T1)) + ε = eps(real(eltype(op))) return abs(yAx - conj(xAty)) < (abs(yAx) + ε) * ε^(1 / 3) end @@ -58,14 +48,14 @@ Tcrssym = lyapop(as,her=true); Tcrssyminv = invlyapop(as,her=true); # define some T/H-Lyapunov operators for upper triangular arguments -Ttur = tulyapop(ur) -Ttaur = tulyapop(transpose(ur)) -Ttuc = tulyapop(uc) -Ttauc = tulyapop(transpose(uc)) -Thur = hulyapop(ur) -Thaur = hulyapop(ur') -Thuc = hulyapop(uc) -Thauc = hulyapop(uc') +Ttur = tulyaplikeop(ur,adj = true) +Ttaur = tulyaplikeop(ur,adj = false) +Ttuc = tulyaplikeop(uc,adj = true) +Ttauc = tulyaplikeop(uc,adj = false) +Thur = hulyaplikeop(ur,adj = true) +Thaur = hulyaplikeop(ur,adj = false) +Thuc = hulyaplikeop(uc,adj = true) +Thauc = hulyaplikeop(uc,adj = false) Tcc = lyapop(ac);