scipy.sparse.linalg.isolve depends on cdotc and zdotc, which require wrappers on Mac OS X (Trac #1321) #1847

scipy-gitbot opened this Issue Apr 25, 2013 · 6 comments


None yet

1 participant

Original ticket on 2010-10-28 by trac user commandant2k, assigned to @wnbell.

The FORTRAN reverse communication routines for iterative solution of linear systems in the scipy.sparse.linalg.isolve module each depend on the BLAS routines cdotc and zdotc to compute complex inner products in single and double precision, respectively. On Mac OS X installations that use the Accelerate or vecLib frameworks for BLAS and LAPACK, these routines are known to have issues when called from FORTRAN. The cause is an ABI mismatch between the functions compiled into the Mac OS X frameworks and the executable code produced by gfortran.

This problem exists in other modules of scipy, especially ARPACK in scipy.sparse.linalg.eigen. A workaround was implemented that defines local wrapper functions to call cblas_cdotc_sub and cblas_zdotc_sub, which provide subroutine versions of the BLAS functions that avoid the ABI issue. However, this accommodation was not made in scipy.sparse.linalg.isolve.

Attached is a patch to SVN revision 6856 that adds the wrapper functions to the isolve module in the same fashion as the wrappers from scipy.lib.blas and scipy.sparse.linalg.eigen. If neither the Accelerate nor vecLib frameworks are used, a dummy wrapper is in place to directly call the BLAS functions. The FORTRAN reverse communication routines have been modified to call the wrapper routine rather than the BLAS routines.

The patch has been tested on scipy 0.8.0 and SVN 4104bf3 on an Intel Macintosh. More testing should be done on PPC Macs and non-Macintosh systems to confirm that the dummy wrapper behaves as expected on other platforms.

Observe that the issues with cdotc and zdotc are only observed when attempting to solve linear systems involving single- or double-precision complex unknowns; for real-valued systems, the module does not invoke the problematic routines.

Attachment added by trac user commandant2k on 2010-10-28: scipy-r6856-isolve-veclib.patch

@pv wrote on 2011-08-18

The patch looks OK to me -- don't have a OSX at hand, so couldn't test it. Some code duplication with the stuff in ARPACK, but not too bad. The veclib_cabi_c.c.src however needs a #include <complex.h> on the top.

Milestone changed to 0.10.0 by @pv on 2011-08-18

@rgommers wrote on 2011-08-19

It doesn't apply cleanly, and doesn't work for me once I resolve the conflicts. I have the relevant parts (without duplication) in Your comment on gh-1997 seemed to favor trying to reuse the files in ARPACK/FWRAPPERS. Duplicate files would actually make the fix simpler and keep isolve and arpack fully independent.

@pv wrote on 2011-08-19

Avoiding code duplication would be better, but if it seems difficult to do given the number of build systems we have, then it's better just to copy the files over.

@rgommers wrote on 2011-08-25

Fixed in a similar way in b75d42e. Thanks a lot for this patch!

@patvarilly patvarilly pushed a commit to patvarilly/scipy that referenced this issue May 1, 2013
@pv pv BUG: special/specfun: fix bug in MTU12
Off-by-one error in the order up to which the Bessel functions were
computed (the arrays are accessed up to KM+1). Doesn't seem to affect
results much as the series seems mostly to converge before reaching the

Covered by test_mathieu_modsem2 (see Trac #1847)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment