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

Code before:
        }
    } 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_s(T, bign, halfn);
        // Bottom right block
        i2 = sqrtm_recursion_s(&T[halfn*(bign + 1)], bign, otherhalfn);

        // Solve the Sylvester equation for the top right block
        strsyl_("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_d(double* T, npy_intp bign, npy_intp n)
{
    int i, i1 = 0, i2 = 0, info = 0, j, halfn, otherhalfn, int1 = 1, intbign = (int)bign;
    double scale = 0.0, a, b, c, d, alpha, theta, mu;
    if (n == 1)
    {
        // Scalar block
        T[0] = sqrt(T[0]);
        return 0;
    } else if (n == 2) {
        a = T[0];
        c = T[1];
        b = T[bign];
        d = T[bign + 1];

Review comment (suggestion):
If I delete the `info = 1;` here and above, the full testsuite still passes locally (random deletion analysis).

Code after:
        }
    } 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_s(T, bign, halfn);
        // Bottom right block
        i2 = sqrtm_recursion_s(&T[halfn*(bign + 1)], bign, otherhalfn);

        // When n == 3, the Sylvester equation is 1x2 or 2x1. If both
        // sides have a zero diagonal, strsyl_ 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.0f) && ((T[intbign + 1] == 0.0f) || (T[2*(intbign + 1)] == 0.0f))) {
            T[intbign] = INFINITY;
            T[2*intbign] = INFINITY;
            info = 1;
        } else if (sylvester21 && (T[2*(intbign + 1)] == 0.0f) && ((T[0] == 0.0f) || (T[intbign + 1] == 0.0f))) {
            T[2*intbign] = INFINITY;
            T[2*intbign + 1] = INFINITY;
            info = 1;
        } else {
            strsyl_("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