C++，LAPACK：使用链接库时出现“dsyev”的未定义引用错误

Question

我正在编写一个程序，需要使用LAPACK库子例程dsyev对矩阵进行对角化。首先，我使用了一个测试代码来确保我正确地链接了BLAS和LAPACK库。 测试代码如下：

// LAPACK test code
//compile with: g++ -Wall -L/lib64 -llapack -lblas LAPACK_testcode.cpp -o LAPACK_testcode.x

#include <iostream>
#include <vector>

using namespace std;

extern "C" void dgetrf_(int* dim1, int* dim2, double* a, int* lda, int* ipiv, int* info);
extern "C" void dgetrs_(char *TRANS, int *N, int *NRHS, double *A, int *LDA, int *IPIV, double *B, int *LDB, int *INFO );

int main()
{
char trans = 'N';
int dim = 2;
int nrhs = 1;
int LDA = dim;
int LDB = dim;
int info;

vector<double> a, b;

a.push_back(1);
a.push_back(1);
a.push_back(1);
a.push_back(-1);

b.push_back(2);
b.push_back(0);

int ipiv[3];

dgetrf_(&dim, &dim, &*a.begin(), &LDA, ipiv, &info);
dgetrs_(&trans, &dim, &nrhs, & *a.begin(), &LDA, ipiv, & *b.begin(), &LDB, &info);


std::cout << "solution is:";
std::cout << "[" << b[0] << ", " << b[1] << ", " << "]" << std::endl;
std::cout << "Info = " << info << std::endl;

return(0);

这段代码编译并运行没有问题，所以我不认为问题是与链接库有关。接下来我从Intel的网站上运行了一个dsyev示例C程序：

https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dsyev_ex.c.htm

而不是在这里发布代码，请查看超链接。 运行此程序（我将其保存为C ++程序，并确保LAPACK函数原型具有extern“C”），会产生以下错误：

dsyev_example.cpp: In function ‘int main()’:
dsyev_example.cpp:95:74: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
dsyev_example.cpp:95:74: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
dsyev_example.cpp:99:72: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
dsyev_example.cpp:99:72: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
dsyev_example.cpp:106:49: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
dsyev_example.cpp:108:72: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
/tmp/ccpoVl5n.o: In function `main':
dsyev_example.cpp:(.text+0x8b): undefined reference to `dsyev'
dsyev_example.cpp:(.text+0xf7): undefined reference to `dsyev'
collect2: error: ld returned 1 exit status

接着我尝试了我的程序，但是出现了类似的错误：

#include <cstdlib>
#include <iostream>
#include <iomanip>
#include <stdio.h>
#include "timer.h"

using namespace std;

const int N = 10;   // Matrix is of size N x N

// ******  NEED TO LINK BLAS AND LAPACK LIBRARIES IN ORDER TO COMPILE  ******
// g++ -Wall -L/lib64 -llapack -lblas matrix_multiply.cpp -o matrix_multiply.x

/* DSYEV prototype */
extern "C"
{
    void dsyev( char* jobz, char* uplo, int* n, double* a, int* lda,
        double* w, double* work, int* lwork, int* info );
}

/* DGEMM prototype */
extern "C"
{
    void dgemm(const char *transa, const char *transb, int *l, int *n, 
        int *m, double *alpha, const void *a, int *lda, void *b, 
        int *ldb, double *beta, void *c, int *ldc);
}

/* Auxiliary routines prototypes */
extern "C"
{
    void print_matrix( char* desc, int m, int n, double* a, int lda );
}

void initialize (double a [N][N]);
void print (int m, int n, double a [N][N]);

int main()
{
// declaration of matrix 
double A [N][N];            

// call initialize to initialize matrix here
initialize (A);

cout << "The initial matrix A is:\n";
// call print to print matrices here
print (N, N, A);

// Declaring Function Input for DSYEV
char jobz = 'V';    // 'V':  Compute eigenvalues and eigenvectors.
char uplo = 'U';    // 'U':  Upper triangle of A is stored.
int n = N;          // The order of the matrix A,  with n >= 0.
int lda = N;        // The leading dimension of the array A,  LDA >= max(1,N).
int lwork = N * N;  // The length of the array work, lwork >= max(1,3*N-1).
int info = 1;       // = 0:  successful exit
double w[N];        // w = array of eigenvalues in ascending order
double work[N];     // work = 

// Convert 2d array into 1d array for use in dsyev
double b [N*N];
for (int q = 0; q < n; q++)
{
    for (int t = 0; t < n; t++)
    {
        b[q * n + t] = A[q][t];
    }
}

// Print 1d array
for (int i = 0; i < N * N; i++)
{
    cout << setw (8) << i << " ";
}

dsyev(&jobz, &uplo, &n, b, &lda, w, work, &lwork, &info);

/* Check for convergence with dsyev */
if( info > 0 ) 
{
    cout << "The algorithm failed to compute eigenvalues." << endl;
    exit(1);
}

/* Print eigenvalues */
print_matrix( "Eigenvalues", 1, n, w, 1 );
/* Print eigenvectors */
print_matrix( "Eigenvectors (stored columnwise)", n, n, b, lda );

/*
timespec before, after;
get_time(&before);

// O(N^3) ALGORITHM FOR TWO-INDEX SIMILARITY TRANSFORMATION
for (int j = 0; j < N; j++)
{
    for (int p = 0; p < N; p++)
    {
        int sum = 0;

        for (int q = 0; q < N; q++)
        {
            sum += F[p][q] * T[q][j];
        }

        G[p][j] = sum;
    }
} 

for (int j = 0; j < N; j++)
{
    for (int i = 0; i < N; i++)
    {
        int sum = 0;

        for (int p = 0; p < N; p++)
        {
            sum += T[p][i] * G[p][j];
        }

        H[i][j] = sum;
    }
} 

get_time(&after);   
timespec time_diff;
diff(&before,&after,&time_diff);
// Time in seconds
double time_s = time_diff.tv_sec + (double)(time_diff.tv_nsec)/1.0e9;
cout << "\nThe time for the operation was " << time_s << endl;


// call print to print result here
print ();
*/

cout << "\n\n";

return 0;
}

// initialize sets all values in the matrix to desired form
void initialize (double a [N][N])
{    
for (int i = 0; i < N; i++)
{
    for (int j = 0; j < N; j++)
    {
        a [i][j] = 1.0;
        a [j][i] = 1.0;

        if (i < 5)
        {
            a [i][i] = 1.0 + (0.1 * i);
        }

        if (i > 5)
        {
            a [i][i] = 2 * (i + 1.0) - 1.0;
        }
    }
}
}

// prints the matrix
void print (int m, int n, double a [N][N])
{
for (int i = 0; i < m; i++)
{
    cout << endl;

    for (int j = 0; j < n; j++)
        cout << setw (8) << a [i][j];
}
cout << endl;
}

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, int m, int n, double* a, int lda ) 
{
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+j*lda] );
        printf( "\n" );
    }
}

我又一次遇到了以下错误：

matrix_multiply.cpp: In function ‘int main()’:
matrix_multiply.cpp:91:45: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
matrix_multiply.cpp:93:68: warning: deprecated conversion from string constant to ‘char*’ [-Wwrite-strings]
/tmp/ccVXf4OM.o: In function `main':
matrix_multiply.cpp:(.text+0x229): undefined reference to `dsyev'
collect2: error: ld returned 1 exit status

LAPACK和BLAS库都位于/lib64目录下（作为liblapack.so和libblas.so）。在两个发布的程序顶部都有被注释掉的命令。考虑到测试程序已经成功运行，我认为问题不是链接库？如果有人知道问题的原因，请指点一下。谢谢！

- Leigh K

Answer 1

注意extern函数声明的区别：

extern "C" void dgetrf_(int* dim1, ... ); // OK

extern "C"
{
    void dsyev( ... ); // KO
}

dsyev的声明缺少下划线_。这个下划线来自于将函数名从LAPACK源代码翻译成库文件对象中函数名的命名约定。例如，让我们看一下gfortran的命名约定：

默认情况下，过程名是小写Fortran名称，后面加一个下划线(_); 使用-fno-underscoring不会附加下划线，而使用-fsecond-underscore会附加两个下划线。取决于目标系统和调用约定，可能还需要对过程进行修饰;例如，在32位Windows上使用stdcall时，会附加一个@符号后跟一个整数号码。要更改调用约定，请参阅GNU Fortran编译器指令。

因此，命名约定可能取决于编译器、编译器选项和操作系统。因此，如果在C++代码中使用给定的命名约定声明函数，则该C++代码可能不可移植：在不同的计算机上或使用不同的编译器编译可能失败。

为了避免这种问题，存在一个名为LAPACKE的接口。这个接口负责命名约定。最后，你只需要在c++中通过extern "C" { #include <lapacke.h>}包含lapacke.h，链接到lapacke并在lapack函数名称前加上LAPACKE_。这里是一个使用dsyev()的简单示例C代码：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <lapacke.h>


int main(void){

    int n=200;

    srand(time(NULL));

    // allocate the matrix, unpacked storage
    double** A=malloc(n*sizeof(double*));
    if(A==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    A[0]=malloc(n*n*sizeof(double));
    if(A[0]==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    int i;
    for(i=1;i<n;i++){
        A[i]=&A[0][n*i];
    }

    //populte the matrix, only the lower diagonal part
    int j;
    for(i=0;i<n;i++){
        for(j=0;j<i;j++){
            A[i][j]=rand()/((double)RAND_MAX);

        }
        A[i][i]=rand()/((double)RAND_MAX)+42.;
    }

    //saving the matrix, because zheev() changes it.
    double complex** As=malloc(n*sizeof(double complex*));
    if(As==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    As[0]=malloc(n*n*sizeof(double complex));
    if(As[0]==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    for(i=1;i<n;i++){
        As[i]=&As[0][n*i];
    }
    for(i=0;i<n;i++){
        for(j=0;j<i;j++){
            As[i][j]=A[i][j];
        }
        As[i][i]=A[i][i];
    }
    //transpose part
    for(i=0;i<n;i++){
        for(j=i+1;j<n;j++){
            As[i][j]=(A[j][i]);
        }
    }

    //let's get the eigenvalues and the eigenvectors!
    double* w=malloc(n*sizeof(double));
    if(w==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    int lda = n;
    int info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'L', n, A[0], lda,  w);
    if(info<0){
        fprintf(stderr,"argument %d has illegal value\n",-info);
    }
    if(info>0){
        fprintf(stderr,"algorithm failed to converge... bad condition number ?\n");
    }

    //printing the eigenvalues...
    for(i=0;i<n;i++){
        printf("%d %g\n",i,w[i]);
    }

    //let's check that the column i of A is now a right eigenvectors corresponding to the eigenvalue w[i]...
    int l=42;

    double *res=malloc(n*sizeof(double));
    if(res==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
    for(i=0;i<n;i++){
        res[i]=0;
        for(j=0;j<n;j++){
            res[i]+=As[i][j]*A[j][l];
        }
        res[i]-=w[l]*A[i][l];
    }
    double norm2res=0;
    double norm2vec=0;
    for(i=0;i<n;i++){
        norm2res+=res[i]*res[i];
        norm2vec+=A[i][l]*A[i][l];
    }
    printf("the norm of the eigenvector is %g\n",sqrt(norm2vec));
    printf("||Ax-\\lambda*x||_2/||x||_2 is about %g\n",sqrt(norm2res/norm2vec));
    //free the matrix
    free(A[0]);
    free(A);

    free(As[0]);
    free(As);

    free(w);
    free(res);
    return 0;
}

这是由gcc main.c -o main -llapacke -llapack -lblas -lm -Wall编译的。

- francis

---

太棒了！非常感谢您提供如此详尽的回复！ - Leigh K