CBLAS/LAPACKE extension for 64bit integer

[mkrainiuk]

Intro

CBLAS/LAPACKE wrappers support 32bit integer and 64bit integer with the same function names located in different libraries but this approach does not allow mixing both libraries in one environment because of the symbols conflict.

This FR is for discussing potential solution to avoid this problem and have stable work in any environment for the applications/libraries that has dependency on CBLAS/LAPACKE built with 64bit integer or any other projects with similar API.

Problem

Similar to Fortran API where integer with no size specified is used, CBLAS/LAPACKE API integer type can be selected by special build option during compilation.
Example for CBLAS:

lapack/CBLAS/include/cblas.h

Lines 20 to 24 in 655e588

	#ifdef WeirdNEC
	#define CBLAS_INT int64_t
	#else
	#define CBLAS_INT int32_t
	#endif

Since C/Fortran symbols do not reflect the data type in the name in general (like C++), linker cannot detect if the resolved symbols has correct Integer type. As the result symbols mismatch could cause unexpected behavior, like incorrect results, data corruption, and segmentation fault.

Example:

A plugin "X" uses CBLAS API with 64bit integer type because it works with number of elements that does not fit to 32bit. Another plugin "Y" loads "libcblas.so" CBLAS library built with 32bit integer as a dependency in the complex application environment, As first loaded "libcblas.so" library will be used for the first plugin "X" which will cause segmentation fault during application execution.

Related Work

• Julia Programming Language is working with other math libraries in order to allow mixing API for different integer types in one environment: ABI conflicts due to 64-bit libopenblas.so JuliaLang/julia#4923
• OpenBLAS introduced SYMBOLSUFFIX build option based on the Julia request: Standardize ILP64 SONAME and symbols suffix OpenMathLib/OpenBLAS#646
• BLIS has good discussion for naming but it's waiting for standardization of the approach (which this FR attempts to do) to implement it in the library: BLIS should allow simultaneously exporting both 32- and 64-bit variants of BLAS/CBLAS flame/blis#43
• FlexiBLAS also has good analysis of the naming, converted to internal issue: Suffixed 64-bit integer symbols names mpimd-csc/flexiblas#12
• Intel oneMKL added Fortran-like BLAS C API extension with _64 suffix for 64bit integer to default LP64 interface library, _64 extension for Fortran-like LAPACK C API and CBLAS/LAPACKE API (link)
• cuBLAS 12.0 extends the cuBLAS API to support 64-bit integer, extended API has a _64 suffix in the name (link)
• Some work for libraries naming was already done for Reference-LAPACK: Allow a possible installation of index-64 library alongside standard index-32 library? #461
• Apple is using a different convention in Accelerate for LAPACK ILP64. Julia Computing wrappers for those routines on macOS: Add support for ILP64 Accelerate JuliaLinearAlgebra/libblastrampoline#113

Proposal

• Define a suffix for API with 64bit integer support so that users can control integer type and right API call in the code and avoid any runtime conflicts

  • 64bit integer API was selected to be extended with the suffix because 32bit integer API is default one for the most of the libraries. So it could be easier to change only code where 64bit integer API is required, and keep default API unchanged.

• Suffix considerations

  • Julia uses OpenBLAS with 64_ suffix for all C and Fortran symbols.
    Please note, from the API perspective it will be two different suffixes for C and Fortran APIs:
    • _64 for Fortran API, because e.g., snrm2/snrm2_64 functions will be converted to snrm2_/snrm2_64_ symbols (default gfortran/ifort compilers behavior on Unix)
    • 64_ for C API in order to get same symbol suffix, because e.g. cblas_snrm2/cblas_snrm264_ functions will be converted to cblas_snrm2/cblas_snrm264_ symbol (default gcc/icc/clang compilers behavior)
  • FlexiBLAS considered adding _64, it will improve readability for function names with numbers, e.g.: for API DGGES3 it is better to have DGGES3_64 instead of DGGES364.
  • Intel oneMKL and NVIDIA cuBLAS use _64 suffix

• Proposed suffix for CBLAS/LAPACK API: _64
  Please note, _64 must be at the end of the function name for the symbol identification simplicity, including cases like LAPACKE_dgeqrf_work, ILP64 API will be LAPACKE_dgeqrf_work_64, not LAPACKE_dgeqrf_64_work

Example:

// Current API
float      cblas_snrm2(const CBLAS_INT N, const float *X, const CBLAS_INT incX);
lapack_int LAPACKE_dgeev( int matrix_layout, char jobvl, char jobvr, lapack_int n, double* a, lapack_int lda, ... );

// 64bit fixed width integer API
float   cblas_snrm2_64(const int64_t N, const float *X, const int64_t incX); 
int64_t LAPACKE_dgeev_64( int matrix_layout, char jobvl, char jobvr, int64_t n, double* a, int64_t lda, ... );

int main() {
int64_t N, incX;
int32_t n, lda, ...;
...
cblas_snrm2_64(N, X, incX);
LAPACKE_dgeev(matrix_layout, jobvl, jobv, n, a, lda, ...);
...
}

