Shot in the dark: Anyone ever use CLAPACK routines?

Bruce Labitt bruce.labitt at myfairpoint.net
Sat May 22 12:22:09 EDT 2010


snip
> Not too much to report.  I even re-compiled ATLAS & LAPACK allowing gcc 
> & gfortran and got my example code to build.  Same problem with the 9x9 
> matrix.  The 2x2 double complex matrix svd worked! 
>
> I do have to say the interface to LAPACK is much better than CLAPACK.  
> (C or C++ calling FORTRAN)  I can go back to my old (bad) habits of 
> using bits of C++ to help make the code easier to follow.
>
> In numpy/scipy the code computed the svd with no problem.
>
> Not too much activity at the lapack-forum.  :(  Since the academic term 
> ended recently at UT Knoxville, maybe everyone is on vacation...
>
> -Bruce
>
>   
Well, I have something to try...  If one actually RTFM or in this case 
the readme, one finds the following:

# As a final point, we must stress that there is a difference in the definition
# of a two-dimensional array in Fortran and C.
# A two-dimensional Fortran array declared as
#
#    DOUBLE PRECISION A(LDA, N)
#
# is a contiguous piece of LDA X N double-words of memory, stored in
# column-major order: elements in a column are contiguous, and elements
# within a row are separated by a stride of LDA double-words.
#
# In C, however, a two-dimensional array is in row-major order.  Further, the
# rows of a two-dimensional C array need not be contiguous.  The array 
#
#    double A[LDA][N];
#
# actually has LDA pointers to rows of length N.  
# These pointers can in principle be anywhere in memory. Passing such a
# two-dimensional C array to a CLAPACK routine will almost surely give
# erroneous results.
#
# Instead, you must use a one-dimensional C array of size LDA X N
# double-words (or else malloc the same amount of space).
# We recommend using the following code to get the array CLAPACK will be 
# expecting:
#
#    double *A;
#    A = malloc( LDA*N*sizeof(double) );

Umm, yeah, I can do that!  My original code had that...  For some reason I got rid of it.
So a SEG FAULT makes some sense for my 9x9.  Because in C the array (matrix) does not have 
to be contiguous in memory.  But FORTRAN is assuming a column major contiguous array.  So it
seems likely I could access unallocated memory, and hence a seg fault.  

OK, I can relax now...  Saga continues on Monday...  




More information about the gnhlug-discuss mailing list