In answer to Julie, in the following the code that produce eigenvalues nulls.
I also follow the advice.
Installed gfortran and build lapack (fortran) e lapacke.
About lapacke I didn't find any documentations.
Could you provide me an example of simple solution of a complex double precision solution of a linear system?
This is just to understand how to use lapacke (it seems that I do not have modify the array storage -rows storage vs columns storage- etc.)
Thank you
Daniele
- Code: Select all
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "f2c.h"
#include "clapack.h"
#include "blaswrap.h"
#define dim 2
int main()
{
FILE *fp;
doublereal re,im;
int i,j,count;
/* Definitions for subtroutine */
char jobvl; // 'N' left eigv non computed; 'V' left eigv computed
char jobvr; // 'N' right eigv non computed; 'V' right eigv computed
integer n; // Order of the matrix
integer lda; // Leading deimension of array
doublereal *a; // Input array
doublereal *wr; // Output eigenvalues
doublereal *wi; // Output eigenvalues
doublereal *vl; // Stores, if required, left eigv, or not referenced
integer ldvl; // Leading dimension for VL
doublereal *vr; // Stores, if required, righr eigv, or not referenced
integer ldvr; // Leading dimension dor VR
integer lwork; // Auxiliary data
doublereal *work; // "
double *rwork; // "
integer info; // "
/* Definition parameters */
jobvl = 'N';
jobvr = 'V';
n = dim;
lda = dim;
wr = (doublereal *)malloc(sizeof(doublereal)*n);
if(wr==NULL)
{
printf("Error stop wr\n");
exit(1);
}
wi = (doublereal *)malloc(sizeof(doublereal)*n);
if(wi==NULL)
{
printf("Error stop wi\n");
exit(1);
}
ldvl = 1;
vl = (doublereal *)malloc(sizeof(doublereal)*ldvl*n);
if(vl==NULL)
{
printf("Error stop vl\n");
exit(1);
}
ldvr = n;
vr = (doublereal *)malloc(sizeof(doublereal)*ldvr*n);
if(vr==NULL)
{
printf("Error stop vr\n");
exit(1);
}
lwork = 8*n;
work = (doublereal *)malloc(sizeof(doublereal)*lwork);
if(work==NULL)
{
printf("Error stop work\n");
exit(1);
}
a = (doublereal *)malloc(sizeof(doublereal)*n*n);
if(a==NULL)
{
printf("Error stop a\n");
exit(1);
}
fp = fopen("in.dat","r");
count=0;
for(i=0; i<n*n; i++)
{
fscanf(fp,"%25.20lf", &re);
a[count] = re;
count++;
}
fclose(fp);
dgeev_(&jobvl, &jobvr, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info);
printf("INFO = %d \n",info);
for(i=0; i<n; i++)
{
printf(" EigenValue n. %d: Real = %le Imag = %le \n",i,wr[i], wi[i]);
}
count = 0;
for(i=0; i<n; i++)
{
printf(" EigenVecto n. %d:",i);
for(j=0; j<n; j++)
{
printf(" %25.20lf ",vr[count]);
count++;
}
printf("\n");
}
}