Skip to content

Commit

Permalink
Adding some more tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@andreasvarga
andreasvarga committed Oct 16, 2023
1 parent 25ed44d commit 3786a30
Show file tree
Hide file tree
Showing 2 changed files with 56 additions and 14 deletions.
8 changes: 4 additions & 4 deletions src/meoperators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -210,7 +210,7 @@ function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.TransposeMap{
m = size(LT.lmap.A,1)
T1 = promote_type(T, eltype(y))
Y = reshape(convert(AbstractVector{T1}, y), m, m)
temp = LT.lmap.isig ==1 ? Y+transpose(Y) : Y-transpose(Y)
temp = LT.lmap.isig ==1 ? Y+transpose(Y) : (LT.lmap.adj ? transpose(Y)-Y : Y-transpose(Y))
x[:] = LT.lmap.adj ? (temp*LT.lmap.A)[:] : (transpose(LT.lmap.A)*temp)[:]
return x
end
Expand All @@ -219,7 +219,7 @@ function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.AdjointMap{T,
m = size(LT.lmap.A,1)
T1 = promote_type(T, eltype(y))
Y = reshape(convert(AbstractVector{T1}, y), m, m)
temp = LT.lmap.isig ==1 ? Y+transpose(Y) : Y-transpose(Y)
temp = LT.lmap.isig ==1 ? Y+transpose(Y) : (LT.lmap.adj ? transpose(Y)-Y : Y-transpose(Y))
x[:] = LT.lmap.adj ? (temp*conj(LT.lmap.A))[:] : (LT.lmap.A'*temp)[:]
return x
end
Expand All @@ -237,7 +237,7 @@ function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.TransposeMap{
m = size(LT.lmap.A,1)
T1 = promote_type(T, eltype(y))
Y = reshape(convert(AbstractVector{T1}, y), m, m)
temp = LT.lmap.isig ==1 ? Y+transpose(Y) : Y-transpose(Y)
temp = LT.lmap.isig ==1 ? Y+transpose(Y) : (LT.lmap.adj ? transpose(Y)-Y : Y-transpose(Y))
x[:] = LT.lmap.adj ? (temp*LT.lmap.A)[:] : (transpose(LT.lmap.A)*temp)[:]
return x
end
Expand All @@ -246,7 +246,7 @@ function LinearMaps._unsafe_mul!(x::AbstractVector, LT::LinearMaps.AdjointMap{T,
m = size(LT.lmap.A,1)
T1 = promote_type(T, eltype(y))
Y = reshape(convert(AbstractVector{T1}, y), m, m)
temp = LT.lmap.isig ==1 ? Y+Y' : Y-Y'
temp = LT.lmap.isig ==1 ? Y+Y' : (LT.lmap.adj ? Y'-Y : Y-Y')
x[:] = LT.lmap.adj ? (temp*LT.lmap.A)[:] : (LT.lmap.A'*temp)[:]
return x
end
Expand Down
62 changes: 52 additions & 10 deletions test/test_iterative.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,13 +45,13 @@ 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)))

X, info = tulyapci(U,Q,adj = true, reltol=1.e-14,maxiter=1000);
@test norm(U*transpose(X) + X*transpose(U)- Q)/norm(X) < reltol


# complex case
U = triu(rand(Ty,n,n)+im*rand(Ty,n,n)); U[1,1] = 0; U[n,n] = 0;
X0 = triu(rand(Ty,n,n)+im*rand(Ty,n,n)); X0[1,1] = 0; X0[n,n] = 0
Expand Down Expand Up @@ -87,47 +87,89 @@ using IterativeSolvers
A = rand(Ty,m,n);
X0 = rand(Ty,n,m);
C = Matrix(Symmetric(A*X0 + transpose(X0)*transpose(A)))
X, info = tlyapci(A, C, adj = false, reltol=1.e-14, maxiter=1000);
@time X, info = tlyapci(A, C, adj = false, reltol=1.e-14, maxiter=1000);
@test norm(A*X + transpose(X)*transpose(A) - C)/norm(X) < reltol

Y = A*X0; C = Y - transpose(Y);
@time X, info = tlyapci(A, C, -1; adj = false, reltol=1.e-14, maxiter=1000);
@test norm(A*X - transpose(X)*transpose(A) - C)/norm(X) < reltol


X0 = rand(Ty,m,n);
C = Matrix(Symmetric(A*transpose(X0)+X0*transpose(A)))
X, info = tlyapci(A, C, adj = true, reltol=1.e-14, maxiter=1000);
@time X, info = tlyapci(A, C, adj = true, reltol=1.e-14, maxiter=1000);
@test norm(A*transpose(X)+X*transpose(A) - C)/norm(X) < reltol

Y = A*transpose(X0); C = Y - transpose(Y);
@time X, info = tlyapci(A, C, -1; adj = true, reltol=1.e-14, maxiter=1000);
@test norm(A*transpose(X)-X*transpose(A) - C)/norm(X) < reltol


# complex case
A = rand(Ty,m,n)+im*rand(Ty,m,n);
X0 =rand(Ty,n,m)+im*rand(Ty,n,m);
C = Matrix(Symmetric(A*X0 + transpose(X0)*transpose(A)))
X, info = tlyapci(A, C, adj = false, reltol=1.e-14, maxiter=1000);
@time X, info = tlyapci(A, C, adj = false, reltol=1.e-14, maxiter=1000);
@test norm(A*X + transpose(X)*transpose(A) - C)/norm(X) < reltol

Y = A*X0; C = Y - transpose(Y);
@time X, info = tlyapci(A, C, -1; adj = false, reltol=1.e-14, maxiter=1000);
@test norm(A*X - transpose(X)*transpose(A) - C)/norm(X) < reltol


X0 =rand(Ty,m,n)+im*rand(Ty,m,n);
C = Matrix(Symmetric(A*transpose(X0)+X0*transpose(A)))
X, info = tlyapci(A, C, adj = true, reltol=1.e-14, maxiter=1000);
@time X, info = tlyapci(A, C, adj = true, reltol=1.e-14, maxiter=1000);
@test norm(A*transpose(X)+X*transpose(A) - C)/norm(X) < reltol

Y = A*transpose(X0); C = Y - transpose(Y);
@time X, info = tlyapci(A, C, -1; adj = true, reltol=1.e-14, maxiter=1000);
@test norm(A*transpose(X)-X*transpose(A) - C)/norm(X) < reltol


# H-Lyapunov
# real case
A = rand(Ty,m,n);
X0 = rand(Ty,n,m);
C = Matrix(Hermitian(A*X0 + X0'*A'))
X, info = hlyapci(A, C, adj = false, reltol=1.e-14, maxiter=1000);
@time X, info = hlyapci(A, C, adj = false, reltol=1.e-14, maxiter=1000);
@test norm(A*X +X'*A' - C)/norm(X) < reltol

Y = A*X0; C = Y - Y';
@time X, info = hlyapci(A, C, -1; adj = false, reltol=1.e-14, maxiter=1000);
@test norm(A*X - X'*A' - C)/norm(X) < reltol

X0 = rand(Ty,m,n);
C = Matrix(Hermitian(A*X0'+X0*A'))
X, info = hlyapci(A, C, adj = true, reltol=1.e-14, maxiter=1000);
@time X, info = hlyapci(A, C, adj = true, reltol=1.e-14, maxiter=1000);
@test norm(A*X'+X*A' - C)/norm(X) < reltol

Y = A*X0'; C = Y - Y';
@time X, info = hlyapci(A, C, -1; adj = true, reltol=1.e-14, maxiter=1000);
@test norm(A*X'-X*A' - C)/norm(X) < reltol


# complex case
A = rand(Ty,m,n)+im*rand(Ty,m,n);
X0 =rand(Ty,n,m)+im*rand(Ty,n,m);
C = Matrix(Hermitian(A*X0 + X0'*A'))
X, info = hlyapci(A, C, adj = false, reltol=1.e-14, maxiter=1000);
@time X, info = hlyapci(A, C, adj = false, reltol=1.e-14, maxiter=1000);
@test norm(A*X +X'*A' - C)/norm(X) < reltol

Y = A*X0; C = Y - Y';
@time X, info = hlyapci(A, C, -1; adj = false, reltol=1.e-14, maxiter=1000);
@test norm(A*X - X'*A' - C)/norm(X) < reltol


X0 = rand(Ty,m,n)+im*rand(Ty,m,n);
C = Matrix(Hermitian(A*X0'+X0*A'))
X, info = hlyapci(A, C, adj = true, reltol=1.e-14, maxiter=1000);
@time X, info = hlyapci(A, C, adj = true, reltol=1.e-14, maxiter=1000);
@test norm(A*X'+X*A' - C)/norm(X) < reltol

Y = A*X0'; C = Y - Y';
@time X, info = hlyapci(A, C, -1; adj = true, reltol=1.e-14, maxiter=1000);
@test norm(A*X'-X*A' - C)/norm(X) < reltol

end
end

Expand Down

2 comments on commit 3786a30

@andreasvarga
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register()

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/93549

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v2.3.0 -m "<description of version>" 3786a302a3a06fffc486c55f5fb8ff2848032860
git push origin v2.3.0

Please sign in to comment.