I am writing a program in which I need to diagonalize a matrix using the library subroutine dsyev in LAPACK. First, I used a test code to make sure that I was linking the BLAS and LAPACK libraries correctly. That code is:
// 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);
}
This code compiled and ran without issue, so I don't expect the problem to be an issue with linking the libraries. I next ran a dsyev example C program from the Intel website:
Instead of posting the code here, simply view the hyperlink. Running this program (I saved it as a C++ program and made sure the LAPACK function prototypes had extern "C") yields the following errors:
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
Then I tried my program, which had similar errors:
#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" );
}
}
And once again I get the following error:
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
Both the LAPACK and BLAS libraries are located in /lib64 (as liblapack.so and libblas.so). The command being used is commented in at the top of both posted programs. Considering the test program ran successfully, I would think that the issue is NOT linking the libraries? If anyone has any clues as to what might be the cause of the problem, I would very much appreciate the pointer. Thanks!
Notice the difference between the declaration of extern functions:
extern "C" void dgetrf_(int* dim1, ... ); // OK
extern "C"
{
void dsyev( ... ); // KO
}
The declaration of dsyev is missing an underscore _. Where does this underscore come from? This comes from the naming convention that translates the name of the function in the source code of LAPACK to the name of the function in the object file of the library. For instance, let's take a look at the 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.
As a result, the naming convention may depend on compiler, compiler options and operating system. Consequently, if the functions are decleared in the C++ code by using a given naming convention, the C++ code may not be portable: compiling it on a different computer or with a different compiler can fail.
To avoid this kind of issues, there exists an interface, named LAPACKE. This interface takes care of the naming convention. Finally, all you need to do is include lapacke.hby extern "C" { #include <lapacke.h>} (in c++), link to lapacke and prepend LAPACKE_ to the the name of the functions of lapack. Here is a simple sample C code using dsyev():
#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;
}
It is compiled by gcc main.c -o main -llapacke -llapack -lblas -lm -Wall