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

}

此代碼編譯並沒有問題跑了，所以我不認爲這個問題是與鏈接庫的問題。我從旁邊英特爾網站跑了dsyev例如C程序：

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

在這裏張貼的代碼代替，只需查看超鏈接。運行這個程序（我保存它作爲一個C++程序，並確保該LAPACK函數原型有外部的「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函數聲明的區別：

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的界面。這個接口負責命名約定。最後，你需要做的只是包括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