diff --git a/matfree/eig.py b/matfree/eig.py
index c6bc9a9..0c40b37 100644
--- a/matfree/eig.py
+++ b/matfree/eig.py
@@ -28,7 +28,7 @@ def svd(Av: Callable, v0: Array, *parameters):
         U, S, Vt = linalg.svd(B, full_matrices=False)
 
         # Combine orthogonal transformations
-        return u @ U, S, Vt @ v.T
+        return (u @ U).T, S, Vt @ v.T
 
     return svd
 
diff --git a/tests/test_eig/test_svd_partial.py b/tests/test_eig/test_svd_partial.py
index 98c958b..ac9deec 100644
--- a/tests/test_eig/test_svd_partial.py
+++ b/tests/test_eig/test_svd_partial.py
@@ -4,62 +4,48 @@
 from matfree.backend import linalg, np, testing
 
 
-@testing.fixture()
 @testing.parametrize("nrows", [10])
 @testing.parametrize("ncols", [3])
-@testing.parametrize("num_significant_singular_vals", [3])
-def A(nrows, ncols, num_significant_singular_vals):
-    """Make a positive definite matrix with certain spectrum."""
-    # 'Invent' a spectrum. Use the number of pre-defined eigenvalues.
-    n = min(nrows, ncols)
-    d = np.arange(n) + 10.0
-    d = d.at[num_significant_singular_vals:].set(0.001)
-    return test_util.asymmetric_matrix_from_singular_values(d, nrows=nrows, ncols=ncols)
-
-
-def test_equal_to_linalg_svd(A):
+def test_equal_to_linalg_svd(nrows, ncols):
     """The output of full-depth SVD should be equal (*) to linalg.svd().
 
     (*) Note: The singular values should be identical,
     and the orthogonal matrices should be orthogonal. They are not unique.
     """
-    nrows, ncols = np.shape(A)
-    num_matvecs = min(nrows, ncols)
-
+    d = np.arange(1.0, 1.0 + min(nrows, ncols))
+    A = test_util.asymmetric_matrix_from_singular_values(d, nrows=nrows, ncols=ncols)
     v0 = np.ones((ncols,))
-    v0 /= linalg.vector_norm(v0)
+    num_matvecs = min(nrows, ncols)
 
     bidiag = decomp.bidiag(num_matvecs)
     svd = eig.svd_partial(bidiag)
-    U, S, Vt = svd(lambda v, p: p @ v, v0, A)
+    Ut, S, Vt = svd(lambda v, p: p @ v, v0, A)
 
     U_, S_, Vt_ = linalg.svd(A, full_matrices=False)
 
     tols_decomp = {"atol": 1e-5, "rtol": 1e-5}
-    assert np.allclose(S, S_)  # strict "allclose"
-
-    assert np.shape(U) == np.shape(U_)
-    assert np.shape(Vt) == np.shape(Vt_)
-    assert np.allclose(U @ U.T, U_ @ U_.T, **tols_decomp)
+    assert np.allclose(S, S_, **tols_decomp)
+    assert np.allclose(Ut.T @ Ut, U_ @ U_.T, **tols_decomp)
     assert np.allclose(Vt @ Vt.T, Vt_ @ Vt_.T, **tols_decomp)
 
 
-def test_shapes_as_expected(A):
+@testing.parametrize("nrows", [10])
+@testing.parametrize("ncols", [3])
+@testing.parametrize("num_matvecs", [0, 2, 3])
+def test_shapes_as_expected(nrows, ncols, num_matvecs):
     """The output of full-depth SVD should be equal (*) to linalg.svd().
 
     (*) Note: The singular values should be identical,
     and the orthogonal matrices should be orthogonal. They are not unique.
     """
-    nrows, ncols = np.shape(A)
-    num_matvecs = min(nrows, ncols) // 2
-
+    d = np.arange(1.0, 1.0 + min(nrows, ncols))
+    A = test_util.asymmetric_matrix_from_singular_values(d, nrows=nrows, ncols=ncols)
     v0 = np.ones((ncols,))
-    v0 /= linalg.vector_norm(v0)
 
     bidiag = decomp.bidiag(num_matvecs)
     svd = eig.svd_partial(bidiag)
-    U, S, Vt = svd(lambda v, p: p @ v, v0, A)
+    Ut, S, Vt = svd(lambda v, p: p @ v, v0, A)
 
-    assert U.shape == (nrows, num_matvecs)
+    assert Ut.shape == (num_matvecs, nrows)
     assert S.shape == (num_matvecs,)
     assert Vt.shape == (num_matvecs, ncols)