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

Code before:
        assert_array_equal(sqrtm(A), A)

    def test_negative_strides(self):
        rng = np.random.default_rng(1738151906092738)
        A = rng.uniform(size=(3, 5, 5))
        A_negneg_orig = A[:, ::-1, ::-1]
        A_negneg_copy = A[:, ::-1, ::-1].copy()
        assert_allclose(sqrtm(A_negneg_orig), sqrtm(A_negneg_copy))

        A_posneg_orig = A[:, :, ::-1]
        A_posneg_copy = A[:, :, ::-1].copy()
        assert_allclose(sqrtm(A_posneg_orig), sqrtm(A_posneg_copy))

        A_negpos_orig = A[:, ::-1, :]
        A_negpos_copy = A[:, ::-1, :].copy()
        assert_allclose(sqrtm(A_negpos_orig), sqrtm(A_negpos_copy))


class TestFractionalMatrixPower:
    def test_round_trip_random_complex(self):
        rng = np.random.default_rng(1234)
        for p in range(1, 5):
            for n in range(1, 5):
                M_unscaled = (rng.standard_normal((n, n)) +
                              1j * rng.standard_normal((n, n)))
                for scale in np.logspace(-4, 4, 9):
                    M = M_unscaled * scale
                    M_root = fractional_matrix_power(M, 1/p)
                    M_round_trip = np.linalg.matrix_power(M_root, p)
                    assert_allclose(M_round_trip, M)

    def test_round_trip_random_float(self):
        # This test is more annoying because it can hit the branch cut;
        # this happens when the matrix has an eigenvalue
        # with no imaginary component and with a real negative component,
        # and it means that the principal branch does not exist.
        rng = np.random.default_rng(1234)
        for p in range(1, 5):
            for n in range(1, 5):
                M_unscaled = rng.standard_normal((n, n))
                for scale in np.logspace(-4, 4, 9):
                    M = M_unscaled * scale
                    M_root = fractional_matrix_power(M, 1/p)
                    M_round_trip = np.linalg.matrix_power(M_root, p)
                    assert_allclose(M_round_trip, M)

    def test_larger_abs_fractional_matrix_powers(self):
        rng = np.random.default_rng(1234)
        for n in (2, 3, 5):
            for i in range(10):
                M = rng.standard_normal((n, n)) + 1j * rng.standard_normal((n, n))

Review comment (bug):
I note that while this line fails before and passes after the patch, the actual warning behavior that was the original source of the matching ticket is not a problem for me on `main` at the moment locally on ARM Mac (which Ilayn had also observed over in the issue I believe).

It looks like the issue is reported for x86_64 Linux + OpenBLAS. However, even on that platform the behavior is the same for me--I only see a failure at the `isinf` line, not with regard to the warning behavior, when the source patch is reverted.

@aravindh-krishnamoorthy did you want to take a look?

Code after:
        assert_array_equal(sqrtm(A), A)

    def test_negative_strides(self):
        rng = np.random.default_rng(1738151906092738)
        A = rng.uniform(size=(3, 5, 5))
        A_negneg_orig = A[:, ::-1, ::-1]
        A_negneg_copy = A[:, ::-1, ::-1].copy()
        assert_allclose(sqrtm(A_negneg_orig), sqrtm(A_negneg_copy))

        A_posneg_orig = A[:, :, ::-1]
        A_posneg_copy = A[:, :, ::-1].copy()
        assert_allclose(sqrtm(A_posneg_orig), sqrtm(A_posneg_copy))

        A_negpos_orig = A[:, ::-1, :]
        A_negpos_copy = A[:, ::-1, :].copy()
        assert_allclose(sqrtm(A_negpos_orig), sqrtm(A_negpos_copy))

    @pytest.mark.parametrize('dtyp', [np.float32, np.float64,
                                      np.complex64, np.complex128])
    def test_gh24508(self, dtyp):
        # Check that sqrtm raises a warning when the input is singular
        # And for this example, the top row is filled with infs.
        a = np.array([[0, 2, 2], [0, 0, 0], [0, 0, 0]], dtype=dtyp)
        with pytest.warns(LinAlgWarning, match="Matrix is singular."):
            res = sqrtm(a)
            assert np.isinf(res[0, 1:]).all()


class TestFractionalMatrixPower:
    def test_round_trip_random_complex(self):
        rng = np.random.default_rng(1234)
        for p in range(1, 5):
            for n in range(1, 5):
                M_unscaled = (rng.standard_normal((n, n)) +
                              1j * rng.standard_normal((n, n)))
                for scale in np.logspace(-4, 4, 9):
                    M = M_unscaled * scale
                    M_root = fractional_matrix_power(M, 1/p)
                    M_round_trip = np.linalg.matrix_power(M_root, p)
                    assert_allclose(M_round_trip, M)

    def test_round_trip_random_float(self):
        # This test is more annoying because it can hit the branch cut;
        # this happens when the matrix has an eigenvalue
        # with no imaginary component and with a real negative component,
        # and it means that the principal branch does not exist.
        rng = np.random.default_rng(1234)
        for p in range(1, 5):
            for n in range(1, 5):
                M_unscaled = rng.standard_normal((n, n))
                for scale in np.logspace(-4, 4, 9):