Repository: scipy/scipy
File: scipy/sparse/linalg/tests/test_interface.py (Python)

Code before:

        Data is generated with the real dtype `data_dtype` and operated on by the
        linear operator `op`.

        If `complex_data` is set to `True`, complex data is instead generated
        by combining randomised real and imaginary components, each generated
        with `data_dtype`.

        If `check_operators` is set to `True`, equivalence is checked between
        `matmat` and `@`, and between `rmatvec` and the composition of `.H` and `@`.

        If `check_dot` is set to `True`, equivalence is checked between
        `matmat` and `.dot`,
        and between `rmatmat` and the composition of `.H` and `.dot`.
        """
        rng = np.random.default_rng()
        k = rng.integers(2, 100)

        U = rng.standard_normal(size=(op.shape[-1], k), dtype=data_dtype)
        if complex_data:
            U = U + (1j * rng.standard_normal(size=(op.shape[-1], k), dtype=data_dtype))

        V = rng.standard_normal(size=(op.shape[-2], k), dtype=data_dtype)
        if complex_data:
            V = V + (1j * rng.standard_normal(size=(op.shape[-2], k), dtype=data_dtype))

        op_U = op.matmat(U)
        opH_V = op.rmatmat(V)

        if check_operator:
            assert_allclose(op_U, op @ U)
            assert_allclose(opH_V, op.H @ V)

        if check_dot:
            assert_allclose(op_U, op.dot(U))
            assert_allclose(opH_V, op.H.dot(V))

        op_U_H = np.conj(op_U).T
        UH = np.conj(U).T

        op_U_H_V = np.matmul(op_U_H, V)
        UH_opH_V = np.matmul(UH, opH_V)

        rtol = 1e-12 if np.finfo(data_dtype).eps < 1e-8 else 5e-4
        assert_allclose(op_U_H_V, UH_opH_V, rtol=rtol)

    @pytest.mark.parametrize("args", square_args_list)
    def test_identity_square(self, args):
        """Simple identity operator on square matrices"""
        def identity(x):
            return x

Review comment (suggestion):
```suggestion
        U = rng.standard_normal(size=(op.shape[-1], k), dtype=data_dtype)
        V = rng.standard_normal(size=(op.shape[-2], k), dtype=data_dtype)

        if complex_data:
            U = U + (1j * rng.standard_normal(size=(op.shape[-1], k), dtype=data_dtype))
            V = V + (1j * rng.standard_normal(size=(op.shape[-2], k), dtype=data_dtype))
```

Code after:
        This check verifies the equivalence of the forward and adjoint computation,
        using `matmat` and `rmatmat` respectively, on randomised data.

        Data is generated with the real dtype `data_dtype` and operated on by the
        linear operator `op`.

        If `complex_data` is set to `True`, complex data is instead generated
        by combining randomised real and imaginary components, each generated
        with `data_dtype`.

        If `check_operators` is set to `True`, equivalence is checked between
        `matmat` and `*` and `@`,
        and between `rmatvec` and the composition of `.H` and `@`.

        If `check_dot` is set to `True`, equivalence is checked between
        `matmat` and `.dot`,
        and between `rmatmat` and the composition of `.H` and `.dot`.
        """
        rng = np.random.default_rng(42)
        k = rng.integers(2, 100)

        U = rng.standard_normal(size=(op.shape[-1], k), dtype=data_dtype)
        V = rng.standard_normal(size=(op.shape[-2], k), dtype=data_dtype)
        if complex_data:
            U = U + (1j * rng.standard_normal(size=(op.shape[-1], k), dtype=data_dtype))
            V = V + (1j * rng.standard_normal(size=(op.shape[-2], k), dtype=data_dtype))

        op_U = op.matmat(U)
        opH_V = op.rmatmat(V)

        if check_operators:
            assert_allclose(op_U, op * U)
            assert_allclose(op_U, op @ U)
            assert_allclose(opH_V, op.H * V)
            assert_allclose(opH_V, op.H @ V)

        if check_dot:
            assert_allclose(op_U, op.dot(U))
            assert_allclose(opH_V, op.H.dot(V))

        op_U_H = np.conj(op_U).T
        UH = np.conj(U).T

        op_U_H_V = np.matmul(op_U_H, V)
        UH_opH_V = np.matmul(UH, opH_V)

        rtol = 3e-12 if np.finfo(data_dtype).eps < 1e-8 else 6e-4
        assert_allclose(op_U_H_V, UH_opH_V, rtol=rtol)

    @pytest.mark.parametrize("args", square_args_list)
    def test_identity_square(self, args):