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 (5.28 KiB) Viewed 13174 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 **********************************************/
}
