From 63e4ca3bb5193c09d28dd2987ee8ae156c9752c9 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Fri, 28 Jan 2022 12:30:28 +0100
Subject: [PATCH] add test for unitary similarity trasnform

---
 src/matrix_comps.jl       |  6 +++---
 test/test_matrix_comps.jl | 10 ++++++++++
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/src/matrix_comps.jl b/src/matrix_comps.jl
index bc2c5c5fd..5a94598e4 100644
--- a/src/matrix_comps.jl
+++ b/src/matrix_comps.jl
@@ -586,11 +586,11 @@ D̃ = D
 If `unitary=true`, `T` is assumed unitary and the matrix adjoint is used instead of the inverse.
 """
 function similarity_transform(sys::ST, T; unitary=false) where ST <: AbstractStateSpace
-    Tf = factorize(T)
     if unitary
-        A = Tf'sys.A*T
-        B = Tf'sys.B
+        A = T'sys.A*T
+        B = T'sys.B
     else
+        Tf = factorize(T)
         A = Tf\sys.A*T
         B = Tf\sys.B
     end
diff --git a/test/test_matrix_comps.jl b/test/test_matrix_comps.jl
index 428561ab9..fc9592bea 100644
--- a/test/test_matrix_comps.jl
+++ b/test/test_matrix_comps.jl
@@ -19,6 +19,16 @@ Ab,Bb,Cb,T = ControlSystems.balance_statespace(A,B,C)
 @test Cb*T ≈ C
 
 @test sysb.A ≈ Ab
+@test similarity_transform(sysb, T) ≈ sys
+
+
+U = svd(randn(2,2)).U
+syst = similarity_transform(sys, U, unitary = true)
+Ab,Bb,Cb,Db = ssdata(syst)
+@test Ab ≈ U'A*U
+@test Bb ≈ U'B
+@test Cb ≈ C*U
+
 
 @test ControlSystems.balance_transform(A,B,C) ≈ ControlSystems.balance_transform(sys)