From: Laurent El Shafey <laurent.elshafey@id...>  20120419 20:00:07

Hello, I'm experiencing the following problem when interfacing Blitz++ arrays with a LAPACK function (dgesdd). To sum up, I'm getting a segmentation fault when using the Cpointers to the memory blocks of the 2D blitz::Array (which are supposed by construction to be contiguous in memory), whereas everything works fine when I'm directly using Cstyle arrays of the same size allocated with the new operator. The problem has been found while developing this library: http://idiap.github.com/bob/. The data containers chosen for this library are Blitz++ arrays. LAPACK functions require Cstyle arrays to be passed. Therefore, I'm passing the pointers to the memory blocks of the blitz::Array to the LAPACK function (using the Blitz++ data() method). Obviously, the blitz::Array should be contiguous in memory. Therefore, I've not set a SIMD width, and I'm not processing the arrays before passing them to the LAPACK function (In the attached example, arrays are created/allocated as in blitz::Array<double,2> A(3,3), for one of them initialized, and then the data() pointers are passed to LAPACK). As the code to reproduce the bug is a bit long, I've attached it and you can also find it at the end of this mail. Only 2D blitz::Array are leading to the problem. If you compile this example (g++ test.cc pthread llapack g o test), and run it, you will get the following segmentation fault: Program received signal SIGSEGV, Segmentation fault. 0xb726c19b in ATL_dgemvT_a1_x1_b0_y1 () from /usr/lib/libblas.so.3gf If you replace the pointers to the memory blocks of the 2D blitz::Array by memory blocks of the same size allocated with the new operator (Just uncomment the lines starting by "//UNCOMMENT" in the file attached), then everything is working fine. Strangely, the problem only occurs when using Ubuntu 32 bits (this has been checked on both Ubuntu 10.04 and 12.04 beta), and the latest checkout of Blitz++ from the mercurial repository. If Ubuntu 64 bits OR the latest (old) release Blitz++ 0.9 is used instead, everything is working fine, and there is no segmentation fault. This looks like a data alignment problem, but I don't manage to understand the problem completely, and hence, to find a nice fix. Thanks in advance for your help. With kind regards, Laurent p.s. 1: Of course, the dependence LAPACK needs to be installed (Ubuntu package liblapack3gf). p.s. 2: For blitz++, I've installed it through the following ppa which seems to be up to date: https://launchpad.net/~biometrics/+archive/blitz/+packages < SOURCE CODE TO REPRODUCE THE BUG > #include <blitz/array.h> #include <algorithm> #include <iostream> // Declaration of the external LAPACK function (Divide and conquer SVD) extern "C" void dgesdd_( const char *jobz, const int *M, const int *N, double *A, const int *lda, double *S, double *U, const int* ldu, double *VT, const int *ldvt, double *work, const int *lwork, int *iwork, int *info); int main() { const int N=3; blitz::Array<double,2> A(N,N); A = 0.8147, 0.9134, 0.2785, 0.9058, 0.6324, 0.5469, 0.1270, 0.0975, 0.9575; blitz::Array<double,1> S(N); blitz::Array<double,2> U(N,N); blitz::Array<double,2> Vt(N,N); // Initialises LAPACK variables const char jobz = 'A'; // Get All left singular vectors int info = 0; const int M = N; const int lda = N; const int ldu = N; const int ldvt = N; // Integer (workspace) array, dimension (8*min(M,N)) const int l_iwork = 8*N; int *iwork = new int[l_iwork]; // Initialises C arrays from blitz array (contiguous memory) double* A_lapack = A.data(); //UNCOMMENT A_lapack = new double[N*N]; //UNCOMMENT for(size_t i=0; i<N*N; ++i) A_lapack[i] = A(i); double *S_lapack = S.data(); double *U_lapack = Vt.data(); //UNCOMMENT U_lapack = new double[N*N]; double *VT_lapack = U.data(); //UNCOMMENT VT_lapack = new double[N*N]; // Call the LAPACK function: // A/ Query the optimal size of the working array const int lwork_query = 1; double work_query; dgesdd_( &jobz, &M, &N, A_lapack, &lda, S_lapack, U_lapack, &ldu, VT_lapack, &ldvt, &work_query, &lwork_query, iwork, &info ); // B/ Compute const int lwork = static_cast<int>(work_query); double *work = new double[lwork]; dgesdd_( &jobz, &N, &N, A_lapack, &lda, S_lapack, U_lapack, &ldu, VT_lapack, &ldvt, work, &lwork, iwork, &info ); // TODO: Check info variable std::cout << S << std::endl; // Free memory delete [] iwork; delete [] work; //UNCOMMENT delete [] A_lapack; //UNCOMMENT delete [] U_lapack; //UNCOMMENT delete [] VT_lapack; } 