Completion of g77 ABI wrappers for MacOSX

[michaelwimmer]

The MacOSX Accelerate framework adheres to g77 ABI conventions (this is stated for example here, although the solution presented there is incomplete): This implies that

• every function returning a single precision floating point value (Fortran REAL) in fact returns a double precision
  value (Fortran DOUBLE PRECISION)
• every function returning a complex number is in fact internally a subroutine with the result as the first argument, i.e. COMPLEX FUNCTION X(A, B, C) is in fact SUBROUTINE X(RESULT, A, B, C)

Trying to link gfortran generated code with this different ABI results either in crashes (typically for complex) or wrong results (especially for REAL). Scipy already contained wrappers that solved these ABI issues for some BLAS/LAPACK routines (in particular C,ZDOTC, C,ZDOTU), but these ABI wrappers were rather incomplete.

As a result, the scipy tests involving single precision failed often on OSX machines (see e.g. issue #2248).

In this pull request:

• the ABI wrappers for the basic lapack/blas wrappers in scipy.linalg and scipy.lib are completed
• the ABI wrappers for the Fortran routines in scipy.sparse.linalg.isolve are completed
• test cases for single precision for the iterative solvers are added (there were no before, so the failure of the iterative solvers on OSX was unnoticed)
• the ABI wrappers for ARPACK are completed
• single precision in ARPACK has been enabled again
• pv's patch for ARPACK is included (see http://forge.scilab.org/index.php/p/arpack-ng/issues/1259/)

With these changes, all relevant tests pass again for my MacOSX 10.8 machine that showed many failures before these patches.

[michaelwimmer]

The Travis-CI failure is actually interesting. Note that it happens for one matrix in the non-symmetric solvers in both single and double precision. It only happens for the python 2.6 build, presumably because the random matrix tested is different in all builds.

Is building in Travis-CI determinsitic, i.e. would I get the same failure in a second build?

I have since a short while a suspicion that there is another bug in ARPACK, and it would be helpful to get the failing matrix ...

Anyways, the failure has nothing to do with my changes, as far as I can tell. It would have been there also before.

[pv]

I can reproduce the Travis CI failure, on Python 2.6 on Ubuntu. It doesn't occur on master branch, so it's probably caused by one of the changes here. It does not occur on Python 2.7, so something is subtly going wrong.

(Future advice: avoid merging the master branch to your own branch while working on a feature, and it's better to have one PR per one feature added.)

[michaelwimmer]

Thanks, pv, for the input. Can you tell me what numpy version you are using? That should enable me to reproduce that error on my ubuntu laptop.

Is the error deterministic for you (i.e. reproducible on successive runs)?

About the format: should I disentangle the PR?

[pv]

Numpy version was 1.7.1 for both. The error is deterministic.

The matrix itself is here: https://dl.dropboxusercontent.com/u/5453551/eigs-fail.npz
However, for this matrix, eigs gives sometimes ArpackNoConvergence. This seems to imply that things depend on the starting vector --- which in this case is determined internally by Arpack. Indeed, if you run the test alone, it doesn't fail (the other tests change Arpack's random generator internal state). I don't yet understand why the Python version matters, and why the changes in this PR would matter.

EDIT: ok, the order the tests are run in is different for different Python versions, so the starting vectors should also differ. The eigenvalues returned here are correct, but not those satisfying the which='LR' condition. The issue may just be that Arpack, being a Krylov method, cannot give perfect guarantees that it converges to the right set of eigenvalues. If so, the test should fix the starting vector v0 for this combination of which and input matrix.

About the PR: no reorganization is needed for now, I can sort it out when merging.

[rgommers]

This sounds promising!

Is the expansion of .c.src templates needed? I haven't checked in detail if something is different for the different typecodes, but at first sight it seems like a lot of duplication.

[rgommers]

Works for me on OS X 10.8 (without -ff2c).

[rgommers]

One error on OS X 10.6 (repeated 7 times):

======================================================================
ERROR: Failure: ImportError (dlopen(/Users/rgommers/Code/scipy/scipy/sparse/linalg/isolve/_iterative.so, 2):     Symbol not found: _clanhf_
  Referenced from: /Users/rgommers/Code/scipy/scipy/sparse/linalg/isolve/_iterative.so
  Expected in: dynamic lookup
)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/nose-1.1.2-py2.6.egg/nose/loader.py", line 390, in loadTestsFromName
    addr.filename, addr.module)
  File "/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/nose-1.1.2-py2.6.egg/nose/importer.py", line 39, in importFromPath
    return self.importFromDir(dir_path, fqname)
  File "/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/nose-1.1.2-py2.6.egg/nose/importer.py", line 86, in importFromDir
    mod = load_module(part_fqname, fh, filename, desc)
  File "/Users/rgommers/Code/scipy/scipy/linalg/tests/test_interpolative.py", line 33, in <module>
    from scipy.sparse.linalg import aslinearoperator
  File "/Users/rgommers/Code/scipy/scipy/sparse/linalg/__init__.py", line 108, in <module>
    from .isolve import *
  File "/Users/rgommers/Code/scipy/scipy/sparse/linalg/isolve/__init__.py", line 6, in <module>
    from .iterative import *
  File "/Users/rgommers/Code/scipy/scipy/sparse/linalg/isolve/iterative.py", line 7, in <module>
    from . import _iterative
