Memory error in complex eigensolver (

[Carsten_F_]

SUMMARY:

===================

Deleting objects after successful calls to lapack_zhptrd, lapack_dstebz, lapack_zstein, and lapack_upmtr crashes with an error in free().

DESCRIPTION

============================

I'm writing an eigenvalue equation solver for complex hermitian matrices. Code is attached. It uses the sequence lapack_zhptrd, lapack_dstebz, lapack_zstein, and lapack_upmtr to 

1) generate the symmetric matrix T = Q^H A Q, with a unitary matrix Q. (lapack_zhptrd)

2.) Get m smallest eigenvalues for T (lapack_dstebz)

3.) Get m corresponding eigenvectors VT for T (lapack_zstein)

4.) Apply Q VT to get the m eigenvectors for A

Using the code below, it runs for dim(A) = N = 2000, M = 583 eigenvalues. increasing M to 584 results in a segfault. i tracked it down to the point where the temporary array tau gets cleaned up. this throws an exception in malloc.c as far as i understand.

Can anyone help?

CODE:

===================

file eigensolver.cpp

#include <iostream>
#include <cmath>
#include <complex>
#include "mkl_lapacke.h"
//
int eigensystem(MKL_Complex16 * A, double* lambda, MKL_Complex16* B, MKL_INT n, MKL_INT m);
//
int main() {

// Dimension
MKL_INT N = 2000;
// Number of solutions sought.
MKL_INT M = 584;

std::complex<double> * A = new std::complex<double>[N*(N+1)/2];
MKL_Complex16 * B = new MKL_Complex16[N*N];
double * lambda = new double;
// Dummy parameters.
double const onsite = -0.5;
double const re_offsite = 0.2;
double const im_offsite = 0.1;

for (int i = 0; i < N/2; ++i) {
A[N*i] = onsite;
A[N*i+1] = std::complex<double>(re_offsite, im_offsite);
A[N*i+3] = std::complex<double>(2.0 * re_offsite, -3.0*im_offsite);
}

eigensystem((MKL_Complex16*)A, lambda, B, N, M);

std::cout << "\n***************************************\n";
std::cout << "* Eigenvalues report *\n";
std::cout << "***************************************\n";
std::cout << 0 << "\t" << lambda[0] << std::endl;
std::cout << M-1 << "\t" << lambda[M-1] << std::endl;
//for (int i = 0; i < M; ++i) {
//std::cout << i << "\t" << lambda << std::endl;
//}

// Cleanup

delete [] A;
delete [] B;
delete [] lambda;

A = NULL;
B = NULL;
lambda = NULL;
}

