Segmentation Fault - clapack / vecLIB

[relain]

Hi,
So i've been trying to use clapack for a finite differences Quantum Mechanics problem, this boils down to eigenvaluing a very sparse non-symmetric matrix. I've been using the dgherd_ and dhseqr_ routines and have been getting good results, however when the order of my matrices (doubles) gets above 27^2 i get a segmentation fault which i simply cannot track down (if it's 26^2 it works out fine...).

The code is available at: http://www.relain.co.uk/stuff/physics/sqhfix.c

I'm compiling the source with the Mac Os 10.3.9 version of gcc and using the clapack and cblas libraries that are included within the power-pc optimised framework called veclib.

I'm running this on a G4 (32bit) with 768mb of ram, am i simply running out of memory? Surely not.

I'm more of a physicist than a computer scientist sooo maybe it's something very basic with my coding? Do i need to be dynamically allocating memory for the arrays?

Thanks for your help,

also if anyone out there is actually an author of lapack i think it's great!

[Julie]

Relain,

Your problem is not linked to LAPACK library. If you remove the LAPACK calls, your code crashes as well.
You should try to solve this problem first.

Once this problem has been solved, here is one way to proceed

- first call the LAPACK routine you want with the lwork parameter set to -1 in input (and all the other correctly set)
The workspace needed is returned in work[0]
- then free(work) and then re-allocate work of size work[0]
- you are ready to call LAPACK's driver

Code: Select all

 
int n,info,ipone=1,lwork;
double *A,*tau,*work;
....
/* workspace query */
work = (double *)malloc(n*sizeof(double)) ;
if (work==NULL){ printf("error of memory allocation\n"); exit(0); }
lwork = -1;
dgehrd_(&n, &one ,&n, A, &n, tau, work, &lwork, &info);
lwork = (int) work[0];
free(work);

/* allocation of the workspace */
work = (double *)malloc(lwork*sizeof(double)) ;
if (work==NULL){ printf("error of memory allocation\n"); exit(0); }

/* call to dgehrd (for example) */
dgehrd_(&n, &one ,&n, A, &n, tau, work, &lwork, &info);
free(work);
.....


Sincerely
Julie Langou

[relain]

worked like a treat, thank you very very much!

Now if only it was faaaaster...

[Julie]

relain,

glad for you that it works fine.

By the way , I went quickly through your code and was just wondering why you are not directly calling DGEEV. Just wondering if you were aware it existed.

Since you are using veclib, there's not much you can do about speed, your can buy a faster processor maybe...

I made some run with VecLib on a Mac G5 [ 2 * 2.5 Ghz with 2G of Ram] for dgeev. Here are my results:

Code: Select all

DGEEV [n=100 JOBVL=N JOBVR=N] 0.012145s
DGEEV [n=500 JOBVL=N JOBVR=N] 1.143880s
DGEEV [n=1000 JOBVL=N JOBVR=N] 9.702637s
DGEEV [n=1500 JOBVL=N JOBVR=N] 36.882315s
DGEEV [n=100 JOBVL=V JOBVR=V] 0.026176s
DGEEV [n=500 JOBVL=V JOBVR=V] 3.368334s
DGEEV [n=1000 JOBVL=V JOBVR=V] 34.054976s
DGEEV [n=1500 JOBVL=V JOBVR=V] 170.101812s

Sincerely
Julie