Page 1 of 1

Segmentation Fault - clapack / vecLIB

PostPosted: Sat Nov 05, 2005 7:40 pm
by relain
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:

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!

PostPosted: Tue Nov 08, 2005 2:30 pm
by Julie

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];

/* 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);

Julie Langou

PostPosted: Tue Nov 08, 2005 7:12 pm
by relain
worked like a treat, thank you very very much!

Now if only it was faaaaster...

PostPosted: Wed Nov 09, 2005 12:12 pm
by Julie

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