Dear Piotr,
Thanks for the response.
In the meantime I have done some investigating myself in this problem, but did not find a conclusive solution yet.
Nonetheless I could confirm that the unnecessary static storage class definitions cause are one of the main problems to the usage of CLAPACK 3.0 in OpenMP. Two different threads running the same CLAPACK routine, could corrupt each others static variables.
I did not investigate further, but imagine that this could lead to the memory corruption I mentioned, there I think that the static variables are indirectly related to a bunch of pointer operations that could corrupt my memory if all of a sudden in the middle of a pointer operation, the needed variable changes.
The memory corruption I am speaking of is given by an error during execution of the program. (the program used routines are zpptrf, zpptri, zhptrf, zhptri) Each run gives a different error, all of them listed bellow.
FATAL: exception not rethrown
*** glibc detected *** ./MatrixElements-pbc.test.out: free(): invalid next size (fast): 0x000000000062c4a0 ***
*** glibc detected *** ./MatrixElements-pbc.test.out: free(): invalid pointer: 0x00002aaaaab06ff0 ***
terminate called without an active exception
*** glibc detected *** ./MatrixElements-pbc.test.out: malloc(): memory corruption: 0x0000000000608130 ***
complex division by zero
Segmentation fault
CLAPACK has been compiled with g++ in the standard way, grabbed straight from the site. (somehow i have troubles with compilation with icc (intel c compiler))
The problems occur in a program that does roughly this. In a loop running over the dimension of square hermitian complex matrices, matrices are filled with data and inverted with a function (i give a rough non working sketch of the function below). This loop is parallelised with openmp. This makes it possible that 2 function calls happen at the same time, and thus corrupt each others data trough the static variables.
This program is compiled with icpc (intel c++) with -openmp.
The program is currently solved a bit, because i have placed all clapack routines in critical omp pragma directives (only one can run at a time).
Before putting the critical directives, I found out that f2c contains a -a flag which changes the default storage class from static to automatic. Would this not be a problem?
I hope this is not to messy and is of any help to you and me both.
With kind regards,
Klaas Vantournhout
(rough sketch of a function)
- Code: Select all
MatrixCH Inverse(MatrixCH A) { // MatrixCH is complex hermitian matrix
lapack::doublecomplex *temp1;
temp1 = (lapack::doublecomplex *) malloc((n*(n+1))/2 * sizeof(lapack::doublecomplex)); // here n is the size of A
// PUT HERE DATA FROM A IN temp1 basically conversion from std::complex
// to lapack::doublecomplex
std2lapack(A,temp1);
long int N = n, info=0;
char UPLO = 'U';
lapack::zpptrf_(&UPLO, &N, temp1, &info);
if (info == 0) { // zpptrf
lapack::zpptri_(&UPLO, &N, temp1, &info);
if (info != 0) exit(1);
} else if (info > 0) { // zpptrf
// PUT DATA BACK
std2lapack(A,temp1);
long int *ipiv;
ipiv = (long int *) malloc (n * sizeof(long int)); // n size of A
lapack::zhptrf_(&UPLO, &N, temp1, ipiv, &info);
if (info == 0) {
lapack::doublecomplex *work;
work = (lapack::doublecomplex *) malloc (n * sizeof(lapack::doublecomplex));
lapack::zhptri_(&UPLO, &N, temp1, ipiv, work, &info);
if (info != 0) exit(1);
free(work); work = NULL; // C variable malloc <-> free
} else if (info != 0) exit(1);
free(ipiv); ipiv = NULL; // C variable malloc <-> free
} else if (info < 0) exit (1);
MatrixCH temp;
// put data from temp1 in temp
lapack2std(temp,temp1);
free(temp1); temp1 = NULL; // C variable malloc <-> free
temp.CleanUp();
return temp;
}