diff --git a/control/canonical.py b/control/canonical.py
index c84ac1f8f..9c9a2a738 100644
--- a/control/canonical.py
+++ b/control/canonical.py
@@ -450,9 +450,7 @@ def bdschur(a, condmax=None, sort=None):
         aschur, tschur, condmax)
 
     if sort in ('continuous', 'discrete'):
-
         idxs = np.cumsum(np.hstack([0, blksizes[:-1]]))
-
         ev_per_blk = [complex(eigvals[i].real, abs(eigvals[i].imag))
                       for i in idxs]
 
@@ -470,7 +468,7 @@ def bdschur(a, condmax=None, sort=None):
         permidx = np.hstack([blkidxs[i] for i in sortidx])
         rperm = np.eye(amodal.shape[0])[permidx]
 
-        tmodal = tmodal @ rperm
+        tmodal = tmodal @ rperm.T
         amodal = rperm @ amodal @ rperm.T
         blksizes = blksizes[sortidx]
 
diff --git a/control/tests/canonical_test.py b/control/tests/canonical_test.py
index f822955fc..ecdaa04cb 100644
--- a/control/tests/canonical_test.py
+++ b/control/tests/canonical_test.py
@@ -287,6 +287,8 @@ def test_bdschur_sort(eigvals, sorted_blk_eigvals, sort):
 
     b, t, blksizes = bdschur(a, sort=sort)
     assert len(blksizes) == len(sorted_blk_eigvals)
+    np.testing.assert_allclose(a, t @ b @ t.T)
+    np.testing.assert_allclose(t.T, np.linalg.inv(t))
 
     blocks = extract_bdiag(b, blksizes)
     for block, blk_eigval in zip(blocks, sorted_blk_eigvals):