int eigensystem(MKL_Complex16 * A, double * eValues, MKL_Complex16 * eVectors, MKL_INT n, MKL_INT const number_of_eigenvalues )
{

// Take a copy of A.
MKL_Complex16 * ap(A);

// Will contain error codes returned from MKL routines.
int info = 0;

// ----------------------------------------------------
// Convert matrix to symmetric tridiagonal form.
// T = Qt*A*Q.
//
// Setup working containers.
MKL_Complex16 * tau = new MKL_Complex16;
// Will contain the diagonal part of the symmetrized matrix.
double * d = new double;
// Will contain the off-diagonal part of the symmetrized matrix.
double * e = new double;

// Call LAPACK symmetrization routine.
info = LAPACKE_zhptrd( LAPACK_COL_MAJOR,
'U',
n,
ap,
d,
e,
tau
);

// Check results.

if (info != 0) {
std::cout << "\n*******************************\n ZHPTRD => info = " << info << "\n*******************************\n\n";
}
// ----------------------------------------------------
// Get eigenvalues.
//
// Input parameters.
char range = 'I'; // Select eigenvalues by indices.
char order = 'E'; // Order all eigenvalues in ascending order.
double const vl = 0.0; // Lower limit on selected eigenvalues.
double const vu = 0.0; // Upper limit on selected eigenvalues.
MKL_INT const il = 1; // Index of first selected eigenvalue.
MKL_INT const iu = number_of_eigenvalues; // Index of last selected eigenvalue.
double const abstol = 1e-8; // Absolute tolerance.

// Resize eigenvalues container.

// Output parameters.
MKL_INT * m = new MKL_INT; // Number of eigenvalues found.
MKL_INT * nsplit = new MKL_INT; // Number of diagonal blocks.
MKL_INT * iblock = new MKL_INT; // Block numbers of eigenvalues.
MKL_INT * isplit = new MKL_INT; // Split information.

// Call LAPACK eigenvalues routine.
info = LAPACKE_dstebz(
range,
order,
n,
vl,
vu,
il,
iu,
abstol,
(const double*) d,
(const double*) e,
m,
nsplit,
eValues,
iblock,
isplit);

// Check results.
if (info != 0) {
//throw MathException(LOCATION,"DSTEBZ failed.");
std::cout << "\n*******************************\n DSTEBZ => info = " << info << "\n*******************************\n\n";
}
if (*m != number_of_eigenvalues) {
//throw MathException(LOCATION,"Found fewer eigenvalues than sought.");
std::clog << "\nFound m = " << *m << " eigenvalues, sought " << number_of_eigenvalues << "." <<std::endl;
}

// -----------------------------------------------------
// Get eigenvectors of symmetric matrix T = Qt*A*Q.
//
// Input parameters.
MKL_INT * ifailv = new MKL_INT[*m];

// Since we did not call the eigenvalue routine stebz with range = 'B', we have to set all
// iblock = 1 and isplit[0] = n.
for (int i = 0; i < n; ++i) {
iblock = 1;
}
isplit[0] = n;

// Call LAPACK eigenvector routine.
info = LAPACKE_zstein(LAPACK_COL_MAJOR,
n,
(const double*) d,
(const double*) e,
*m,
(const double*) eValues,
(const MKL_INT*) iblock,
(const MKL_INT*) isplit,
eVectors,
n,
ifailv);

// Check results.
if (info != 0) {
std::cout << "\n*******************************\n ZSTEIN => info = " << info << "\n*******************************\n\n";
if (info > 0) {
for (int i = 0; i < info; ++i) {
std::cout << "Eigenvector #"<< ifailv << "failed to converge." << std::endl;
}
}
}

// ------------------------------------------------------------------------
// Back transformation V -> Q*V to get eigenvalues of original matrix.
// Multiply from left with Q.
//
info = LAPACKE_zupmtr(LAPACK_COL_MAJOR,
'L',
'U',
'N',
n, // !!! This must be the number of ROWS in the matrix!!!
*m, // !!! This must be the number of COLUMNS in the matrix!!!
(const MKL_Complex16*) ap,
(const MKL_Complex16*) tau,
eVectors,
n
);

// Check results.
if (info != 0) {
std::cout << "\n*******************************\n ZUPMTR => info = " << info << "\n*******************************\n\n";
}

std::cout<< "\n\nALL DONE\n" << std::endl;
// All done.
//
// Clean-up
std::clog << "Cleanup" << std::endl;

delete [] iblock;
delete [] isplit;
delete [] nsplit;
delete [] tau;
delete [] d;
delete [] e;
delete m;
delete [] ifailv;

iblock = NULL;
isplit = NULL;
nsplit = NULL;
tau = NULL;
d = NULL;
e = NULL;
m = NULL;
ifailv = NULL;

std::clog << "Return" << std::endl;
// Return.
return info;

}

SYSTEM

=========================

ubuntu 12.04 LTS, Kernel 3.2.0.55, gcc 4.6.4, intel mkl v 11.1 update 1, shipped with Intel Fortran Composer XE 2013 sp1.0.80

$> ulimit -a

core file size (blocks, -c) 0
data seg size (kbytes, -d) unlimited
scheduling priority (-e) 0
file size (blocks, -f) unlimited
pending signals (-i) 157796
max locked memory (kbytes, -l) 64
max memory size (kbytes, -m) unlimited
open files (-n) 1024
pipe size (512 bytes, -p) 8
POSIX message queues (bytes, -q) 819200
real-time priority (-r) 0
stack size (kbytes, -s) 65536
cpu time (seconds, -t) unlimited
max user processes (-u) 157796
virtual memory (kbytes, -v) unlimited
file locks (-x) unlimited