Other Considerations

• Update APIs that do not have dependency on the Integer type, e.g.

  lapack/CBLAS/include/cblas.h

  Line 134 in 655e588

  void cblas_srotg(float *a, float *b, float *c, float *s);

• Do we need to consider other data types that could benefit from this approach?
• Add function mangling to the header so that users can convert names for all functions in the application during compilation, same way as currently it's supported for integer type: 32bit by default, 64bit if special macro is defined.
• Extend Fortran API. It's partially should be done anyway in order to enable extensions for CBLAS/LAPACKE wrappers

[mkrainiuk]

Summon people participated in the discussions mentioned in the Related Work section:
Julia @staticfloat @ViralBShah
FlexiBLAS @grisuthedragon
BLIS @devinamatthews
OpenBLAS @xianyi

[ViralBShah]

For Julia: @amontoison
For OpenBLAS: @martin-frbg

Also related: #824

[martin-frbg]

Yes, I'm lurking... OpenBLAS currently does it by running a big objcopy after the initial build, which is far from ideal. And standardization of the expected ILP64 suffix across (at least) Julia/Go/(Num|Sci)Py is obviously desirable.

[grisuthedragon]

And standardization of the expected ILP64 suffix across (at least) Julia/Go/(Num|Sci)Py is obviously desirable.

Ok, I compiled a bit what @mkrainiuk proposed, some of us already discussed in mpimd-csc/flexiblas#12, and what I discussed with a colleague.

Let us distinguish between suffixed and not suffixed. First, we have the ones without suffixes:

• libblas.so - This contains the 32-bit integer build with no suffixes. This is the default name and maintains compatibility with old codes, makefiles, build systems. So no one with old code will be affected by the change.
• libblas64.so - This contains the 64-bit integer build without suffixes.

I would like to introduce the 32-bit and 64-bit builds with suffixes. This could be for example

• libblas_32.so - The 32-bit integer build with suffixes.
• libblas_64.so - The 64-bit integer build with suffixes.

