clang -fPIC blaswrap.c -o blaswrap.o -DUSE_BLASWRAP -c

blaswrap2.c

/*

Copyright (C) 2011 Jarno Rajahalme
Copyright (C) 2013 Jameson Nash

This is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the
Free Software Foundation; either version 2.1 of the License, or (at
your option) any later version.

This software is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this software; see the file COPYING.  If not, see
<http://www.gnu.org/licenses/>.


Wrapper for Apple libBLAS.dylib and libLAPACK.dylib

At least on the versions of OSX 10.6 so far (up and including 10.6.6)
these libraries are incompatible with 64 bit builds, as some functions
in libBLAS.dylib are not conforming to F2C calling conventions, as
they should.  This breaks them in 64-bit builds on the x86_64
architecture.

Newer gfortran compoilers no longer default to the F2C calling
convention.  These wrappers map the F2C conformant functions in
libBLAS and libLAPACK to the native gfortran calling convention, so
that the libraries can be used with software built for x86_64
architecture.

The wrappers could be made more efficient by calling the cblas functions
directly instead of redirecting through their fortran equivalents. But
that would require more typing on my part.
 
 */

#ifdef HAVE_CONFIG_H
#include <config.h> /* USE_BLASWRAP ? */
#endif

#ifdef USE_BLASWRAP

/*
 * Since this is a wrapper for fortran functions, we do not have prototypes for them.
 */
#pragma GCC diagnostic ignored "-Wmissing-prototypes"

#include <dlfcn.h>

void abort (void);

/*
 * Apple LAPACK follows F2C calling convention,
 * Convert to normal gfortran calling convention
 */

static void (*f2c_lapack_func[])(void); /* forward declaration for the wrapper */

/*
 * LAPACK Wrappers, only need to convert the return value from double to float
 */

typedef double (*F2C_CALL_0)(void);
typedef double (*F2C_CALL_1)(void *a1);
typedef double (*F2C_CALL_2)(void *a1, void *a2);
typedef double (*F2C_CALL_3)(void *a1, void *a2, void *a3);
typedef double (*F2C_CALL_4)(void *a1, void *a2, void *a3, void *a4);
typedef double (*F2C_CALL_5)(void *a1, void *a2, void *a3, void *a4, void *a5);
typedef double (*F2C_CALL_6)(void *a1, void *a2, void *a3, void *a4, void *a5, void *a6);
typedef double (*F2C_CALL_7)(void *a1, void *a2, void *a3, void *a4, void *a5, void *a6, void *a7);
typedef double (*F2C_CALL_8)(void *a1, void *a2, void *a3, void *a4, void *a5, void *a6, void *a7, void *a8);

