Samdp likes to fail with InversionError when theta approaches pole

The issue pops up in CI periodically, like here.

The following code snippet reproduces the issue for me on Linux with PyPI packages:

import numpy as np
from pymortests.algorithms.samdp import conv_diff_1d_fd
from pymor.algorithms.samdp import samdp
from pymor.operators.numpy import NumpyMatrixOperator

def test_samdp(n, m, k, wanted, which, rng):
    A = conv_diff_1d_fd(n, 1, 1)
    Eop = None

    B = rng.standard_normal((n, m))
    C = rng.standard_normal((k, n))

    Aop = NumpyMatrixOperator(A)
    Bva = Aop.source.from_numpy(B.T)
    Cva = Aop.source.from_numpy(C)

    return samdp(Aop, Eop, Bva, Cva, wanted, which=which)

rng = np.random.default_rng(595)
test_samdp(50, 1, 3, 15, 'NR', rng)

@lbalicki, any ideas?

Start by reproducing the failure with the supplied example and inspect pymor.algorithms.samdp.samdp, using pymortests.algorithms.samdp.conv_diff_1d_fd and NumpyMatrixOperator as shown. Trace why the NR case raises InversionError when theta approaches a pole, then add regression coverage demonstrating that the supplied case no longer fails.