Permalink
Browse files

BUG: solve segfaults in sparse.linalg.isolve on OS X. Closes #1321 and

…#1472.

The problem is that Accelerate is clapack but has flapack names.  See #725 for
more detailed description.
  • Loading branch information...
1 parent 1479811 commit b75d42e6930dd43c35438b58063b0b2212544ddd @rgommers rgommers committed Aug 18, 2011
@@ -1,3 +1,4 @@
+# vim:syntax=python
from os.path import join as pjoin
from numscons import GetNumpyEnvironment
@@ -6,7 +6,7 @@ from os.path import join as pjoin, splitext
from numscons import GetNumpyEnvironment
from numscons import CheckF77LAPACK
-from numscons import write_info
+from numscons import write_info, IsAccelerate, IsVeclib
env = GetNumpyEnvironment(ARGUMENTS)
env.Tool('f2py')
@@ -53,4 +53,21 @@ for method in raw_sources:
res = env.FromFTemplate(target, pjoin('iterative', method))
sources.append(res[0])
-env.NumpyPythonExtension('_iterative', source = sources)
+#--------------------------------------------------
+# BLAS wrapper to fix ABI incompatibilities on OS X
+# -------------------------------------------------
+use_c_calling = IsAccelerate(env, "lapack") or IsVeclib(env, "lapack")
+
+if use_c_calling:
+ blas_wrapper = ['veclib_cabi_c.c', 'veclib_cabi_f.f']
+else:
+ blas_wrapper = ['dummy.f']
+
+blas_wrapper = [pjoin('iterative', 'FWRAPPERS', wrapper) \
+ for wrapper in blas_wrapper]
+
+env.AppendUnique(LIBPATH = ['.'])
+veclibwrap_lib = env.DistutilsStaticExtLibrary('_veclibwrap', source=blas_wrapper)
+env.Prepend(LIBS='_veclibwrap')
+
+env.NumpyPythonExtension('_iterative', source=sources)
@@ -113,7 +113,7 @@
<rt> TOL, BNRM2, RHOTOL,
$ <sdsd=s,d,s,d>GETBREAK,
$ <rc=s,d,sc,dz>NRM2
- <_t> ALPHA, BETA, RHO, RHO1, <xdot=sdot,ddot,cdotc,zdotc>
+ <_t> ALPHA, BETA, RHO, RHO1, <xdot=sdot,ddot,wcdotc,wzdotc>
*
* indicates where to resume from. Only valid when IJOB = 2!
INTEGER RLBL
@@ -117,7 +117,7 @@
$ RHOTOL, OMEGATOL, <sdsd=s,d,s,d>GETBREAK,
$ <rc=s,d,sc,dz>NRM2
<_t> ALPHA, BETA, RHO, RHO1, OMEGA, TMPVAL,
- $ <xdot=sdot,ddot,cdotc,zdotc>
+ $ <xdot=sdot,ddot,wcdotc,wzdotc>
* indicates where to resume from. Only valid when IJOB = 2!
INTEGER RLBL
*
@@ -100,7 +100,7 @@
* .. Local Scalars ..
INTEGER MAXIT, R, Z, P, Q, NEED1, NEED2
<_t> ALPHA, BETA, RHO, RHO1,
- $ <xdot=sdot,ddot,cdotc,zdotc>
+ $ <xdot=sdot,ddot,wcdotc,wzdotc>
<rt> <rc=s,d,sc,dz>NRM2, TOL
*
* indicates where to resume from. Only valid when IJOB = 2!
@@ -118,7 +118,7 @@
$ <rc=s,d,sc,dz>NRM2
<_t> ALPHA, BETA, RHO, RHO1, TMPVAL,
- $ <xdot=sdot,ddot,cdotc,zdotc>
+ $ <xdot=sdot,ddot,wcdotc,wzdotc>
* ..
* indicates where to resume from. Only valid when IJOB = 2!
INTEGER RLBL
@@ -0,0 +1,57 @@
+ double complex function wzdotc(n, zx, incx, zy, incy)
+ double complex zx(*), zy(*), z
+ double complex zdotc
+ integer n, incx, incy
+
+ z = zdotc(n, zx, incx, zy, incy)
+ wzdotc = z
+
+ end
+
+ double complex function wzdotu(n, zx, incx, zy, incy)
+ double complex zx(*), zy(*), z, zdotu
+ integer n, incx, incy
+
+ z = zdotu(n, zx, incx, zy, incy)
+ wzdotu = z
+
+ return
+ end
+
+ complex function wcdotc(n, cx, incx, cy, incy)
+ complex cx(*), cy(*), c, cdotc
+ integer n, incx, incy
+
+ c = cdotc(n, cx, incx, cy, incy)
+ wcdotc = c
+
+ return
+ end
+
+ complex function wcdotu(n, cx, incx, cy, incy)
+ complex cx(*), cy(*), c, cdotu
+ integer n, incx, incy
+
+ c = cdotu(n, cx, incx, cy, incy)
+ wcdotu = c
+
+ return
+ end
+
+ complex function wcladiv(x, y)
+ complex x, y, z
+ complex cladiv
+
+ z = cladiv(x, y)
+ wcladiv = z
+ return
+ end
+
+ double complex function wzladiv(x, y)
+ double complex x, y, z
+ double complex zladiv
+
+ z = zladiv(x, y)
+ wzladiv = z
+ return
+ end
@@ -0,0 +1,26 @@
+#include <complex.h>
+#include <vecLib/vecLib.h>
+
+#define WRAP_F77(a) a##_
+void WRAP_F77(veclib_cdotc)(const int *N, const complex float *X, const int *incX,
+const complex float *Y, const int *incY, complex float *dotc)
+{
+ cblas_cdotc_sub(*N, X, *incX, Y, *incY, dotc);
+}
+
+void WRAP_F77(veclib_cdotu)(const int *N, const complex float *X, const int *incX,
+const complex float *Y, const int *incY, complex float* dotu)
+{
+ cblas_cdotu_sub(*N, X, *incX, Y, *incY, dotu);
+}
+
+void WRAP_F77(veclib_zdotc)(const int *N, const double complex *X, const int
+*incX, const double complex *Y, const int *incY, double complex *dotu)
+{
+ cblas_zdotc_sub(*N, X, *incX, Y, *incY, dotu);
+}
+void WRAP_F77(veclib_zdotu)(const int *N, const double complex *X, const int
+*incX, const double complex *Y, const int *incY, double complex *dotu)
+{
+ cblas_zdotu_sub(*N, X, *incX, Y, *incY, dotu);
+}
@@ -0,0 +1,55 @@
+ double complex function wzdotc(n, zx, incx, zy, incy)
+ double complex zx(*), zy(*), z
+ integer n, incx, incy
+
+ call veclib_zdotc(n, zx, incx, zy, incy, z)
+
+ wzdotc = z
+ return
+ end
+
+ double complex function wzdotu(n, zx, incx, zy, incy)
+ double complex zx(*), zy(*), z
+ integer n, incx, incy
+
+ call veclib_zdotu(n, zx, incx, zy, incy, z)
+
+ wzdotu = z
+ return
+ end
+
+ complex function wcdotc(n, cx, incx, cy, incy)
+ complex cx(*), cy(*), c
+ integer n, incx, incy
+
+ call veclib_cdotc(n, cx, incx, cy, incy, c)
+
+ wcdotc = c
+ return
+ end
+
+ complex function wcdotu(n, cx, incx, cy, incy)
+ complex cx(*), cy(*), c
+ integer n, incx, incy
+
+ call veclib_cdotu(n, cx, incx, cy, incy, c)
+
+ wcdotu = c
+ return
+ end
+
+ complex function wcladiv(x, y)
+ complex x, y, z
+
+ call cladiv(z, x, y)
+ wcladiv = z
+ return
+ end
+
+ double complex function wzladiv(x, y)
+ double complex x, y, z
+
+ call zladiv(z, x, y)
+ wzladiv = z
+ return
+ end
@@ -104,7 +104,7 @@
* .. Local Scalars ..
INTEGER I, MAXIT, AV, GIV, H, R, S, V, W, Y,
$ NEED1, NEED2
- <_t> <xdot=sdot,ddot,cdotc,zdotc>
+ <_t> <xdot=sdot,ddot,wcdotc,wzdotc>
<_t> toz
<_t> TMPVAL
<rt> BNRM2, RNORM, TOL,
@@ -416,7 +416,7 @@
<rt=real,double precision,real,double precision>
$ <rc=s,d,sc,dz>NRM2, ONE
PARAMETER ( ONE = 1.0D+0 )
- <_t> <xdot=sdot,ddot,cdotc,zdotc>
+ <_t> <xdot=sdot,ddot,wcdotc,wzdotc>
<_t> TMPVAL
EXTERNAL <_c>AXPY, <_c>COPY, <xdot>, <rc>NRM2, <_c>SCAL
*
@@ -129,7 +129,7 @@
<_t> BETA, GAMMA, GAMMA1, DELTA, EPS, ETA, XI,
$ RHO, RHO1, THETA, THETA1, C1, TMPVAL,
- $ <xdot=sdot,ddot,cdotc,zdotc>,
+ $ <xdot=sdot,ddot,wcdotc,wzdotc>,
$ toz
*
* indicates where to resume from. Only valid when IJOB = 2!
@@ -7,6 +7,24 @@
from glob import glob
from os.path import join
+
+def needs_veclib_wrapper(info):
+ """Returns true if needs special veclib wrapper."""
+ import re
+ r_accel = re.compile("Accelerate")
+ r_vec = re.compile("vecLib")
+ res = False
+ try:
+ tmpstr = info['extra_link_args']
+ for i in tmpstr:
+ if r_accel.search(i) or r_vec.search(i):
+ res = True
+ except KeyError:
+ pass
+
+ return res
+
+
def configuration(parent_package='',top_path=None):
from numpy.distutils.system_info import get_info, NotFoundError
@@ -30,17 +48,28 @@ def configuration(parent_package='',top_path=None):
'QMRREVCOM.f.src',
# 'SORREVCOM.f.src'
]
+
+ if needs_veclib_wrapper(lapack_opt):
+ methods += [join('FWRAPPERS', 'veclib_cabi_f.f'),
+ join('FWRAPPERS', 'veclib_cabi_c.c')]
+ else:
+ methods += [join('FWRAPPERS', 'dummy.f')]
+
+
Util = ['STOPTEST2.f.src','getbreak.f.src']
sources = Util + methods + ['_iterative.pyf.src']
config.add_extension('_iterative',
- sources = [join('iterative',x) for x in sources],
- extra_info = lapack_opt
+ sources=[join('iterative', x) for x in sources],
+ extra_info=lapack_opt,
+ depends=[join('iterative', 'FWRAPPERS', x) for x in
+ ['veclib_cabi_f.f', 'veclib_cabi_c.c', 'dummy.f']]
)
config.add_data_dir('tests')
return config
+
if __name__ == '__main__':
from numpy.distutils.core import setup

0 comments on commit b75d42e

Please sign in to comment.