Repository: scipy/scipy
File: scipy/linalg/src/_matfuncs_sqrtm.c (C)

Code before:
        // Solve the Sylvester equation for the top right block
        dtrsyl_("N", "N", &int1, &halfn, &otherhalfn, T, &intbign, &T[halfn*(intbign + 1)], &intbign, &T[halfn*intbign], &intbign, &scale, &info);

        if (scale != 1.0)
        {
            for (i = 0; i < otherhalfn; i++)
            {
                for (j = 0; j < halfn; j++)
                {
                    T[(halfn + i)*bign + j] *= scale;
                }
            }
        }
    }
    // Propagate the Sylvester solver info up the recursion stack
    if (i1 || i2 || info)
    {
        return 1;
    } else {
        return 0;
    }
}


int
sqrtm_recursion_c(SCIPY_C* T, npy_intp bign, npy_intp n)
{
    int i, i1 = 0, i2 = 0, info = 0, j, halfn, otherhalfn, int1 = 1, intbign = (int)bign;
    float scale = 0.0;
    if (n == 1)
    {
        // Scalar block
        T[0] = csqrtf(T[0]);
        return 0;
    } else if (n == 2) {
        // 2x2 upper triangular special case

        // Check if the 2x2 block is exactly 0.0 and if so do nothing.
        if ((cabsf(T[0]) == 0.0) && (cabsf(T[bign + 1]) == 0.0) && (cabsf(T[bign]) == 0.0)) { return 0; }

        T[0] = csqrtf(T[0]);
        // n + 1 is the next diagonal from T[0]
        T[bign + 1] = csqrtf(T[bign + 1]);
#if defined(_MSC_VER)
        // MSVC does not support complex division
        SCIPY_C tt = CPLX_C(crealf(T[0]) + crealf(T[bign + 1]), cimagf(T[0]) + cimagf(T[bign + 1]));
        float cc = cabsf(tt);
        T[bign] = _FCmulcr(_FCmulcc(T[bign], conjf(tt)), 1 / (cc*cc));
#else
        T[bign] = T[bign] / (T[0] + T[bign + 1]);
#endif

Review comment (suggestion):
FWIW, I can delete this line and the full testsuite still passes locally (random deletion analysis).

```diff
--- a/scipy/linalg/src/_matfuncs_sqrtm.c
+++ b/scipy/linalg/src/_matfuncs_sqrtm.c
@@ -878,7 +878,6 @@ sqrtm_recursion_d(double* T, npy_intp bign, npy_intp n)
             info = 1;
         } else if (sylvester21 && (T[2*(intbign + 1)] == 0.0) && ((T[0] == 0.0) || (T[intbign + 1] == 0.0))) {
             T[2*intbign] = INFINITY;
-            T[2*intbign + 1] = INFINITY;
```

Code after:
            return 0;
        }
    } else {
        // Recursion
        halfn = (int)(n / 2);
        // Don't split the 2x2 block, check if the entry under separation is zero
        if (T[(halfn-1)*bign + halfn] != 0.0) { halfn++; }
        otherhalfn = (int)(n - halfn);

        // Top left block
        i1 = sqrtm_recursion_d(T, bign, halfn);
        // Bottom right block
        i2 = sqrtm_recursion_d(&T[halfn*(bign + 1)], bign, otherhalfn);

        // When n == 3, the Sylvester equation is 1x2 or 2x1. If both
        // sides have a zero diagonal, dtrsyl_ produces garbage instead
        // of inf. Detect this and set INFINITY so it propagates.
        int sylvester12 = (n == 3) && (halfn == 1);
        int sylvester21 = (n == 3) && (halfn == 2);
        if (sylvester12  && (T[0] == 0.0) && ((T[intbign + 1] == 0.0 )|| (T[2*(intbign + 1)] == 0.0))) {
            T[intbign] = INFINITY;
            T[2*intbign] = INFINITY;
            info = 1;
        } else if (sylvester21 && (T[2*(intbign + 1)] == 0.0) && ((T[0] == 0.0) || (T[intbign + 1] == 0.0))) {
            T[2*intbign] = INFINITY;
            T[2*intbign + 1] = INFINITY;
            info = 1;
        } else {
            dtrsyl_("N", "N", &int1, &halfn, &otherhalfn, T, &intbign, &T[halfn*(intbign + 1)], &intbign, &T[halfn*intbign], &intbign, &scale, &info);
            if (scale != 1.0)
            {
                for (i = 0; i < otherhalfn; i++)
                {
                    for (j = 0; j < halfn; j++)
                    {
                        T[(halfn + i)*bign + j] *= scale;
                    }
                }
            }
        }
    }
    // Propagate the Sylvester solver info up the recursion stack
    if (i1 || i2 || info)
    {
        return 1;
    } else {
        return 0;
    }
}