LAPACKE/GNU C++：LAPACKE_zheevx() 函数中的奇怪错误
LAPACKE/GNU C++: Weird bug in LAPACKE_zheevx() function
这里描述一下我在使用LAPACKE的函数LAPACKE_zheevx()时遇到的一个奇怪的bug。计算四分之三的简单测试代码 eigenvalues/vectors （来自英特尔网站的示例）运行良好并提供正确的输出。但是，如果我在源代码中引入任何字符串的声明（例如 std::string OutputFIlename;），编译进行得很顺利，但在执行时我得到分段错误 SIGSEGV ！！！
首先，我将列出有效的代码：
//===========================
#include <iostream>
#include <string>
#include <fstream>
#include <cassert>
#include <stdlib.h>
#include <stdio.h>
#include <complex>
#include <math.h>
#include "Headers_LAPACKE\lapacke.h"
#include "Headers_LAPACKE\lapacke_config.h"
#include "Headers_LAPACKE\lapacke_mangling.h"
#include "Headers_LAPACKE\lapacke_utils.h"
void print_matrix( char* desc, lapack_int m, lapack_int n, lapack_complex_double* a, lapack_int lda );
int main()
{
std::cout << "Start..." << std::endl;
//std::string fn_VALS;
// --------------- LAPACKE --- define variables ----------------------------------------------------------------
//
// Define arguments for LAPACKE_zheevx() routine:
int matrix_layout; // = LAPACK_ROW_MAJOR or LAPACK_COL_MAJOR
char jobz; // It is ="V" to calculate EigenVECTORS and
// ="N" if you don't do it.
char range; // ="A" for ALL values (don't want this),
// ="V" for values in in the half-open interval (vl,vu],
// ="I" for EigenVALUES indexed from il through iu (this is what you want).
char uplo; // ="U" for Upper Triangle of the matrix, or
// ="L" for the Lower triangle of the matrix.
lapack_int n; // the order of matrix "a" to be diagonalized.
lapack_complex_double* a; // complex array of dimension (lda,n).
// On entry it is Hermitian matrix "a".
// uplo="U" means the leading n-by-n upper triangular part of "a" contain the upper triangular part of the Hermitial matrix to be diagonalized.
// uplo="L" means equivalent for the lower triangular part.
// On exit, this content is destroyed.
lapack_int lda; // leading dimension of "a": lda >= max(1,n).
double vl; // taken into account only if range="V": vl<vu. Not referenced if range="I" or range="A".
double vu; // taken into account only if range="V": vl<vu. Not referenced if range="I" or range="A".
lapack_int il; // taken into account only if range="I": il<iu. Indices in ascending order of SMALLEST EigenVALUE to be returned. Not referenced if range="v" or range="A".
lapack_int iu; // taken into account only if range="I": il<iu. Indices in ascending order of LARGEST EigenVALUE to be returned. Not referenced if range="v" or range="A".
double abstol; // The absolute error tolerance for EigenVALUES. If abstop =<0, then EPS*|T| is used. If you get info>0 (some eigenvalues did not converge) then try abstol=2*DLAMCH('S').
lapack_int* m; // total number of EigenVALUES found: 0 =< m =< n, If range="A" then m=n, if range="I" then m = iu-il+1.
double* w; // double precision array of dimension n. On normal exit, the first m elements contain the selected EigenVALUES in ASCENDING ORDER.
lapack_complex_double* z; // double precision array of dimension (ldz, max(1,m)). If jobz="V" and info=0, the first m-columns of z contain normalized EigenVECTORS of a, corresponding to the selected EigenVALUES, with i-th column of z contains the Eigenvectro corresponding to w(i) eigenvalue.
lapack_int ldz; // leading dimension of array z. ldz>=1 and if jobz="V" then ldz >= max(1,n).
// following are used only with LAPACKE_zheevx_work() routine:.
//lapack_complex_double* work; // array of dimension max(1,lwork).
//lapack_int lwork; // lwork=-1 means workspace query. Othewise it has to be length of the array work: lwork=2*n for n>1, and lwork >=1 for N=<1.
//double* rwork; // array dimension is 7*n;
//lapack_int* iwork; // array dimension is 5*n;
lapack_int* ifail; // jobz="V" and info=0: first m elemens of ifail are zero. jobz="V" and info>0: ifail contain indices of the eigenvectors that failed to converge.
lapack_int info; //info=0 means successful exit. info>0 means eigenvectors failed to converge, their indices are in ifail. info<0, info=-i means i-th argument had an illegal value.
//
// ------------------------------------------------------------------------------------------------------------
matrix_layout = LAPACK_ROW_MAJOR;
jobz = 'V'; vl = 0.0; vu =100.0;
il = 1; iu=4; // are ignored now.
range = 'V';
uplo ='U';
n = 4;
lda = n;
ldz = n;
z = new lapack_complex_double [ldz*n];
w = new double [n];
a = new lapack_complex_double [lda*n];
a[0] =lapack_complex_double{6.51,0.0}; a[1] =lapack_make_complex_double(-5.92, 9.53); a[2]=lapack_complex_double{-2.46,2.91}; a[3]=lapack_complex_double{8.84,3.21};
a[4] =lapack_complex_double{0.0,0.0}; a[5] =lapack_make_complex_double(-1.73,0.0); a[6]=lapack_complex_double{6.5,2.09}; a[7]=lapack_complex_double{1.32, 8.81};
a[8] =lapack_complex_double{0.0,0.0}; a[9] =lapack_make_complex_double(0.0,0.0); a[10]=lapack_complex_double{6.90,0.0}; a[11]=lapack_complex_double{-0.59,2.47};
a[12]=lapack_complex_double{0.0,0.0}; a[13]=lapack_make_complex_double(0.0,0.0); a[14]=lapack_complex_double{0.0,0.0}; a[15]=lapack_complex_double{-2.85,0.0};
ifail = new lapack_int [n];
abstol = -1; // set default tolerance for calcuation of EigVals in the assigned interval.
print_matrix( "Entry Matrix A:", n, n, a, lda );
std::cout << std::endl;
info = LAPACKE_zheevx(matrix_layout, jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, ifail);
if (info>0)
{
std::cout << "Error: ZHEEVX failed to compute eigenvalues/vectors.";
exit(1);
}
std::cout << "info = " << info << std::endl;
std::cout << "Number of eigvals found: " << *m << std::endl;
for (int i_e =0; i_e<*m; i_e++)
{
std::cout << "Eigval. " << i_e << " is " << w[i_e] << std::endl;
}
print_matrix( "Selected EigVECTORS (column-wise):", n, n, z, ldz );
std::cout << std::endl;
std::cout << "Done :-) !!!" <<std::endl;
return 0;
}
////////////////////////////////////////////////////////* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, lapack_int m, lapack_int n, lapack_complex_double* a, lapack_int lda )
{
lapack_int i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ )
{
for( j = 0; j < n; j++ )
{
printf( " (%6.2f,%6.2f)", lapack_complex_double_real(a[i*lda+j]), lapack_complex_double_imag(a[i*lda+j]) );
}
printf( "\n" );
}
}
//=======================================
现在，如果在 main() 中删除第 2 行的注释符号 (//)： //std::string fn_VALS;这将变成 std::string fn_VALS；有了这个源代码将编译，但它在运行时失败并出现分段错误 SIGSEGV。
更多信息：
我正在使用 Windows 7 Pro 和 Code::Blocks，LAPACKE headers 和 dll 是在 6/15/2016 下载的。
从控制台：进程返回 -1073741819 (0xC0000005)
来自 Code::Blocks
中的调用堆栈 window
....... main() 调用 LAPACKE_zheevx() [lapacke.dll]
........ LAPACKE_zheevx() 调用 LAPACKE_zheevx_work() [lapacke.dll]
........ LAPACKE_zheevx_work() 调用 zheevx_() [lapack.dll]
请帮忙。
我尝试通过键入以下内容来编译您提供的程序：
g++ main.cpp -o main -llapacke -llapack -lblas -lm -Wall
标志-Wall启用所有警告。其中一个警告很有趣：
warning: ‘m’ is used uninitialized in this function [-Wuninitialized]
指针 lapack_int* m; 未在调用 LAPACKE_zheevx(..., m, ...); 时初始化。因此，它可以指向内存中的任何地方。指针 m 应该承载函数的输出参数，并且 LAPACKE_zheevx() 或后续函数很可能取消引用指针 m。它会导致未定义的行为：根据 m 指向的位置，它可能不被注意或触发分段错误。
你能试试 lapack_int* m[1]; 而不是 lapack_int* m; 吗？
This example from intel 不存在与 m 相同的问题，它被声明为 MKL_INT 并通过参数传递给函数 (&m)。是你的起点吗？
这里描述一下我在使用LAPACKE的函数LAPACKE_zheevx()时遇到的一个奇怪的bug。计算四分之三的简单测试代码 eigenvalues/vectors （来自英特尔网站的示例）运行良好并提供正确的输出。但是，如果我在源代码中引入任何字符串的声明（例如 std::string OutputFIlename;），编译进行得很顺利，但在执行时我得到分段错误 SIGSEGV ！！！
首先，我将列出有效的代码：
//===========================
#include <iostream>
#include <string>
#include <fstream>
#include <cassert>
#include <stdlib.h>
#include <stdio.h>
#include <complex>
#include <math.h>
#include "Headers_LAPACKE\lapacke.h"
#include "Headers_LAPACKE\lapacke_config.h"
#include "Headers_LAPACKE\lapacke_mangling.h"
#include "Headers_LAPACKE\lapacke_utils.h"
void print_matrix( char* desc, lapack_int m, lapack_int n, lapack_complex_double* a, lapack_int lda );
int main()
{
std::cout << "Start..." << std::endl;
//std::string fn_VALS;
// --------------- LAPACKE --- define variables ----------------------------------------------------------------
//
// Define arguments for LAPACKE_zheevx() routine:
int matrix_layout; // = LAPACK_ROW_MAJOR or LAPACK_COL_MAJOR
char jobz; // It is ="V" to calculate EigenVECTORS and
// ="N" if you don't do it.
char range; // ="A" for ALL values (don't want this),
// ="V" for values in in the half-open interval (vl,vu],
// ="I" for EigenVALUES indexed from il through iu (this is what you want).
char uplo; // ="U" for Upper Triangle of the matrix, or
// ="L" for the Lower triangle of the matrix.
lapack_int n; // the order of matrix "a" to be diagonalized.
lapack_complex_double* a; // complex array of dimension (lda,n).
// On entry it is Hermitian matrix "a".
// uplo="U" means the leading n-by-n upper triangular part of "a" contain the upper triangular part of the Hermitial matrix to be diagonalized.
// uplo="L" means equivalent for the lower triangular part.
// On exit, this content is destroyed.
lapack_int lda; // leading dimension of "a": lda >= max(1,n).
double vl; // taken into account only if range="V": vl<vu. Not referenced if range="I" or range="A".
double vu; // taken into account only if range="V": vl<vu. Not referenced if range="I" or range="A".
lapack_int il; // taken into account only if range="I": il<iu. Indices in ascending order of SMALLEST EigenVALUE to be returned. Not referenced if range="v" or range="A".
lapack_int iu; // taken into account only if range="I": il<iu. Indices in ascending order of LARGEST EigenVALUE to be returned. Not referenced if range="v" or range="A".
double abstol; // The absolute error tolerance for EigenVALUES. If abstop =<0, then EPS*|T| is used. If you get info>0 (some eigenvalues did not converge) then try abstol=2*DLAMCH('S').
lapack_int* m; // total number of EigenVALUES found: 0 =< m =< n, If range="A" then m=n, if range="I" then m = iu-il+1.
double* w; // double precision array of dimension n. On normal exit, the first m elements contain the selected EigenVALUES in ASCENDING ORDER.
lapack_complex_double* z; // double precision array of dimension (ldz, max(1,m)). If jobz="V" and info=0, the first m-columns of z contain normalized EigenVECTORS of a, corresponding to the selected EigenVALUES, with i-th column of z contains the Eigenvectro corresponding to w(i) eigenvalue.
lapack_int ldz; // leading dimension of array z. ldz>=1 and if jobz="V" then ldz >= max(1,n).
// following are used only with LAPACKE_zheevx_work() routine:.
//lapack_complex_double* work; // array of dimension max(1,lwork).
//lapack_int lwork; // lwork=-1 means workspace query. Othewise it has to be length of the array work: lwork=2*n for n>1, and lwork >=1 for N=<1.
//double* rwork; // array dimension is 7*n;
//lapack_int* iwork; // array dimension is 5*n;
lapack_int* ifail; // jobz="V" and info=0: first m elemens of ifail are zero. jobz="V" and info>0: ifail contain indices of the eigenvectors that failed to converge.
lapack_int info; //info=0 means successful exit. info>0 means eigenvectors failed to converge, their indices are in ifail. info<0, info=-i means i-th argument had an illegal value.
//
// ------------------------------------------------------------------------------------------------------------
matrix_layout = LAPACK_ROW_MAJOR;
jobz = 'V'; vl = 0.0; vu =100.0;
il = 1; iu=4; // are ignored now.
range = 'V';
uplo ='U';
n = 4;
lda = n;
ldz = n;
z = new lapack_complex_double [ldz*n];
w = new double [n];
a = new lapack_complex_double [lda*n];
a[0] =lapack_complex_double{6.51,0.0}; a[1] =lapack_make_complex_double(-5.92, 9.53); a[2]=lapack_complex_double{-2.46,2.91}; a[3]=lapack_complex_double{8.84,3.21};
a[4] =lapack_complex_double{0.0,0.0}; a[5] =lapack_make_complex_double(-1.73,0.0); a[6]=lapack_complex_double{6.5,2.09}; a[7]=lapack_complex_double{1.32, 8.81};
a[8] =lapack_complex_double{0.0,0.0}; a[9] =lapack_make_complex_double(0.0,0.0); a[10]=lapack_complex_double{6.90,0.0}; a[11]=lapack_complex_double{-0.59,2.47};
a[12]=lapack_complex_double{0.0,0.0}; a[13]=lapack_make_complex_double(0.0,0.0); a[14]=lapack_complex_double{0.0,0.0}; a[15]=lapack_complex_double{-2.85,0.0};
ifail = new lapack_int [n];
abstol = -1; // set default tolerance for calcuation of EigVals in the assigned interval.
print_matrix( "Entry Matrix A:", n, n, a, lda );
std::cout << std::endl;
info = LAPACKE_zheevx(matrix_layout, jobz, range, uplo, n, a, lda, vl, vu, il, iu, abstol, m, w, z, ldz, ifail);
if (info>0)
{
std::cout << "Error: ZHEEVX failed to compute eigenvalues/vectors.";
exit(1);
}
std::cout << "info = " << info << std::endl;
std::cout << "Number of eigvals found: " << *m << std::endl;
for (int i_e =0; i_e<*m; i_e++)
{
std::cout << "Eigval. " << i_e << " is " << w[i_e] << std::endl;
}
print_matrix( "Selected EigVECTORS (column-wise):", n, n, z, ldz );
std::cout << std::endl;
std::cout << "Done :-) !!!" <<std::endl;
return 0;
}
////////////////////////////////////////////////////////* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, lapack_int m, lapack_int n, lapack_complex_double* a, lapack_int lda )
{
lapack_int i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ )
{
for( j = 0; j < n; j++ )
{
printf( " (%6.2f,%6.2f)", lapack_complex_double_real(a[i*lda+j]), lapack_complex_double_imag(a[i*lda+j]) );
}
printf( "\n" );
}
}
//=======================================
现在，如果在 main() 中删除第 2 行的注释符号 (//)： //std::string fn_VALS;这将变成 std::string fn_VALS；有了这个源代码将编译，但它在运行时失败并出现分段错误 SIGSEGV。
更多信息：
我正在使用 Windows 7 Pro 和 Code::Blocks，LAPACKE headers 和 dll 是在 6/15/2016 下载的。 从控制台：进程返回 -1073741819 (0xC0000005) 来自 Code::Blocks
中的调用堆栈 window....... main() 调用 LAPACKE_zheevx() [lapacke.dll]
........ LAPACKE_zheevx() 调用 LAPACKE_zheevx_work() [lapacke.dll]
........ LAPACKE_zheevx_work() 调用 zheevx_() [lapack.dll]
请帮忙。
我尝试通过键入以下内容来编译您提供的程序：
g++ main.cpp -o main -llapacke -llapack -lblas -lm -Wall
标志-Wall启用所有警告。其中一个警告很有趣：
warning: ‘m’ is used uninitialized in this function [-Wuninitialized]
指针 lapack_int* m; 未在调用 LAPACKE_zheevx(..., m, ...); 时初始化。因此，它可以指向内存中的任何地方。指针 m 应该承载函数的输出参数，并且 LAPACKE_zheevx() 或后续函数很可能取消引用指针 m。它会导致未定义的行为：根据 m 指向的位置，它可能不被注意或触发分段错误。
你能试试 lapack_int* m[1]; 而不是 lapack_int* m; 吗？
This example from intel 不存在与 m 相同的问题，它被声明为 MKL_INT 并通过参数传递给函数 (&m)。是你的起点吗？