Search code examples
clapacklapacke

Using LAPACKE_zgetrs with LAPACK_ROW_MAJOR causes illegal memory access

I am trying to solve a linear system using the following code:

#include <stdio.h>
#include <lapacke.h>

int main () {
    lapack_complex_double mat[4];
    lapack_complex_double vec[2];
    lapack_int p[2];

    mat[0] = lapack_make_complex_double(1,0);
    mat[1] = lapack_make_complex_double(1,0);
    mat[2] = lapack_make_complex_double(1,0);
    mat[3] = lapack_make_complex_double(-1,0);

    vec[0] = lapack_make_complex_double(1,0);
    vec[1] = lapack_make_complex_double(1,0);

    LAPACKE_zgetrf(LAPACK_ROW_MAJOR, 2, 2, mat, 2, p);
    LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N', 2, 1, mat, 2, p, vec, 2);

    printf("%g %g\n", lapack_complex_double_real(vec[0]),
        lapack_complex_double_imag(vec[0]));
    return 0;
}

For some reasons, this causes illegal memory access in LAPACKE_zgetrs (as detected by valgrind and by my big program crashing in zgetrs because of "glibc detected corruption or double free"). I did not include this in my SSCCE for brevity, but all LAPACKE routines that return, return 0.

The same code with LAPACK_COL_MAJOR runs and valgrinds flawlessly.

My lapacke, lapack etc. is self-built for Ubuntu 12.04. I used the following settings in the lapack CMake file:

BUILD_COMPLEX       ON
BUILD_COMPLEX16     ON
BUILD_DOUBLE        ON
BUILD_SHARED_LIBS   ON
BUILD_SINGLE        ON
BUILD_STATIC_LIBS   ON
BUILD_TESTING       ON
CMAKE_BUILD_TYPE    Release
LAPACKE             ON
LAPACKE_WITH_TMG    ON

and the rest (the optimized blas/lapack and xblas) off. There were no errors during the build and all tests succeeded.

Where did I mess up?

Edit: I just tried this with Fedora21 and the packaged lapacke. It did not reproduce the error.

Edit 2: While it does not reproduce the memory fails, it produces a wrong solution, namely (1 + 0I, 1 + 0I) for the above input (should be (1,0))

Solution

  • After some more research and overthinking things, I found the culprit:

    Using LAPACK_ROW_MAJOR switches the meaning of the ld* leading dimension parameters. While the leading dimension of a normal Fortran array is the numbers of rows, switching to ROW_MAJOR switches its meaning to the number of columns. So the correct calls (giving correct results) would be:

    LAPACKE_zgetrs(LAPACK_ROW_MAJOR, 'N', 2, 1, mat, 2, p, vec, 1);

    where the second 2 is the number of columns (not rows!) of mat, and the last parameter must equal the number of right hand sides nrhs (not the number of variables!). I isolated this very call because all the other calls in my project dealt with square matrices, so the "wrong" calls do not have any negative effect due to symmetry.

    As usual, if you are skipping columns at the end, the leading dimensions get bigger accordingly, as they would with skipping rows in the normal setting.

    Obviously, this is not mentioned in the Fortran documentations. Unfortunately, I did see no such remark in the Lapacke documentation, which would have saved me a couple of hours of my life. :)