Interfacing GSL with LAPACK

Interfacing GSL with LAPACK

[Mark Hymers]

Hi,

I'm currently working on some neural network software using GSL and need
to do matrix inverses and eigenvalues.  I'm using ATLAS as my BLAS
library for speed.

The code I'm writing is in C++ (just to provide a nice interface to the
neural networks which it creates).

I have written the following code to wrap calls to invert a matrix:


//////////////////////////////////////////////////////////////////////
//
// inverse.cpp: routine to invert a square matrix; uses LAPACK
//
//////////////////////////////////////////////////////////////////////

extern "C" { 
int dgetri_ (int*, double*, int*, int*, double*, int *, int *); 
int dgetrf_ (int*, int*, double*, int*, int*, int*);
}

/**
 * \param matrix the address of a pointer to the data part of a matrix
 * \param size the size of the matrix
 * \return bool if false, the matrix was singular
 * 
 * Inverts the square matrix given to the function.
 *
 * <em> For example: </em> <br>
 * Assuming a is a 4 x 4 gsl_matrix: <br>
 * \code bool worked = matrix_inv(&a->data, (int)a->size1); \endcode
 * If the returned value is true, a will now contain the inverse 
 * of the original matrix.  If false, the matrix was singular.
 **/
bool matrix_inv(double **matrix, int size)
{
	int *info = new int();
	int index[size];
	double work[size*size];

	dgetrf_(&size,&size,&matrix[0][0],&size,index,info);

	if (*info)
	{
  		delete info;
		return false;
	}

	dgetri_(&size,&matrix[0][0],&size,index,work,&size,info);
	delete info;
	return true;
}


If I then call this with a small matrix (say 2,2 in the following case),
it works fine:


///////////////////////////////////////////////////////////////////////
//
// inv.cpp: example of GSL usage and how to use the generic routines we
//          provide
//
///////////////////////////////////////////////////////////////////////


#include "netlab/generic.h"

/*
	Should be compiled using:

	g++ -I../include inv.cpp ../src/lapack/lapack.o \
		-lgsl -lf77blas -lcblas -latlas -lg2c

	Note that you need to have compiled lapack.o first.
*/

int main()
{

	// Initialise a 2x2 matrix with known values
	Matrix2d *a = gsl_matrix_alloc(2,2);
	gsl_matrix_set(a, 0, 0, 1.0);
	gsl_matrix_set(a, 0, 1, 2.0);
	gsl_matrix_set(a, 1, 0, 4.0);
	gsl_matrix_set(a, 1, 1, 3.0);

	clock_t stime = std::clock();
	matrix_inv ( &a->data, (int)a->size1); 

	// You *MUST* free memory once it's finished with
	gsl_matrix_free(a);
	
	return 0;
}



However when I use a 4500x4500 matrix, the program runs fine on my P4
laptop but segmentation faults on our AMD Athlon processing machine (the
code is compiled with -march=pentium4/athlon for the appropriate case): 


///////////////////////////////////////////////////////////////////////
//
// inv2.cpp: example of GSL usage and how to use the generic routines we
//          provide
//
///////////////////////////////////////////////////////////////////////

#include "netlab/generic.h"
#include "netlab/twister.h"

/*
	Should be compiled using:

	g++ -I../include inv.cpp ../src/lapack/lapack.o \
		../lib/twister/out.o \
		-lgsl -lf77blas -lcblas -latlas -lg2c

	Note that you need to have compiled lapack.o 
	and generic.o (in src/lapack src/generic respectively)
	first.
*/

int main()
{

	Matrix2d *a = gsl_matrix_alloc(4500,4500);
	twister(a); // twister simply fills up the matrix with random values

	matrix_inv ( &a->data, (int)a->size1); 

	// You *MUST* free memory once it's finished with
	gsl_matrix_free(a);
	
	return 0;
}

Using gdb, I've tracked this down and it seems to happen with this call
in matrix_inv:
dgetrf_(&size,&size,&matrix[0][0],&size,index,info);

The strange thing is, according to the ASM the SegFault seems to occur
*before* the actual call to the routine.

Am I doing something incredibly stupid here?

Many thanks,

Mark

-- 
Mark Hymers, University of Newcastle Medical School
Intercalating Medical Student (BMedSci)

Re: Interfacing GSL with LAPACK

[Lukas Dobrek]

On Tue, Apr 08, 2003 at 04:59:25PM +0100, Mark Hymers wrote:
> Hi,
> 
> I'm currently working on some neural network software using GSL and need
> to do matrix inverses and eigenvalues.  I'm using ATLAS as my BLAS
> library for speed.
> 
> The code I'm writing is in C++ (just to provide a nice interface to the
> neural networks which it creates).
> 
> I have written the following code to wrap calls to invert a matrix:
> 	int *info = new int();
> 	int index[size];
> 	double work[size*size];

> The strange thing is, according to the ASM the SegFault seems to occur
> *before* the actual call to the routine.

Allocate it on heap not on stack. i.e.  using malloc or new.

Take Care
ld

-- 
"The United States, as the world knows, will never start a war. . . .
we shall also do our part to build a world of peace where the weak 
are safe and the strong are just."
									- John F. Kennedy
£ukasz Dobrek

Re: Interfacing GSL with LAPACK

[Mark Hymers]

On Wed, 09, Apr, 2003 at 09:08:34AM +0200, Lukas Dobrek spoke thus..
> > The strange thing is, according to the ASM the SegFault seems to occur
> > *before* the actual call to the routine.
> 
> Allocate it on heap not on stack. i.e.  using malloc or new.

Ahh.. Thanks - that makes sense.  I'm also switching to using CLAPACK
which makes the function definitions (for me) easier to read.

Many thanks again,

Mark

-- 
Mark Hymers, University of Newcastle Medical School
Intercalating Medical Student (BMedSci)

Re: Interfacing GSL with LAPACK

[Jeff Spirko]

On Tue, Apr 08, 2003 at 04:59:25PM +0100, Mark Hymers wrote:
> 	double work[size*size];
> 	dgetrf_(&size,&size,&matrix[0][0],&size,index,info);

You're declaring work as a 1-dimensional array, but referencing it
as a 2-dimensional array.  matrix[0] is ok as the first element, but
matrix[0][0] tries to use that first element as a pointer, hence the
segfault.

-- 
Jeff Spirko   spirko@lehigh.edu   spirko@yahoo.com   WD3V   |=>

The study of non-linear physics is like the study of non-elephant biology.

All theoretical chemistry is really physics;
and all theoretical chemists know it. -- Richard P. Feynman 

Re: Interfacing GSL with LAPACK

[Mark Hymers]

On Wed, 09, Apr, 2003 at 10:41:06AM -0400, Jeff Spirko spoke thus..
> On Tue, Apr 08, 2003 at 04:59:25PM +0100, Mark Hymers wrote:
> > 	double work[size*size];
> > 	dgetrf_(&size,&size,&matrix[0][0],&size,index,info);
> 
> You're declaring work as a 1-dimensional array, but referencing it
> as a 2-dimensional array.  matrix[0] is ok as the first element, but
> matrix[0][0] tries to use that first element as a pointer, hence the
> segfault.

Thanks.  There were so many mistakes in that code it's unbelievable :-)

Anyways, I've sorted it all now - thanks all for your help.

Mark

-- 
Mark Hymers, University of Newcastle Medical School
Intercalating Medical Student (BMedSci)