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
        `matvec` and `*` and `@`,
        and between `rmatvec` and the composition of `.H` and `*`.

        If `check_dot` is set to `True`, equivalence is checked between
        `matvec` and `.dot`,
        and between `rmatvec` and the composition of `.H` and `.dot`.
        """
        rng = np.random.default_rng(42)

        dtype = np.dtype(data_dtype)
        *batch_shape, M, N = op.shape

        u_broadcast, v_broadcast = generate_broadcastable_shapes(
            2, ndim=3, min=0, max=5, rng=63
        )

        for u_shape, v_shape in (
            # test broadcasting for unbatched RHS
            ((N,), (M,)),
            # as well as with the batch shape of `op` + 3-D broadcast dims
            ((*u_broadcast, *batch_shape, N), (*v_broadcast, *batch_shape, M))
        ):
            u = rng.standard_normal(u_shape, dtype=dtype)
            v = rng.standard_normal(v_shape, dtype=dtype)
            if complex_data:
                u = u + (1j * rng.standard_normal(u_shape, dtype=dtype))
                v = v + (1j * rng.standard_normal(v_shape, dtype=dtype))
            u = xp.asarray(u)
            v = xp.asarray(v)
    
            op_u = op.matvec(u)
            opH_v = op.rmatvec(v)
    
            op_u_H_v = xp.vecdot(op_u, v, axis=-1)
            uH_opH_v = xp.vecdot(u, opH_v, axis=-1)
    
            rtol = 1e-12 if np.finfo(data_dtype).eps < 1e-8 else 1e-5
            atol = 2e-15 if np.finfo(data_dtype).eps < 1e-8 else 1e-5
            xp_assert_close(op_u_H_v, uH_opH_v, rtol=rtol, atol=atol)

    def check_matmat(
        self, xp, op: interface.LinearOperator, data_dtype: str,

Review comment (bug):
Good. It's surprising this escaped attention for so long. I'm not sure if you observed this, but when I just plugged in `rng=rng`, I actually wasn't seeing the failures locally.

Code after:
        """
        This check verifies the equivalence of the forward and adjoint computation,
        using `matvec` and `rmatvec` 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
        `matvec` and `*` and `@`,
        and between `rmatvec` and the composition of `.H` and `*`.

        If `check_dot` is set to `True`, equivalence is checked between
        `matvec` and `.dot`,
        and between `rmatvec` and the composition of `.H` and `.dot`.
        """
        rng = np.random.default_rng(42)

        dtype = np.dtype(data_dtype)
        *batch_shape, M, N = op.shape

        u_broadcast, v_broadcast = generate_broadcastable_shapes(
            2, ndim=3, min=0, max=5, rng=63
        )

        for u_shape, v_shape in (
            # test broadcasting for unbatched RHS
            ((N,), (M,)),
            # as well as with the batch shape of `op` + 3-D broadcast dims
            ((*u_broadcast, *batch_shape, N), (*v_broadcast, *batch_shape, M))
        ):
            u = rng.standard_normal(u_shape, dtype=dtype)
            v = rng.standard_normal(v_shape, dtype=dtype)
            if complex_data:
                u = u + (1j * rng.standard_normal(u_shape, dtype=dtype))
                v = v + (1j * rng.standard_normal(v_shape, dtype=dtype))
            u = xp.asarray(u)
            v = xp.asarray(v)
    
            op_u = op.matvec(u)
            opH_v = op.rmatvec(v)
    
            op_u_H_v = xp.vecdot(op_u, v, axis=-1)
            uH_opH_v = xp.vecdot(u, opH_v, axis=-1)
    
            rtol = 1e-12 if np.finfo(data_dtype).eps < 1e-8 else 1e-5
            atol = 2e-15 if np.finfo(data_dtype).eps < 1e-8 else 1e-5
            xp_assert_close(op_u_H_v, uH_opH_v, rtol=rtol, atol=atol)