GeoStat-Framework / pentapy

A Python toolbox for pentadiagonal linear systems
https://pentapy.readthedocs.io/
MIT License
14 stars 4 forks source link

[Bug] Cython division optimizations are too aggressive for the error handling implemented on the Python side #23

Open MothNik opened 5 months ago

MothNik commented 5 months ago

The Cython compiler options in this line

# cython: language_level=3, boundscheck=True, wraparound=False, cdivision=True

are counterproductive when the following checks are made here

        try:
            return penta_solver1(mat_flat, rhs)
        except ZeroDivisionError:
            warnings.warn("pentapy: PTRANS-I not suitable for input-matrix.")
            return np.full_like(rhs, np.nan)

because cdivision=True explicitly prevents the compiled C-code from ever raising a ZeroDivisionError. This branch of the code was always dead. This is also indicated by the coverage which shows that these lines was never hit. Actually, the code produced nan-values naturally because a single nan will propagate to contaminate the full solution, but the error handling was never triggered.

For the warning to be issued, cdivision=False is required here. Even though the Cython file will be relatively yellow, there is no need to hesitate when setting cdivision=False. Cython does a very good job when it comes to handling the more elaborate division because it helps the branch predictor by making the zero-case the unlikely case. With this, the performance impact is negligible:

ZeroDivision

I will fix this in the same go as #11

MuellerSeb commented 3 months ago

Good catch. If it's not to cost intensive, this seems reasonable to me.