Repository: scipy/scipy
File: scipy/linalg/src/_linalg_inv.hh (C++)

Code before:
inline void invert_slice_diagonal(
    CBLAS_INT N, T *data, SliceStatus& status
) {
    using real_type = typename type_traits<T>::real_type;
    using value_type = typename type_traits<T>::value_type;
    value_type *pdata = reinterpret_cast<value_type *>(data);

    value_type zero(0.), one(1.);
    real_type maxa(0.), maxinva(0.);

    for (CBLAS_INT j=0; j<N; j++) {
        value_type ajj = pdata[j*N + j];

        status.is_singular = (ajj == zero);
        if (status.is_singular) {
            status.lapack_info  = j;
            return;
        }

        value_type inv_ajj = one / ajj;
        pdata[j*N + j] = inv_ajj;

        // condition number
        real_type absa = std::abs(ajj), absinva = std::abs(inv_ajj);

        if(absa != absa) {
            status.is_ill_conditioned = true;
        } else {
            if(absa > maxa) {maxa = absa;}
            if(absinva > maxinva) {maxinva = absinva;}
        }
    }
    status.is_ill_conditioned = status.is_ill_conditioned || (maxa * maxinva > 1./ numeric_limits<real_type>::eps);
    status.rcond = maxa * maxinva;
}


template<typename T>
int
_inverse(PyArrayObject* ap_Am, T* ret_data, St structure, int lower, int overwrite_a, SliceStatusVec& vec_status)
{
    using real_type = typename type_traits<T>::real_type; // float if T==npy_cfloat etc

    npy_intp lower_band = 0, upper_band = 0;
    bool is_symm_or_herm = false, is_symm_not_herm = false;
    char uplo = lower ? 'L' : 'U';
    St slice_structure = St::NONE;
    bool posdef_fallback = true;
    SliceStatus slice_status;

    // --------------------------------------------------------------------

Review comment (suggestion):
Any reason to prefer this over `std::isnan(absa)`?

And in my usual "is this really tested" stance, if I change this to a codepath that can't be hit `spin test -j 8` still passes in my hands on x86-64 Linux:

```diff
-        if(absa != absa) {
+        if(0) {
             status.is_ill_conditioned = true;
```

Code after:
inline void invert_slice_diagonal(
    CBLAS_INT N, T *data, SliceStatus& status
) {
    using real_type = typename type_traits<T>::real_type;
    using value_type = typename type_traits<T>::value_type;
    value_type *pdata = reinterpret_cast<value_type *>(data);

    value_type zero(0.), one(1.);
    real_type maxa(0.), maxinva(0.);

    for (CBLAS_INT j=0; j<N; j++) {
        value_type ajj = pdata[j*N + j];

        status.is_singular = (ajj == zero);
        if (status.is_singular) {
            status.lapack_info  = j;
            return;
        }

        value_type inv_ajj = one / ajj;
        pdata[j*N + j] = inv_ajj;

        // condition number
        real_type absa = std::abs(ajj), absinva = std::abs(inv_ajj);

        if(absa > maxa) {maxa = absa;}
        if(absinva > maxinva) {maxinva = absinva;}
    }
    status.is_ill_conditioned = maxa * maxinva > 1./ numeric_limits<real_type>::eps;
    status.rcond = maxa * maxinva;
}


template<typename T>
int
_inverse(PyArrayObject* ap_Am, T* ret_data, St structure, int lower, int overwrite_a, SliceStatusVec& vec_status)
{
    using real_type = typename type_traits<T>::real_type; // float if T==npy_cfloat etc

    npy_intp lower_band = 0, upper_band = 0;
    bool is_symm_or_herm = false, is_symm_not_herm = false;
    char uplo = lower ? 'L' : 'U';
    St slice_structure = St::NONE;
    bool posdef_fallback = true;
    SliceStatus slice_status;

    // --------------------------------------------------------------------
    // Input Array Attributes
    // --------------------------------------------------------------------
    T* Am_data = (T *)PyArray_DATA(ap_Am);
    int ndim = PyArray_NDIM(ap_Am);              // Number of dimensions