Reputation: 50081
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)
)
Upvotes: 1
Views: 449
Reputation: 50081
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. :)
Upvotes: 0