Segmentation fault when using lapacke_dsyev_row with more than 1300x1300 elements on 64-bit system

[arthur_dent]

Hello,

I have installed MKL package on 64-bit linux and try to calculate eigenvalues of a matrix consisting of more than 1300x1300 elements (actually I need 12kx12k elements). When I compile with the latest version of icc (make libintel64 function=lapacke_dsyev_row) Im always getting segmentation fault if N and LDA number exceed or are equal to 1300.

here's the output:
:/home/myuser/intel/composerxe-2011.5.220/mkl/examples/lapacke> uname -a
Linux chassagne.simap.grenoble-inp.fr 2.6.18-238.19.1.el5 #1 SMP Fri Jul 15 07:31:24 EDT 2011 x86_64 x86_64 x86_64 GNU/Linux
Chassagne:/home/myuser/intel/composerxe-2011.5.220/mkl/examples/lapacke> make libintel64 function=lapacke_dsyev_row

----- Compiling intel_lp64_parallel_intel64_lib ----- lapacke_dsyev_row
icc -03 -shared-intel -mcmodel=large -w -I"/home/myuser/intel/composer_xe_2011_sp1.6.233/mkl/include" \\
./source/lapacke_dsyev_row.c -Wl,--start-group \\
"/home/myuser/intel/composer_xe_2011_sp1.6.233/mkl/lib/intel64/libmkl_intel_lp64.a" \\
"/home/myuser/intel/composer_xe_2011_sp1.6.233/mkl/lib/intel64/libmkl_intel_thread.a" \\
"/home/myuser/intel/composer_xe_2011_sp1.6.233/mkl/lib/intel64/libmkl_core.a" \\
-Wl,--end-group -L"/home/myuser/intel/composer_xe_2011_sp1.6.233/mkl/../compiler/lib/intel64" -liomp5 -lpthread -lm -o _results/intel_lp64_parallel_intel64_lib/lapacke_dsyev_row.out
----- Execution intel_lp64_parallel_intel64_lib ----- lapacke_dsyev_row
export LD_LIBRARY_PATH="/home/myuser/intel/composer_xe_2011_sp1.6.233/mkl/lib/intel64":/home/myuser/intel/composer_xe_2011_sp1.6.233/compiler/lib/intel64:/home/myuser/intel/composer_xe_2011_sp1.6.233/ipp/../compiler/lib/intel64:/home/myuser/intel/composer_xe_2011_sp1.6.233/ipp/lib/intel64:/home/myuser/intel/composer_xe_2011_sp1.6.233/compiler/lib/intel64:/home/myuser/intel/composer_xe_2011_sp1.6.233/mkl/lib/intel64:/home/myuser/intel/composer_xe_2011_sp1.6.233/tbb/lib/intel64//cc4.1.0_libc2.4_kernel2.6.16.21:/opt/intel/fce/10.0.025/lib:/home/myuser/intel/composer_xe_2011_sp1.6.233/mpirt/lib/intel64:/home/myuser/intel/composer_xe_2011_sp1.6.233/mkl/../compiler/lib/intel64; \\
_results/intel_lp64_parallel_intel64_lib/lapacke_dsyev_row.out > _results/intel_lp64_parallel_intel64_lib/lapacke_dsyev_row.res
/bin/sh: line 1: 28206 Erreur de segmentation _results/intel_lp64_parallel_intel64_lib/lapacke_dsyev_row.out > _results/intel_lp64_parallel_intel64_lib/lapacke_dsyev_row.res
make[1]: *** [lapacke_dsyev_row] Erreur 139
make: *** [libintel64] Erreur 2


could you please tell me what is going on?
thanks in advance

[barragan_villanueva_]

Original test works and gives the following results

LAPACKE_dsyev (row-major, high-level) Example Program Results

Eigenvalues
-11.07 -6.23 0.86 8.87 16.09

Eigenvectors (stored columnwise)
-0.30 -0.61 0.40 -0.37 0.49
-0.51 -0.29 -0.41 -0.36 -0.61
-0.08 -0.38 -0.66 0.50 0.40
-0.00 -0.45 0.46 0.62 -0.46
-0.80 0.45 0.17 0.31 0.16

What changes were done in the test by you?

[arthur_dent]

I just changed N value to 1300, the matrix was completed with zero values.

in other example of the same function I read values from a file but it only works for N < 1300.

[Gennady_F_Intel]

how do you allacate the working array? it might caused the segmentation if you just allocate like

N = 1300

double a[N*N] .

--Gennady

[arthur_dent]

yes, it is exactly like in the lapacke example c file. I just changed #define N 5 to #define N 1300. The array is then double a[N*LDA], where LDA = N. How should I allocate it?

thanks

[Todd_R_Intel]

Well it shouldn't be exactly like the definition of the example, because the example allocates it by defining it explicitly with some data provided. You would need allocate enough space.

If you attach your source we may be able to provide you more advice.

Todd

[Gennady_F_Intel · Accepted solution]

it would help you if you dynamically allocate the working array:

double* a = (double*)malloc( LDA*N*sizeof(double));

[arthur_dent]

Here's my .c file. The DYNAMICAL_MATRIX.dat file consists of one column with 2250000 lines (matrix elements for 500-atom system). thank you for your help.

#include
#include
#include "mkl_lapacke.h"

/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda );

/* Parameters */
#define N 1500
#define LDA N

/* Main program */
int main() {
/* Locals */
int i = 0;
MKL_INT n = N, lda = LDA, info;
FILE *fp= fopen("DYNAMICAL_MATRIX.dat", "r");

/* Local arrays */
double w;
double a[LDA*N];

if(!fp)
{
puts("File not found.");
return EXIT_FAILURE;
}
while (feof(fp) == 0)
{
fscanf(fp, "%lf\n", &a);
++i;
}
fclose(fp);

/* Executable statements */
printf( "LAPACKE_dsyev (row-major, high-level) Example Program Results\n" );
/* Solve eigenproblem */
info = LAPACKE_dsyev( LAPACK_ROW_MAJOR, 'V', 'U', n, a, lda, w );
/* Check for convergence */
if( info > 0 ) {
printf( "The algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
/* Print eigenvalues */
print_matrix( "Eigenvalues", 1, n, w, 1 );
/* Print eigenvectors */
print_matrix( "Eigenvectors (stored columnwise)", n, n, a, lda );
exit( 0 );
} /* End of LAPACKE_dsyev Example */

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda ) {
MKL_INT i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %6.2f", a[i*lda+j] );
printf( "\n" );
}
}

[mecej4]

The local array a[] would attempt to consume 9 MB on the stack. If your default stack limit is, say, 8 MB, or if your data file contains more lines of data than N*N, your program would crash.

1. As Gennady suggested, allocate the array on the heap, instead, using malloc() or calloc().

2. Make sure that the number of items read from the file is exactly N², for example, by checking that the loop counter i has this value when end-of-file occurs.

[arthur_dent]

the number of lines is 2250000 when I do wc -l on this file so it seems like there is no empty line. I will allocate the array on the heap. thank you. I will let you know if it worked. Thank you guys.

[arthur_dent]

It worked. Thank you!