Interfacing GSL with LAPACK

[Mark Hymers <mark dot hymers at ncl dot ac dot uk>]

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)