C++，拉帕克 : undefined reference to `dsyev' with libraries linked

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);

这段代码编译和运行没有问题，所以我不认为问题是链接库的问题。接下来，我从英特尔网站运行了一个 dsyev 示例 C 程序:

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

无需在此处发布代码，只需查看超链接即可。运行此程序(我将其保存为 C++ 程序并确保 LAPACK 函数原型(prototype)具有 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)。正在使用的命令在两个已发布程序的顶部都有注释。考虑到测试程序成功运行，我认为问题不在于链接库？如果有人对问题的可能原因有任何线索，我将非常感谢您的指点。谢谢!

Answer 1

注意外部函数声明的区别:

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

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

dsyev的声明缺少下划线 _ .这个下划线是从哪里来的？这来自将 LAPACK 源代码中的函数名称翻译为库的目标文件中的函数名称的命名约定。例如，让我们看一下 naming convention of gfortran :

By default, the procedure name is the lower-cased Fortran name with an appended underscore (_); using -fno-underscoring no underscore is appended while -fsecond-underscore appends two underscores. Depending on the target system and the calling convention, the procedure might be additionally dressed; for instance, on 32bit Windows with stdcall, an at-sign @ followed by an integer number is appended. For the changing the calling convention, see see GNU Fortran Compiler Directives.

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

为了避免此类问题，存在一个名为 LAPACKE 的接口(interface).该接口(interface)负责命名约定。最后，您需要做的就是包括 lapacke.h通过 extern "C" { #include <lapacke.h>} (在 C++ 中)，链接到 lapacke 并添加 LAPACKE_以lapack的功能名称。这是使用 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