Compile and link line

======================

g++ -ggdb eigensolver.cpp -I/opt/intel/mkl/include -L/opt/intel/mkl/lib/intel64 -lmkl_intel_lp64 -lmkl_core -lmkl_intel_thread -lmkl_sequential -L/opt/intel/lib -liomp5 -lpthread -lm -o eigensolver

STDERR output

=====================

$> ./eigensolver

*** glibc detected *** ./eigensolver: munmap_chunk(): invalid pointer: 0x000000000083d260 ***
======= Backtrace: =========
/lib/x86_64-linux-gnu/libc.so.6(+0x7eb96)[0x7f8d38955b96]
./eigensolver[0x4013b8]
./eigensolver[0x400db3]
/lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xed)[0x7f8d388f876d]
./eigensolver[0x400b39]
======= Memory map: ========
00400000-00402000 r-xp 00000000 00:14 3420923 /home/cf/git/quantumsource/c++/eigensolver
00601000-00602000 r--p 00001000 00:14 3420923 /home/cf/git/quantumsource/c++/eigensolver
00602000-00603000 rw-p 00002000 00:14 3420923 /home/cf/git/quantumsource/c++/eigensolver
0083c000-008a5000 rw-p 00000000 00:00 0 [heap]
7f8d18000000-7f8d18090000 rw-p 00000000 00:00 0
7f8d18090000-7f8d1c000000 ---p 00000000 00:00 0
7f8d20000000-7f8d2007b000 rw-p 00000000 00:00 0
7f8d2007b000-7f8d24000000 ---p 00000000 00:00 0
7f8d24000000-7f8d2408c000 rw-p 00000000 00:00 0
7f8d2408c000-7f8d28000000 ---p 00000000 00:00 0
7f8d28000000-7f8d28021000 rw-p 00000000 00:00 0
7f8d28021000-7f8d2c000000 ---p 00000000 00:00 0
7f8d2ef5d000-7f8d2f377000 rw-p 00000000 00:00 0
7f8d2f377000-7f8d2f378000 ---p 00000000 00:00 0
7f8d2f378000-7f8d2f778000 rw-p 00000000 00:00 0
7f8d2f778000-7f8d2f779000 ---p 00000000 00:00 0
7f8d2f779000-7f8d2fb79000 rw-p 00000000 00:00 0
7f8d2fb79000-7f8d2fb7a000 ---p 00000000 00:00 0
7f8d2fb7a000-7f8d2ff7a000 rw-p 00000000 00:00 0
7f8d2ff7a000-7f8d319f9000 r-xp 00000000 08:05 2499559 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_avx.so
7f8d319f9000-7f8d31bf8000 ---p 01a7f000 08:05 2499559 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_avx.so
7f8d31bf8000-7f8d31c08000 r--p 01a7e000 08:05 2499559 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_avx.so
7f8d31c08000-7f8d31c10000 rw-p 01a8e000 08:05 2499559 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_avx.so
7f8d31c10000-7f8d37fa4000 rw-p 00000000 00:00 0
7f8d37fa4000-7f8d3809f000 r-xp 00000000 08:05 4460092 /lib/x86_64-linux-gnu/libm-2.15.so
7f8d3809f000-7f8d3829e000 ---p 000fb000 08:05 4460092 /lib/x86_64-linux-gnu/libm-2.15.so
7f8d3829e000-7f8d3829f000 r--p 000fa000 08:05 4460092 /lib/x86_64-linux-gnu/libm-2.15.so
7f8d3829f000-7f8d382a0000 rw-p 000fb000 08:05 4460092 /lib/x86_64-linux-gnu/libm-2.15.so
7f8d382a0000-7f8d382b8000 r-xp 00000000 08:05 4457015 /lib/x86_64-linux-gnu/libpthread-2.15.so
7f8d382b8000-7f8d384b7000 ---p 00018000 08:05 4457015 /lib/x86_64-linux-gnu/libpthread-2.15.so
7f8d384b7000-7f8d384b8000 r--p 00017000 08:05 4457015 /lib/x86_64-linux-gnu/libpthread-2.15.so
7f8d384b8000-7f8d384b9000 rw-p 00018000 08:05 4457015 /lib/x86_64-linux-gnu/libpthread-2.15.so
7f8d384b9000-7f8d384bd000 rw-p 00000000 00:00 0
7f8d384bd000-7f8d384d2000 r-xp 00000000 08:05 4460185 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f8d384d2000-7f8d386d1000 ---p 00015000 08:05 4460185 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f8d386d1000-7f8d386d2000 r--p 00014000 08:05 4460185 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f8d386d2000-7f8d386d3000 rw-p 00015000 08:05 4460185 /lib/x86_64-linux-gnu/libgcc_s.so.1
7f8d386d3000-7f8d386d5000 r-xp 00000000 08:05 4460099 /lib/x86_64-linux-gnu/libdl-2.15.so
7f8d386d5000-7f8d388d5000 ---p 00002000 08:05 4460099 /lib/x86_64-linux-gnu/libdl-2.15.so
7f8d388d5000-7f8d388d6000 r--p 00002000 08:05 4460099 /lib/x86_64-linux-gnu/libdl-2.15.so
7f8d388d6000-7f8d388d7000 rw-p 00003000 08:05 4460099 /lib/x86_64-linux-gnu/libdl-2.15.so
7f8d388d7000-7f8d38a8c000 r-xp 00000000 08:05 4456972 /lib/x86_64-linux-gnu/libc-2.15.so
7f8d38a8c000-7f8d38c8c000 ---p 001b5000 08:05 4456972 /lib/x86_64-linux-gnu/libc-2.15.so
7f8d38c8c000-7f8d38c90000 r--p 001b5000 08:05 4456972 /lib/x86_64-linux-gnu/libc-2.15.so
7f8d38c90000-7f8d38c92000 rw-p 001b9000 08:05 4456972 /lib/x86_64-linux-gnu/libc-2.15.so
7f8d38c92000-7f8d38c97000 rw-p 00000000 00:00 0
7f8d38c97000-7f8d38d79000 r-xp 00000000 08:05 270211 /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7f8d38d79000-7f8d38f78000 ---p 000e2000 08:05 270211 /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7f8d38f78000-7f8d38f80000 r--p 000e1000 08:05 270211 /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7f8d38f80000-7f8d38f82000 rw-p 000e9000 08:05 270211 /usr/lib/x86_64-linux-gnu/libstdc++.so.6.0.16
7f8d38f82000-7f8d38f97000 rw-p 00000000 00:00 0
7f8d38f97000-7f8d3907e000 r-xp 00000000 08:05 2495293 /opt/intel/composer_xe_2013_sp1.0.080/compiler/lib/intel64/libiomp5.so
7f8d3907e000-7f8d3927d000 ---p 000e7000 08:05 2495293 /opt/intel/composer_xe_2013_sp1.0.080/compiler/lib/intel64/libiomp5.so
7f8d3927d000-7f8d39288000 rw-p 000e6000 08:05 2495293 /opt/intel/composer_xe_2013_sp1.0.080/compiler/lib/intel64/libiomp5.so
7f8d39288000-7f8d392b0000 rw-p 00000000 00:00 0
7f8d392b0000-7f8d3a011000 r-xp 00000000 08:05 2499545 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_intel_thread.so
7f8d3a011000-7f8d3a211000 ---p 00d61000 08:05 2499545 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_intel_thread.so
7f8d3a211000-7f8d3a212000 r--p 00d61000 08:05 2499545 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_intel_thread.so
7f8d3a212000-7f8d3a347000 rw-p 00d62000 08:05 2499545 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_intel_thread.so
7f8d3a347000-7f8d3a34a000 rw-p 00000000 00:00 0
7f8d3a34a000-7f8d3b611000 r-xp 00000000 08:05 2499567 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_core.so
7f8d3b611000-7f8d3b811000 ---p 012c7000 08:05 2499567 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_core.so
7f8d3b811000-7f8d3b850000 r--p 012c7000 08:05 2499567 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_core.so
7f8d3b850000-7f8d3b85a000 rw-p 01306000 08:05 2499567 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_core.so
7f8d3b85a000-7f8d3b878000 rw-p 00000000 00:00 0
7f8d3b878000-7f8d3bd89000 r-xp 00000000 08:05 2499556 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_intel_lp64.so
7f8d3bd89000-7f8d3bf88000 ---p 00511000 08:05 2499556 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_intel_lp64.so
7f8d3bf88000-7f8d3bf8a000 r--p 00510000 08:05 2499556 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_intel_lp64.so
7f8d3bf8a000-7f8d3bf96000 rw-p 00512000 08:05 2499556 /opt/intel/composer_xe_2013_sp1.0.080/mkl/lib/intel64/libmkl_intel_lp64.so
7f8d3bf96000-7f8d3bf9b000 rw-p 00000000 00:00 0
7f8d3bf9b000-7f8d3bfbd000 r-xp 00000000 08:05 4460093 /lib/x86_64-linux-gnu/ld-2.15.so
7f8d3c014000-7f8d3c172000 rw-p 00000000 00:00 0
7f8d3c172000-7f8d3c173000 ---p 00000000 00:00 0
7f8d3c173000-7f8d3c18a000 rw-p 00000000 00:00 0
7f8d3c198000-7f8d3c1bd000 rw-p 00000000 00:00 0
7f8d3c1bd000-7f8d3c1be000 r--p 00022000 08:05 4460093 /lib/x86_64-linux-gnu/ld-2.15.so
7f8d3c1be000-7f8d3c1c0000 rw-p 00023000 08:05 4460093 /lib/x86_64-linux-gnu/ld-2.15.so
7ffff40c6000-7ffff40e8000 rw-p 00000000 00:00 0 [stack]
7ffff41ff000-7ffff4200000 r-xp 00000000 00:00 0 [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0 [vsyscall]
Aborted (core dumped)

GDB backtrace

===========================

(gdb) bt
#0 0x00007ffff474c425 in __GI_raise (sig=<optimized out>)
at ../nptl/sysdeps/unix/sysv/linux/raise.c:64
#1 0x00007ffff474fb8b in __GI_abort () at abort.c:91
#2 0x00007ffff478a39e in __libc_message (do_abort=2,
fmt=0x7ffff4894748 "*** glibc detected *** %s: %s: 0x%s ***\n")
at ../sysdeps/unix/sysv/linux/libc_fatal.c:201
#3 0x00007ffff4794b96 in malloc_printerr (action=3,
str=0x7ffff4894838 "munmap_chunk(): invalid pointer", ptr=<optimized out>)
at malloc.c:5039
#4 0x00000000004013b8 in eigensystem (A=0x7ffff1f5a010, eValues=0x603010,
eVectors=0x7fffee250010, n=2000, number_of_eigenvalues=584)
at eigensolver.cpp:202
#5 0x0000000000400db3 in main () at eigensolver.cpp:29

