Linking: undefined reference to `sgees_'

Post here if you have a question about linking your program with LAPACK or ScaLAPACK library

Linking: undefined reference to `sgees_'

Postby epolyach » Tue Jun 11, 2013 9:01 am

Hi,

I am rewriting MATLAB code to C. The problem I've encountered is to calculate eigenvectors and eigenvalues of large (1000x1000) matrix. The MATLAB function is "eig" actully uses LAPACK, so I am trying to reproduce <eig> results. I can reproduce <eig> results using GSL, but it is slow.

According to the rules, the function I need should be "LAPACKE_sgees". The compiler finds the library /usr/lib/liblapacke.a, but complains as follows:

Code:
// lapack2
LAPACKE_sgees(LAPACK_ROW_MAJOR, 'V', 'N', select, NRB, A_lapack, NRB, &sdim, wr, wi, vs, NRB );


martin:~/Desktop/WORK/radial/Stability/C_files> make
g++ aniz_field.o my_useful_f.o -o aniz_field.exe -L/usr/lib -lgsl -lgslcblas -lm -llapack -llapacke -lblas -gfortran
/usr/lib/liblapacke.a(lapacke_sgees_work.o): In function `LAPACKE_sgees_work':
lapacke_sgees_work.c:(.text+0x1ee): undefined reference to `sgees_'
lapacke_sgees_work.c:(.text+0x320): undefined reference to `sgees_'
collect2: ld returned 1 exit status
make: *** [aniz_field.exe] Error 1

I have succeeded to get occasionally some results using a "sgees_", although I found no explanation why this is the function I need.

Code:
// lapack1
// sgees_(&c5, &c4, select, &c1, A_lapack, &c1, &sdim, wr, wi, vs, &c1, sWORK, &c2, bWORK, &info);

However, the result differs from the results obtained with MATLAB and GSL, see figure.

evec.png
evec.png (5.28 KiB) Viewed 8846 times


Could you tell me how to overcome linking errors and what is the reason for the discrepancy? I give the Makefile and the code below.

Here is the Makefile
Code: Select all
 CC=g++
LAPACK_PATH = /usr/lib/lapack
LIB= -L/usr/lib -lgsl -lgslcblas -lm -llapack -llapacke  -lblas -gfortran
INCLUDE= -I/usr/include/lapack
aniz_field.exe: aniz_field.o my_useful_f.o
   $(CC) $^ -o $@ $(LIB) 

.c.o:
   $(CC) $(INCLUDE) -c -g $^


Here is the code
Code: Select all
/*****************************************************************************
  File Name      : "aniz_field.c"
  Contents       : PP81 matrix method for calculation of eigenmodes
                      :
  Last version   : 2013.VI.04. 12:59
  Version No.    : 1.0
*****************************************************************************/

#define   FILE_INI   "aniz_field.cfg"

#define NA        0
#define NB        101         
#define NE        (NA+NB)
#define L1MAX     15
#define NALPHA    50

#define NPOT     502
#define BASIS     3
#define NRB        1001

#define NR        10001
#define NTH       3001
#define   TH_EPS    .5e-5
#define R_EPS     1e-14


/***************************************************************/

#include <stdio.h>
#include <iostream> 
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/stat.h>

#include <lapacke.h>
#include <lapacke_utils.h>
#include <lapacke_mangling.h>
#include <lapacke_config.h>

#include <gsl/gsl_math.h>
//#include <gsl/gsl_errno.h>
//#include <gsl/gsl_roots.h>
#include <gsl/gsl_complex_math.h>
//#include <gsl/gsl_linalg.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
//#include <gsl/gsl_spline.h>
#include <gsl/gsl_eigen.h>
//#include <gsl/gsl_blas.h>

#include "my_useful_f.h"

/***************************************************************/
/* Some "good" functions... */

#define SIG(x)   ( ( (x)< 0 ) ? (-1):(1) )
#define ABS(x)   ( ( (x)< 0 ) ? (-x):(x) )

#define MAX(a,b) ( ( (a)>(b) ) ? (a):(b) )
#define MIN(a,b) ( ( (a)<(b) ) ? (a):(b) )

#define SQR(x)   (    (x) * (x) )
#define POW3(x)  ( SQR(x) * (x) )

#define Pi     3.14159265358979323846

/***************************************************************/

//  global vars

// gsl_function        F;
// gsl_spline          *spline[NALPHA];

float               tmp;
FILE              *inp, *out;
int             q, l,
               Nwx, Nwy;
               
float            L_T, chi,
               w_x_min, w_x_max,
               w_y_min, w_y_max,
               N;

float              rp[NPOT], Psi[NPOT], dPsi[NPOT],
               rho[NPOT], drho[NPOT];

float            Fcal[NRB][NRB],
               r_chi[NRB], Sr_chi[NRB], g_chi[NRB];
               
float            **chi_b, *lambda_b;               

char             dirname[10], BUFF[1024];


/***************************************************************/

void fsave1s(char * name, float *array, int Nx){

   FILE *fp1;
   if((fp1=fopen(my_strcat(dirname, name), "w"))==NULL) {
      printf("Error opening output file. \n");
   }
   for(int i=0; i<Nx; i++)
      fprintf(fp1, "%.7E\n", array[i]);      
   fclose(fp1);
}
/***************************************************************/

void fsave2s(char * name, float **array, int Nx, int Ny){

   FILE *fp1;
   // "wb+"; fwrite(&array[i][j], sizeof(float), 1, fp1);
   if((fp1=fopen(my_strcat(dirname, name), "w"))==NULL) {
      printf("Error opening output file. \n");
   }
   for(int i=0; i<Nx; i++)
      for(int j=0; j<Ny; j++)
         fprintf(fp1, "%.7E\n", array[i][j]);      
   fclose(fp1);
}
/***************************************************************/

void fsave2s1(char * name, float *array, int Nx, int Ny){

   FILE *fp1;
   if((fp1=fopen(my_strcat(dirname, name), "w"))==NULL) {
      printf("Error opening output file. \n");
   }
   for(int i=0; i<Nx; i++)
      for(int j=0; j<Ny; j++)
         fprintf(fp1, "%.7E\n", array[i*Ny + j]);      
   fclose(fp1);
}

/***************************************************************/

void init_chi(void){

   double    dr;
   int i,j;
   
   chi_b = new float *[NRB];
   for(i = 0; i <NALPHA; i++)
       chi_b[i] = new float[NRB];
    lambda_b = new float [NALPHA];
   
   dr = 1./(NRB-1);
   for(i=0;i<NRB;i++)
      r_chi[i] = dr*i;
   for(i=2;i<NRB;i+=2)
      Sr_chi[i] = dr*2./3.;
   for(i=1;i<NRB-1;i+=2)
      Sr_chi[i] = dr*4./3.;
   Sr_chi[0] = dr;
   Sr_chi[NRB-1] = dr;
   
   for(i=0;i<NRB;i++)
      for(j=0;j<NRB;j++)
         Fcal[i][j] = pow(MIN(r_chi[i], r_chi[j]),l)/pow(MAX(r_chi[i], r_chi[j]),l+1);
   Fcal[0][0] = 1.;
}

/***************************************************************/

void basis_exp(void){

   int i;
   
   for(i=0;i<NRB;i++)
      g_chi[i] = -exp(-r_chi[i]/L_T);

}

/***************************************************************/

void read_cfg(void){

    inp = fopen(FILE_INI,"r");
    if (inp==NULL)
    {
        printf("Current directory must have %s file. Stop.\n",  FILE_INI);
        exit;
    }
   
    fgets(BUFF, 1024, inp);
    sscanf(BUFF,"%d %d %E %E", &l, &q, &L_T, &chi);
    fgets(BUFF, 1024, inp);
    sscanf(BUFF,"%d %d", &Nwx, &Nwy);
    fgets(BUFF, 1024, inp);
    sscanf(BUFF,"%E %E", &w_x_min, &w_x_max);
    fgets(BUFF, 1024, inp);
    sscanf(BUFF,"%E %E", &w_y_min, &w_y_max);
    fgets(BUFF, 1024, inp);
    sscanf(BUFF,"%E", &N);
   
    printf("Model q = %d, LT = % .4E , harmonics l = %d \n", q, L_T, l);     
   
    for(int i=0;i<NPOT;i++)
    {
        fgets(BUFF, 1024, inp);
        sscanf(BUFF,"%E %E %E %E %E\n", &rp[i], &Psi[i], &dPsi[i], &rho[i], &drho[i]);
        // printf(" r = % .4E, Psi = % .4E , dPsi = % .4E \n", rp[i], Psi[i], dPsi[i]);
    }
    fclose(inp);
}

/***************************************************************/

void calc_basis(void){

   gsl_matrix         *A =  gsl_matrix_alloc(NRB,NRB);    
   gsl_vector_complex *eval = gsl_vector_complex_alloc (NRB);
   gsl_matrix_complex *evec = gsl_matrix_complex_alloc (NRB, NRB);
   gsl_eigen_nonsymmv_workspace *workspace = gsl_eigen_nonsymmv_alloc(NRB);
    gsl_complex z, eval_i;
   int i, j;
   
   float       A_lapack[NRB*NRB];

   double aij;
   
   for(int ia = 0; ia < NRB; ia++)
      for(int ja = 0; ja < NRB; ja++){
         aij = -Fcal[ia][ja] * SQR(r_chi[ja]) * g_chi[ja] *Sr_chi[ja]/(2*l+1);
         A_lapack[ia+ja*NRB] = aij;
         gsl_matrix_set(A,ia,ja,aij);
      }
   // printf("hello 1.\n");
    gsl_eigen_nonsymmv(A, eval, evec, workspace);
   // printf("hello 2.\n");
    gsl_matrix_free(A);   
    gsl_eigen_nonsymmv_free (workspace);     
    gsl_eigen_nonsymmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC);
       
   for (i = 0; i < NALPHA; i++)
   {
      eval_i = gsl_vector_complex_get (eval, i);
      lambda_b[i] = 1./GSL_REAL(eval_i);
      // printf("eval (%d) %E \n", i, lambda_b[i]);
      gsl_vector_complex_view evec_i = gsl_matrix_complex_column (evec, i);
      for (j = 0; j < NRB; j++){
         z = gsl_vector_complex_get(&evec_i.vector, j);
         chi_b[i][j] = GSL_REAL(z);
      }
   }   
   gsl_vector_complex_free(eval);
   gsl_matrix_complex_free(evec);

   fsave1s("lambda_b", lambda_b, NALPHA);   
   fsave2s("chi_b", chi_b, NALPHA, NRB);   
   
   float *wr = new float [NRB];
   float *wi = new float [NRB];
   float *vs = new float [NRB*NRB];   
   float *sWORK = new float [3*NRB];
   int *bWORK;
   LAPACK_S_SELECT2 select;
   
   char c4,c5;
   int c1,c2,c3,info,sdim;
   c1 = NRB;
   c2 = NRB*3;
   
   c4 = 'N';
   c5 = 'V';

   // lapack1
   // sgees_(&c5, &c4, select, &c1, A_lapack, &c1, &sdim, wr, wi, vs, &c1, sWORK, &c2, bWORK, &info);

   // lapack2
   LAPACKE_sgees(LAPACK_ROW_MAJOR, 'V', 'N', select, NRB, A_lapack, NRB, &sdim, wr, wi, vs, NRB );
   

   if (!info)
      printf("lapack finished ok.\n");
   else
      printf("lapack has problems (%d).\n", info);
    // zgeev_(&c4, &c4,&c1, AT, &c1, b, DUMMY, &c3, DUMMY, &c3, sWORK, &c2, bWORK, &ok);     

   for(i=0;i<NALPHA;i++)
      wi[i] = 1./wr[i];
      
   fsave1s("lambda_lapack_r", wi, NALPHA);   
   //fsave1s("lambda_lapack_i", wi, NALPHA);   
   fsave2s1("chi_lapack", vs, NRB, NRB);   
   
}

/***************************************************************/

LAPACK_S_SELECT2 select(float x, float y)
{
   LAPACK_S_SELECT2 a;   
   return *a;      
}

/***************************************************************/

int main(){

//  local vars
    int             i, status;
       
/***** read cfg file with equilibrium data *********************/

   read_cfg();
   status=sprintf(dirname, "%.0f-%d/", L_T*100000, l); 
   printf("dirname = %s\n", dirname);
   status = mkdir(dirname,0755);
   if (status)
   {
     printf("Unable to create directory\n");
     exit(1);
   }
/***** choosing basis chi[][] and prepare splines **************/

   init_chi();
   
   switch (l)
   {
      case 0:
          puts("Not implemented. Stop.");
           return(0);
          break;
         
      default:
         switch (BASIS)
         {
            case 0: // Bessel
               ;
               break;
            
            case 1: // rho
               ;
               break;
            
            case 2: // rho'/Phi'
               ;
               break;
            
            case 3: // custom - exp(-r/L_T)
               basis_exp();
               printf("Basis functions: custom (4)...");
               break;
            
            default:
               ;
         }
   }
   
   calc_basis();
   printf(" calculated.\n");
   

/***** circular orbits, phase space ****************************/



/***** frequencies Om1, Om2, vpha, Sab *************************/



/***** matrix M(\alpha, \beta, w_x, w_y) ***********************/



/***** det(M) **************************************************/



/***** free space **********************************************/



}


epolyach
 
Posts: 2
Joined: Tue Jun 11, 2013 8:20 am

Re: Linking: undefined reference to `sgees_'

Postby admin » Tue Jun 11, 2013 3:50 pm

Try to put the libraries in order at linking time:
Code: Select all
 -lm -llapacke -llapack -lblas -gfortran


Hope it helps
admin
Site Admin
 
Posts: 468
Joined: Wed Dec 08, 2004 7:07 pm

Re: Linking: undefined reference to `sgees_'

Postby epolyach » Tue Jun 11, 2013 7:22 pm

Thanks a lot for the advise! I have linked safely. Why it is?

However, the results with

Code: Select all
LAPACKE_sgees(LAPACK_ROW_MAJOR, 'V', 'N', select, NRB, A_lapack, NRB, &sdim, wr, wi, vs, NRB );



are of the same accuracy as with

Code: Select all
sgees_(&c5, &c4, select, &c1, A_lapack, &c1, &sdim, wr, wi, vs, &c1, sWORK, &c2, bWORK, &info);


I have also checked whether it is because of using single precision (float) instead of double. It is not.
epolyach
 
Posts: 2
Joined: Tue Jun 11, 2013 8:20 am


Return to Linking Problem

Who is online

Users browsing this forum: No registered users and 1 guest