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

Code before:
        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

        op = interface.LinearOperator(
            shape=args.shape, dtype=args.op_dtype,
            matvec=identity, rmatvec=identity
        )

        self.check_matvec(op, data_dtype=args.data_dtype, complex_data=args.complex)
        self.check_matmat(op, data_dtype=args.data_dtype, complex_data=args.complex)
    
    @pytest.mark.parametrize("args", all_args_list)
    def test_identity_nonsquare(self, args):
        """Identity operator with zero-padding on non-square matrices"""
        def mv(x):
            # handle column vectors too
            # (`LinearOperator` handles reshape in post-processing)
            x = x.flatten()

            match np.sign(x.shape[0] - args.shape[-2]):

Review comment (suggestion):
```suggestion
        rtol = 3e-12 if np.finfo(data_dtype).eps < 1e-8 else 6e-4
```

Code after:
        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):
        """Simple identity operator on square matrices"""
        def identity(x):
            return x

        op = interface.LinearOperator(
            shape=args.shape, dtype=args.op_dtype,
            matvec=identity, rmatvec=identity
        )

        self.check_matvec(op, data_dtype=args.data_dtype, complex_data=args.complex)
        self.check_matmat(op, data_dtype=args.data_dtype, complex_data=args.complex)
    
    @pytest.mark.parametrize("args", all_args_list)
    def test_identity_nonsquare(self, args):
        """Identity operator with zero-padding on non-square matrices"""
        def mv(x):
            # handle column vectors too
            # (`LinearOperator` handles reshape in post-processing)