Skip to content

Commit

Permalink
Updating some operator definitions
Browse files Browse the repository at this point in the history
  • Loading branch information
@andreasvarga
andreasvarga committed Oct 16, 2023
1 parent 81e6f0b commit 25ed44d
Show file tree
Hide file tree
Showing 9 changed files with 260 additions and 270 deletions.
2 changes: 1 addition & 1 deletion docs/make.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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",
Expand Down
11 changes: 5 additions & 6 deletions docs/src/meoperators.md
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down
3 changes: 2 additions & 1 deletion src/MatrixEquations.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
66 changes: 43 additions & 23 deletions src/iterative_methods.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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`
Expand All @@ -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
Expand All @@ -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`
Expand All @@ -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
Expand Down
102 changes: 51 additions & 51 deletions src/lyapunov.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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.
Expand All @@ -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)
Expand Down
Loading

0 comments on commit 25ed44d

Please sign in to comment.