Regarding the rules for suffixes for each symbol, my suggestion is this. We look at all symbols from the point of view of their native interface. That is, before the compiler does any name mangling. This means that the normal BLAS and LAPACK routines are treated in the Fortran context, and CBLAS and LAPACK from a C perspective. This keeps the approach invariant even under strange compilers and ABI definition, like using IBM/Sun compilers or macOS(like @ViralBShah reported there: mpimd-csc/flexiblas#12 (comment) )

As for the suffixes themselves, I thought of the following:

• BLAS/LAPACK: DGEMM -> DGEMM_32 and DGEMM_64. The addition of the underscore is important so that routines like DGEQP3 do not end in DGEQP332 and DGEQP364.
• CBLAS: cblas_dgemm -> cblas_dgemm_32 and cblas_dgemm_64
• LAPACKE: LAPACKE_dgeqrf -> LAPACKE_dgeqrf_64. But in case of the _work routines, we could either go to LAPACKE_GEQRF_WORK -> LAPACKE_GEQRF_64_WORK or LAPACKE_GEQRF_WORK -> LAPACKE_GEQRF_WORK_64. I prefer the latter for simplicity when writing macros.

[martin-frbg]

I do not like the additional complexity of having a libblas64.so with non-suffixed symbols, but I guess it is inevitable for backwards compatibility ? Fully agreed on the leading underscore, and on not treating the "work" interfaces as something special - after all, they are (or may be) using ILP64 internally themselves, not just doing some work for another function that happens to be ILP64.
However, do we really need the _32 suffixes too, or could we just assume that any non-suffixed symbol is an implied _32 ? (I have to admit that enforcing suffixes would solve a potential problem with purely internal functions in OpenBLAS without having to mark them as hidden)

Cc @rgommers for NumPy and SciPy as the discussion seems to be gravitating towards here (we had something similar planned on a smaller scope) Maybe @kortschak or @vladimir-ch for Gonum if there is interest ?

[kortschak]

Gonum has a focus on floats, so I'm not sure how much interest there is in this from us.

[martin-frbg]

Thanks - discussion is entirely about large array addressing not fixed-point math.

[grisuthedragon]

I do not like the additional complexity of having a libblas64.so with non-suffixed symbols, but I guess it is inevitable for backwards compatibility ?

That allows to recompile projects with -fdefault-integer-8 or -i8 without changing any code (as long as it is purely Fortran)

Fully agreed on the leading underscore, and on not treating the "work" interfaces as something special - after all, they are (or may be) using ILP64 internally themselves, not just doing some work for another function that happens to be ILP64. However, do we really need the _32 suffixes too, or could we just assume that any non-suffixed symbol is an implied _32 ? (I have to admit that enforcing suffixes would solve a potential problem with purely internal functions in OpenBLAS without having to mark them as hidden)

The _32 suffixes is only for completeness, so I have no problem to neglect this variant.

[grisuthedragon]

@mkrainiuk
I saw you added the first steps towards the _64 in 4c2236a. Doing some tests show up with the following problems:

• the _64 is built although BUILD_INDEX64=OFF
• the _64 is built together with the symbols without the suffix. Thus, linking 32 bit and 64 bit interfaces is impossible again.
• I do not know if the way through the preprocessor and the command line is the best. Let's consider the following situation: We have a maximum command line length between 32kB and 2MB. In the case of 32kb and around 2000 BLAS/LAPACK routines, we have 16 characters for each definition in average. Having -DDGEQP3=DGEQP3_64 we get 19 characters (include a whitespace to separate them. This leads to strange situations in some platforms. I think a better way would be to provide a header-file containing the translation.

In general, it looks as it fits the suggestions we made here.

[mkrainiuk]

Thanks for the great feedback!

Hi @grisuthedragon

• the _64 is built although BUILD_INDEX64=OFF

BUILD_INDEX64 and BUILD_INDEX64_EXT_API are completely independent options:

• The main idea of BUILD_INDEX64_EXT_API option is to build additional set of symbols with "_64" suffix for ILP64 only that can be mixed with standard API for LP64 in one application or in one library.
• BUILD_INDEX64 option changes the integer type from 32-bit to 64-bit for the standard API and the library name from cblas to cblas64, but in some cases having different library names can't guarantee that there won't be symbol conflicts when both libraries are loaded to the same env. So I guess the long term solution could be migrating to "_64" API for ILP64 and dropping BUILD_INDEX64 option completely so that the standard API will be always for LP64.

• the _64 is built together with the symbols without the suffix. Thus, linking 32 bit and 64 bit interfaces is impossible again.

It's possible when the standard API is built for LP64 and the extended _64 API is built for ILP64, but it requires some code modifications for using _64 in an application:

int32_t m1, n1, lda1, incx1, incy1;
int64_t m2, n2, lda2, incx2, incy2;
...
cblas_dgemm(CblasColMajor, CblasNoTrans, m1, n1, ... );     //LP64 API
cblas_dgemm_64(CblasColMajor, CblasNoTrans, m2, n2, ...);   //ILP64 API

• I do not know if the way through the preprocessor and the command line is the best. Let's consider the following situation: We have a maximum command line length between 32kB and 2MB. In the case of 32kb and around 2000 BLAS/LAPACK routines, we have 16 characters for each definition in average. Having -DDGEQP3=DGEQP3_64 we get 19 characters (include a whitespace to separate them. This leads to strange situations in some platforms. I think a better way would be to provide a header-file containing the translation.

This is a good point, I considered new header file, but it requires manual update for any new function, so for CBLAS since the number of functions is relatively small I decided to generate the macro on the fly instead. I agree that this approach won't work for LAPACKE, so I'm looking for another solution that also can make this Fortran symbol renaming automatically during the build.

[grisuthedragon]

@mkrainiuk

BUILD_INDEX64 option changes the integer type from 32-bit to 64-bit for the standard API and the library name from cblas to cblas64, but in some cases having different library names can't guarantee that there won't be symbol conflicts when both libraries are loaded to the same env. So I guess the long term solution could be migrating to "_64" API for ILP64 and dropping BUILD_INDEX64 option completely so that the standard API will be always for LP64.

From my opinion, as mentioned before, we need three build variants. The standard build, without any special options, suffixes or similar one, that leads to the LP64 variant, as you mentioned. Then an ILP64 built without suffixes, (and 64 or _ilp64 added to the library name). This is required to rebuilt applications with -i8 or -fdefault-integer8 without code changes. And finally, the ilp64 build with suffixes, that, as you mentioned, should be the standard for ilp64 usage.
As far as I know from Julia (@ViralBShah) they are using the suffixed ilp64 library internally and thus it makes not problem to load other code that is linked against a LP64 BLAS. From this point of view, having the non-suffixed and the suffixed symbols in one library would be bad.

but in some cases having different library names can't guarantee that there won't be symbol conflicts when both libraries are loaded to the same env.

That's something nobody can guarantee and it is up to the programmer.

This is a good point, I considered new header file, but it requires manual update for any new function, so for CBLAS since the number of functions is relatively small I decided to generate the macro on the fly instead. I agree that this approach won't work for LAPACKE, so I'm looking for another solution that also can make this Fortran symbol renaming automatically during the build.

Since the set of functions changes only slowly, one can provide such a header file and eventually a script, which generates a new header from all the sources. One advantage of such a header file would be that one can provide it to the user as well to allow an easy migration between LP64 and ILP64 mode.

[langou]

I agree with @martin-frbg that the set of functions in LAPACK and LAPACKE changes slowly. For LAPACKE, at this point, there is no script that I am aware of and we are writing the LAPACKE layer functions by hand. This is possible because we only add a few routines at each release. Adding a layer related to LP64 / ILP64 / etc variants by hand would not be too much of an overkill. There is an exponential growth here though, but that would work. I am not asking for more work by hand but we are already generating S, C, D, Z by hand, and while not ideal that work-ish. The point is that LAPACK is slowly growing. @martin-frbg is correct.

[martin-frbg]

Different Martin but of course I agree with him :)