[SciPy-user] Confusing BLAS/LAPACK situation on Mac OS X

[Arnd Baecker]

On Sat, 28 Jan 2006, Robert Kern wrote:

> Brian Granger wrote:
>
> > What about:
> >
> > /System/Library/Frameworks/Accelerate.framework/Frameworks/vecLib.framework/
> >
> > Versions/Current/libLAPACK.dylib
> >
> > and
> >
> > Headers/clapack.h
> >
> > Also, Apple's docs on Accelerate claim that they include clapack.
> >
> > Are these incomplete/unusable for some reason?
> > It seems silly to build ATLAS if it is already there.

It might also be a bad choice, as ATLAS can be slower than
other optimized LAPACK routines, e.g. MKL or ACML.
Do any benchmarks of `Accelerate` vs. ATLAS exist?

> Look at the symbols in libBLAS.dylib and libLAPACK.dylib. You can see the CBLAS
> versions in libBLAS.dylib, but there are only the FORTRAN versions in
> libLAPACK.dylib. You are right that the docs claim (correctly) that
> Accelerate.framework contains CLAPACK. However, I goofed and referred to things
> incorrectly. CLAPACK is simply an implementation of FORTRAN LAPACK in C. It uses
> FORTRAN column-major arrays. The flapack module links against these. However,
> ATLAS also provides row-major versions of some(?) LAPACK functions to match the
> row-major CBLAS routines. The clapack module wraps these. These are not included
> in Accelerate.framework, although the row-major CBLAS versions are.

If I understand things correctly, creating a matrix with (e.g)
  mat = zeros((Ny, Nx), "d")
gives row-major (i.e. C style) ordering.

Calling (for example) scipy.linalg.eig()
will determine the ordering via `scipy.linalg.basic.get_lapack_funcs`
and return the "Fortran code for leading array with column major order"
(i.e. flapack) and in all other cases clapack is preferred over flapack.

Now, clapack is not available, e.g. when
ATLAS is replaced by  Accelerate or MKL (presumably ACML as well),
or some other Fortran only Lapack implementation.
In such a case a copy of the original array will be made while
calling `eig` on the level of f2py (at least that's my experience so far -
Pearu, please correct me if this statement is too general).

Therefore, depending on the available LAPACK implementation,
a copy of the array will be made or not.
I am not sure if this is optimal, or if one
should start straight away with fortran type column major order.
For Numeric I usually used the trick
  mat = transpose(zeros((Nx, Ny), "d"))
- note the change order for Nx and Ny - and filled the matrix elements
just via
  mat[ny,nx] = ...

With numpy
  mat = zeros((Ny, Nx), fortran=1)
should do the job (I have not tested this yet).

If all the above is correct, then the solution I would use
for myself is to set `fortran=1` for all arrays which
will be used by some LAPACK routine.
Then no unnecessary copies of (presumably large) arrays will take
place on any machine and one could stop worrying about flapack
vs. clapack ;-).

Not sure if that is the solution for everyone - so I am happy
to learn about any drawbacks....

Just my 2cent,

Arnd