ImportError: dlopen(/Users/rgommers/Code/scipy/scipy/sparse/linalg/isolve/_iterative.so, 2): Symbol not found: _clanhf_
  Referenced from: /Users/rgommers/Code/scipy/scipy/sparse/linalg/isolve/_iterative.so
  Expected in: dynamic lookup

[michaelwimmer]

@pv: Indeed, it seems we hit a very unlucky seed. In fact, but repeatedly calling eigs() a few hundred times on the matrix you posted, I eventually do get the same (wrong) result even for python 2.7 and vanilla scipy 0.12 as shipped in Ubuntu 12.04. As you noted, the eigenvalues/vectors are correct, it's just that the method didn't converge to the one desired.

Fixing the starting vector is one way of course, although it wouldn't prevent such a situation with an unlucky seed, for example when the random test matrix changes, for example due to changes in numpy, etc.

[michaelwimmer]

@rgommers:

• why I removed the *.c.src files: no particular reason, I simply fixed scipy.linalg first (which didn't have *c.src files),
  and then copied to scipy.lib.
  It would probably make sense anyways to have one instance of the wrapper files at a central position. The
  problem was that I simply didn't know where to put it, and how to make sure that this central wrapper is
  compiled before the modules that need it. (I simply do not understand the build process very well). Any ideas in
  that direction?
• Good catch on OSX 10.6: clanhf was only introduced in Lapack 3.2, seems 10.6 has an older version. This is easy
  to fix, it is actually not used at all (I created the wrapper file using a script from Lapack, it contains all possibly
  problematic functions, regardless if they are used by scipy or not).
  I'll fix it once I heard from you an advice about the possibility of a central wrapper (see above).

[pv]

@michaelwimmer: Numpy is committed to keeping the random sequence (for a given seed) the same, so a test with a fixed starting vector v0 should be well-defined.

[michaelwimmer]

@pv : An alternative to a fixed starting vector would be to retry the computation if the eigenvalues and -vectors are in principle correct (as is the case here), but the eigenvalues are not the desired ones. eval_evec() could be easily changed to do that.

The problem I see with a fixed starting vector is that if even with numpy having the same random numbers, if different python versions have a different order of tests (as you said was the case with python 2.6), we could run into the same problem in the future.

[pv]

@michaelwimmer: the different order of tests affects only Arpack's internal random generator, which I believe is not used if a starting vector is specified.

However, retrying the computation is also an OK solution, and might be better.

[pv]

@rgommers: does it work now on 10.6?

The wrappers under scipy.lib are deprecated and will be removed at some point, so this part of the duplication will go away. There are some incompatibilities between them and those under scipy.linalg; I'm not sure if this PR changes some of the scipy.lib routine signatures.

One possible place to put the unified wrappers could be under scipy/_build_tools/src and add a function to scipy/_build_tools/__init__.py for retrieving the file names of the source files in question. This can be done later in a separate PR, though.

[rgommers]

This now works for 10.6 and 10.8. Haven't been able to check 10.5, but I guess it should work there too. If not we'll find out soon enough.

[michaelwimmer]

Sorry, was away for a few days. Unfortunately I'm also unable to find a machine with 10.5 to test on ...

So, is there anything left for me to do in this PR?

@pv about the unified wrappers: Right now, the wz/cdotc/u wrappers for ARPACK and iterative are functions (as expected from the fortran routines), whereas the wrappers of the same name in linalg and lib are routines. This also has to be "unified". But that's not a problem by itself.

[pv]

Looks ready to go to me. I also looked through the f2py changes and didn't see anything out of place.

[pv]

Rebased and merged in 62883c9. Thanks a lot!

[sturlamolden]

This does not look like g77 ABI:

https://developer.apple.com/library/mac/documentation/Accelerate/Reference/BLAS_Ref/Reference/reference.html#//apple_ref/c/func/SDOT

[sturlamolden]

From vBLAS.h

/*
   =================================================================================================
   Prototypes for FORTRAN BLAS
   ===========================
   These are prototypes for the FORTRAN callable BLAS functions.  They are implemented in C for
   Mac OS, as thin shims that simply call the C BLAS counterpart.  These routines should never be
   called from C, but need to be included here so they will get output for the stub library.  It
   won't hurt to call them from C, but who would want to since you can't pass literals for sizes?
   FORTRAN compilers are typically MPW tools and use PPCLink, so they will link with the official
   vecLib stub from Apple.
   =================================================================================================
*/
/*
 *  SDOT()
 *  
 *  Availability:
 *    Mac OS X:         in version 10.0 and later in vecLib.framework
 *    CarbonLib:        not in Carbon, but vecLib is compatible with CarbonLib
 *    Non-Carbon CFM:   in vecLib 1.0.2 and later
 */
extern float 
SDOT(
  const int *    /* N */,
  const float *  /* X */,
  const int *    /* incX */,
  const float *  /* Y */,
  const int *    /* incY */)                                        AVAILABLE_MAC_OS_X_VERSION_10_0_AND_LATER;

This is not g77 ABI.

It seems this PR duplicates a CBLAS wrapper already present in Accelerate.

[michaelwimmer]

@sturlamolden: Well, Apple was in fact even inconsistent with their own design choices. LAPACK all adheres to the g77 convention. BLAS is a bit of a mixture, with SDOT actually not adhering to this convention. (If you google a bit for it, this was reported elsewhere). However, this is all taken into account in the wrappers now, to my knowledge.

[sturlamolden]

Testing on OSX 10.9 is seems the C prototype for SDOT in vBLAS.h is incorrect. In reality SDOT returns a double with -m64.