[mecej4]

It may be more profitable for you to note that, for N = 2000, M = 583 (in fact, even for M=50 with the same N) the program gives

[bash] ZSTEIN => info = -6[/bash]

When such an error occurs, you should probably terminate the program instead of continuing the calculations and calling ZUMPTR with invalid results from the earlier calculation. 

In other words, look into the failure of ZSTEIN instead of the failure of ZUMPTR.

[Carsten_F_]

Hi mecej4.

Thank for your reply.

In other words, look into the failure of ZSTEIN instead of the failure of ZUMPTR.

That's even more weird,  because I don't get that message. Here's the complete log from my run (except the backtrace).

$> ./eigensolver

ALL DONE

Cleanup
*** glibc detected *** ./eigensolver: munmap_chunk(): invalid pointer: 0x0000000002529260 ***
======= Backtrace: =========

[Gennady_F_Intel]

pls, set                  double * lambda = new double;

instead of             double * lambda = new double;

see from documentation for  dstebz 

w

REAL for sstebz

DOUBLE PRECISION for dstebz.

Array, DIMENSION at least max(1, n). The computed eigenvalues, stored in w(1) to w(m).

here is output of your example after this correction:

ALL DONE

Cleanup
Return

***************************************
*          Eigenvalues report         *
***************************************
0       -2.70582
583     0
Press any key to continue . . .

regards, Gennady