#define F2C_LAPACK_CALL_8(name) \
  float name (void *a1, void *a2, void *a3, void *a4, void *a5, void *a6, void *a7, void *a8) \
  {                                                                     \
    return ((F2C_CALL_8)f2c_lapack_func[f2c_ ## name]) (a1, a2, a3, a4, a5, a6, a7, a8); \
  }

#define F2C_LAPACK_CALL_7(name) \
  float name (void *a1, void *a2, void *a3, void *a4, void *a5, void *a6, void *a7) \
  {                                                                     \
    return ((F2C_CALL_7)f2c_lapack_func[f2c_ ## name]) (a1, a2, a3, a4, a5, a6, a7); \
  }

#define F2C_LAPACK_CALL_6(name) \
  float name (void *a1, void *a2, void *a3, void *a4, void *a5, void *a6) \
  {                                                                     \
    return ((F2C_CALL_6)f2c_lapack_func[f2c_ ## name]) (a1, a2, a3, a4, a5, a6); \
  }

#define F2C_LAPACK_CALL_5(name) \
  float name (void *a1, void *a2, void *a3, void *a4, void *a5)         \
  {                                                                     \
    return ((F2C_CALL_5)f2c_lapack_func[f2c_ ## name]) (a1, a2, a3, a4, a5); \
  }

#define F2C_LAPACK_CALL_4(name) \
  float name (void *a1, void *a2, void *a3, void *a4)                   \
  {                                                                     \
    return ((F2C_CALL_4)f2c_lapack_func[f2c_ ## name]) (a1, a2, a3, a4); \
  }

#define F2C_LAPACK_CALL_3(name) \
  float name (void *a1, void *a2, void *a3)                          \
  {                                                                  \
    return ((F2C_CALL_3)f2c_lapack_func[f2c_ ## name]) (a1, a2, a3); \
  }

#define F2C_LAPACK_CALL_2(name) \
  float name (void *a1, void *a2)                                \
  {                                                              \
    return ((F2C_CALL_2)f2c_lapack_func[f2c_ ## name]) (a1, a2); \
  }

#define F2C_LAPACK_CALL_1(name) \
  float name (void *a1)                                      \
  {                                                          \
    return ((F2C_CALL_1)f2c_lapack_func[f2c_ ## name]) (a1); \
  }

#define F2C_LAPACK_CALL_0(name) \
  float name (void)                                        \
  {                                                        \
    return ((F2C_CALL_0)f2c_lapack_func[f2c_ ## name]) (); \
  }

#define F2C_LAPACK_CALL_NONE(name)

#define F2C_LAPACK_CALL(name, args) F2C_LAPACK_CALL_ ## args (name)

#define ENUM_ITEM(name, args)                   \
  f2c_ ## name, 

#define NAME_TO_STRING_CASE(name, args)         \
  case f2c_ ## name: return #name;

#define DEFINE_LAPACK_ENUM(name, list)  \
  typedef enum {                        \
    list(ENUM_ITEM)                     \
  } name;                               \
  static const char*                    \
  f2c_ ## name ## _name (name n) {      \
    switch (n) {                        \
      list(NAME_TO_STRING_CASE)         \
    default: return "";                 \
    }                                   \
  }                                     \
  list(F2C_LAPACK_CALL)

#define DEFINE_BLAS_ENUM(name, list)    \
  typedef enum {                        \
    list(ENUM_ITEM)                     \
  } name;                               \
  static const char*                    \
  f2c_ ## name ## _name(name n) {       \
    switch (n) {                        \
      list(NAME_TO_STRING_CASE)         \
    default: return "";                 \
    }                                   \
  }

/*
 * Lapack functions that need the return value converted from double to float
 * Simple stub is used for functions with all (0-6) args in registers,
 * individual stub types are used for longer argument counts.
 */
#define LAPACK_LIST(_)  \
  _(sasum_,3) \
  _(scasum_,3) \
  _(scnrm2_,3) \
  _(sdot_,5) \
  _(sdsdot_,6) \
  _(snrm2_,3) \
  _(LAPACK_COUNT,NONE)
//  _(scabs1,1)

/*
 * These need more complex wrappers, as arguments need to be handled
 */
#define BLAS_LIST(_)    \
  _(cdotu_,6)           \
  _(zdotu_,6)           \
  _(cdotc_,6)           \
  _(zdotc_,6)           \
  _(BLAS_COUNT,NONE)

DEFINE_BLAS_ENUM(blas, BLAS_LIST)

DEFINE_LAPACK_ENUM(lapack, LAPACK_LIST)

/*
 * Function pointer arrays, indexed by the enums
 */
static void (*f2c_blas_func[f2c_BLAS_COUNT])(void) = { 0 };
static void (*f2c_lapack_func[f2c_LAPACK_COUNT])(void) = { 0 };

/*
 * Initialization: This is called by ld before any of the wrappers are called
 * Get the function pointers to the wrapped functions in Apple vecLib
 */
typedef void (*BLASParamErrorProc)(
   const char *funcName,
   const char *paramName,
   const int *paramPos,
   const int *paramValue);
extern void SetBLASParamErrorProc (BLASParamErrorProc ErrorProc);
void BLASParamErrorProcNULL(
   const char *funcName,
   const char *paramName,
   const int *paramPos,
   const int *paramValue) { }

static void * apple_vecLib = 0;

__attribute__((constructor))
static void initVecLibWrappers (void)
{
  apple_vecLib = dlopen ("/System/Library/Frameworks/vecLib.framework/Versions/A/vecLib", RTLD_LOCAL | RTLD_NOLOAD | RTLD_FIRST);
  if (0 == apple_vecLib)
    abort ();

  int i;
  for (i = 0; i < f2c_LAPACK_COUNT; i++)
    if (0 == (f2c_lapack_func[i] = dlsym(apple_vecLib, f2c_lapack_name(i))))
      abort ();  
  for (i = 0; i < f2c_BLAS_COUNT; i++)
    if (0 == (f2c_blas_func[i] = dlsym(apple_vecLib, f2c_blas_name(i))))
      abort ();  
  SetBLASParamErrorProc( BLASParamErrorProcNULL );
}

__attribute__((destructor))
static void finiVecLibWrappers (void)
{
  if (apple_vecLib)
    dlclose (apple_vecLib);
  apple_vecLib = 0;
  SetBLASParamErrorProc(0);
}


/*
 * BLAS wrappers, F2C convention passes returned complex as an extra first
 * argument
 */
typedef struct { float r, i; } complex;
typedef struct { double r, i; } doublecomplex;

typedef void (*F2C_BLAS_CALL_6)(void *c, void *a1, void *a2, void *a3, void *a4, void *a5);

#define F2C_BLAS_CALL(type, name) \
type name (void *a1, void *a2, void *a3, void *a4, void *a5) \
{ \
  type cplx; \
  ((F2C_BLAS_CALL_6)f2c_blas_func[f2c_ ## name]) (&cplx, a1, a2, a3, a4, a5); \
  return cplx; \
}

F2C_BLAS_CALL(complex, cdotu_)
F2C_BLAS_CALL(doublecomplex, zdotu_)
F2C_BLAS_CALL(complex, cdotc_)
F2C_BLAS_CALL(doublecomplex, zdotc_)

#endif /* USE_BLASWRAP */