Hi,

Firstly, sorry about the large amount of code in this question. OK, so I'm making a C++ wrapper for the gsl_matrix struct and its associated functions. To make my approach a bit easier to extend to other kinds of gsl structs, I have a template base class called gsl_base , which has a set of useful functions for getting at and checking the status of the underlying gsl struct when using the C++ wrapper class. So, gsl_base looks like this:

namespace gsl{

template< typename T >
class gsl_base
{
public:

    /// Get a pointer to the underlying GSL struct
    ///
    /// Will not throw, can return NULL
    inline T* ptr(){                    return M_pGSLData;  }

    /// Get a const pointer to the underlying GSL struct
    ///
    /// Will not throw, can return NULL
    inline const T* const_ptr() const{  return M_pGSLData;  }

    /// Check if it's OK to try and manipulate the object
    ///
    /// Will not throw
    inline bool hasValue() const{       return M_pGSLData != 0;  }

    /// Check if it's not OK to try and manipulate the object
    ///
    /// Will not throw
    inline bool isNull() const{         return M_pGSLData == 0; }

protected:

    gsl_base() : M_pGSLData( NULL ){}
    gsl_base( T* original ) : M_pGSLData( original ){}
    ~gsl_base(){}

    /// Set the value of the underlying GSL struct
    ///
    /// Will not throw
    inline void set_ptr( T* p ){        M_pGSLData = p;     }

    T* M_pGSLData;

};

}

And to make a simple C++ wrapper for gsl_matrix I inherit from gsl_base like this:

#pragma once

#include <stdexcept>
#include <cassert>
#include <iostream>

#include "../Common/number.h"
#include "../Common/gslstructbase.h"

#include <gsl/gsl_matrix.h>

namespace gsl
{

class realMatrix : public gsl::gsl_base< gsl_matrix > {
	
	public:
	
	typedef real				value_type;
	typedef size_t				size_type;
	typedef real* 				iterator;
	typedef const real*		const_iterator;
	typedef real&				reference;
	typedef const real&		const_reference;
	typedef real*				pointer;
	typedef const real*		const_pointer;
	typedef ptrdiff_t			difference_type;
	
	struct M_matrix_element_type{
		realMatrix::size_type row;
		realMatrix::size_type col;
		
		M_matrix_element_type( realMatrix::size_type r, realMatrix::size_type c ) :
		row( r ), col( c ) {}
	};
	
	typedef M_matrix_element_type element_type;

	/// Constructors
	realMatrix();
	realMatrix( size_type r, size_type c ) throw ( std::bad_alloc );
	~realMatrix(){}
	
	/// Print the addresses of the elements in the matrix (for debugging only!)
	/// to be removed when problems with memory allocation are fixed
	void printAddresses() const
	{
		std::cout << std::endl;
		for ( size_t i = 0; i < M_uRows; ++i ){
			for( size_t j = 0; j < M_uCols; ++j )
				std::cout << M_pStart + i*M_uCols + j << " ";
			std::cout << std::endl;
		}
	}

	private:
	
	pointer M_pStart;
	pointer M_pFinish;
	size_type M_uRows;
	size_type M_uCols;
	
	/// Set all private fields to correspond to a gsl_matrix structs values
	/// m must not be NULL
	///
	/// Will not throw
	inline void M_set_all_fields( const gsl_matrix* m )
	{
		assert( m != 0 );
		M_pStart = m->data;
		M_pFinish = M_pStart + m->size1 * m->size2;
		M_uRows = m->size1;
		M_uCols = m->size2;
	}
};

}

real is a typedef for double defined in ../Common/number.h . And the corresponding cpp file for reamMatrix is quite simple and looks like:

#include "matrix.h"

#include "../Utils/mathutils.h"

////////////////////////////////////////////////////////////

namespace gsl{

////////////////////////////////////////////////////////////

realMatrix::realMatrix() : gsl_base(0), M_pStart(0), M_pFinish(0), M_uRows(0), M_uCols(0) {}	

////////////////////////////////////////////////////////////

realMatrix::realMatrix( realMatrix::size_type r, realMatrix::size_type c )
throw( std::bad_alloc ) : gsl_base( gsl_matrix_alloc( atLeastOne( r ), atLeastOne( c ) ) )
{
	if ( isNull() )
		throw std::bad_alloc();
		
	M_set_all_fields( this->const_ptr() );
}

////////////////////////////////////////////////////////////

}	// End of gsl namespace

////////////////////////////////////////////////////////////

There are a bunch of other functions for doing various matrixy things, but these are the minimum needed to show the issue. OK, I also have an auto-test that was supposed to check a boolean equals that I had defined for gsl::realMatrix , but that's all been stripped away now and all it does is instantiate two matrices and print the addresses of their elements. (for diagnosing why I was getting segfaults and crash dumps). So, my test looks like:

void MatrixTestSuite::MatricesAreEqual()
{
    const int iRows = 10;
    const int iCols = 5;

    gsl::realMatrix m_1( iRows, iCols );
	 std::cout << "m_1 addresses = ";
	 m_1.printAddresses();
    gsl::realMatrix mat_2( iRows, iCols );
	 std::cout << "m_1 addresses = ";
	 m_1.printAddresses();
	 std::cout << "m_2 addresses = ";
	 mat_2.printAddresses();
}

That's it. Instantiate a matrix, print its element addresses. Instantiate a second matrix, print the address of the first matrix's elements again and then the second's elements. The results are not what I expected. I get something like this:

m_1 addresses = 
0x81bfd70 0x81bfd78 0x81bfd80 0x81bfd88 0x81bfd90 
0x81bfd98 0x81bfda0 0x81bfda8 0x81bfdb0 0x81bfdb8 
0x81bfdc0 0x81bfdc8 0x81bfdd0 0x81bfdd8 0x81bfde0 
0x81bfde8 0x81bfdf0 0x81bfdf8 0x81bfe00 0x81bfe08 
0x81bfe10 0x81bfe18 0x81bfe20 0x81bfe28 0x81bfe30 
0x81bfe38 0x81bfe40 0x81bfe48 0x81bfe50 0x81bfe58 
0x81bfe60 0x81bfe68 0x81bfe70 0x81bfe78 0x81bfe80 
0x81bfe88 0x81bfe90 0x81bfe98 0x81bfea0 0x81bfea8 
0x81bfeb0 0x81bfeb8 0x81bfec0 0x81bfec8 0x81bfed0 
0x81bfed8 0x81bfee0 0x81bfee8 0x81bfef0 0x81bfef8 
m_1 addresses = 
0x1 0x9 0x11 0x19 0x21 
0x29 0x31 0x39 0x41 0x49 
0x51 0x59 0x61 0x69 0x71 
0x79 0x81 0x89 0x91 0x99 
0xa1 0xa9 0xb1 0xb9 0xc1 
0xc9 0xd1 0xd9 0xe1 0xe9 
0xf1 0xf9 0x101 0x109 0x111 
0x119 0x121 0x129 0x131 0x139 
0x141 0x149 0x151 0x159 0x161 
0x169 0x171 0x179 0x181 0x189 
m_2 addresses = 
0x81bff08 0x81bff10 0x81bff18 0x81bff20 0x81bff28 
0x81bff30 0x81bff38 0x81bff40 0x81bff48 0x81bff50 
0x81bff58 0x81bff60 0x81bff68 0x81bff70 0x81bff78 
0x81bff80 0x81bff88 0x81bff90 0x81bff98 0x81bffa0 
0x81bffa8 0x81bffb0 0x81bffb8 0x81bffc0 0x81bffc8 
0x81bffd0 0x81bffd8 0x81bffe0 0x81bffe8 0x81bfff0 
0x81bfff8 0x81c0000 0x81c0008 0x81c0010 0x81c0018 
0x81c0020 0x81c0028 0x81c0030 0x81c0038 0x81c0040 
0x81c0048 0x81c0050 0x81c0058 0x81c0060 0x81c0068 
0x81c0070 0x81c0078 0x81c0080 0x81c0088 0x81c0090

As you can see, the elements of the first matrix are invalidated by the instantiation of the second!?!? This is bugging the hell out of me, so any suggestions would be really gratefully received! Also, the command that's being used to compile this code is: g++ -c "gsl++/Source/gslpp/Matrix/matrix.cpp" -g -pedantic -Wall -o ../../Matrix_matrix.o -I. -I/usr/local/include -I/usr/include Cheers.

Comments
Nice work.

What if you try std::cout << [B]M_pGSLData->data[/B] + i*[B]M_pGSLData->tda[/B] + j << " "; instead?

Edited 5 Years Ago by m4ster_r0shi: n/a

Hi,

Can you please check sizeof(real) and sizeof(double) is same or not, may be lame check but only invalid memory access is the reason that making your first matrix corrupted.

On a first glance, code looks so far fine, only potential invalid memory access in the above code is while giving output, so lets check the size of those two.

Second is, can you use valgrind to see where the invalid memory access is happening, valgrind should point you exactly where you are accessing the invalid memory.

it comes free with most linux distro.

Another thing I want to point out is: Why you are allocating memory in the intializatin list, its very dangerous and could produce potential memory leak which is very hard to debug, consider the following code:

class A{
     AnotherObject *a, *b;
public:
     A(): a(new AnotherObject()), b(new AnotherObject()){ //What will happen if initialization of a is successful but b fails, potential memory leak.
     }
}

Solution of this is smart pointers. if you store value in smart pointer, lets say you you store a in the smart pointer, even b fails, smart pointer will take care of memory deallocation of a.

Last but not least, I would suggest not to define any Exception specification except throw(), it doesn't do any good.

Edited 5 Years Ago by alwaysLearning0: n/a

I tried to recreate the problem but I couldn't. Does this work ok for you?

#include <iostream>
#include <stdexcept>
#include <cassert>
#include <gsl/gsl_matrix.h>

typedef double real;

namespace gsl{

template< typename T >
class gsl_base
{
public:

    /// Get a pointer to the underlying GSL struct
    ///
    /// Will not throw, can return NULL
    inline T * ptr() { return M_pGSLData; }

    /// Get a const pointer to the underlying GSL struct
    ///
    /// Will not throw, can return NULL
    inline const T * const_ptr() const { return M_pGSLData; }

    /// Check if it's OK to try and manipulate the object
    ///
    /// Will not throw
    inline bool hasValue() const { return M_pGSLData != 0; }

    /// Check if it's not OK to try and manipulate the object
    ///
    /// Will not throw
    inline bool isNull() const { return M_pGSLData == 0; }

protected:

    gsl_base() : M_pGSLData( NULL ) {}
    gsl_base( T * original ) : M_pGSLData( original ) {}
    ~gsl_base() {}

    /// Set the value of the underlying GSL struct
    ///
    /// Will not throw
    inline void set_ptr( T * p ){ M_pGSLData = p; }

    T * M_pGSLData;

};

}

namespace gsl
{

class realMatrix : public gsl::gsl_base< gsl_matrix >
{
public:

	typedef real              value_type;
	typedef size_t             size_type;
	typedef real *              iterator;
	typedef const real *  const_iterator;
	typedef real &             reference;
	typedef const real & const_reference;
	typedef real *               pointer;
	typedef const real *   const_pointer;
	typedef ptrdiff_t    difference_type;

	struct M_matrix_element_type
	{
		realMatrix::size_type row;
		realMatrix::size_type col;

		M_matrix_element_type( realMatrix::size_type r, realMatrix::size_type c ) :
            row( r ), col( c ) {}
	};

	typedef M_matrix_element_type element_type;

	/// Print the addresses of the elements in the matrix (for debugging only!)
	/// to be removed when problems with memory allocation are fixed
	void printAddresses() const
	{
		std::cout << std::endl;
		for ( size_t i = 0; i < M_uRows; ++i ){
			for( size_t j = 0; j < M_uCols; ++j )
				std::cout << M_pStart + i*M_pGSLData->tda + j << " ";
				//std::cout << M_pGSLData->data + i*M_pGSLData->tda + j << " ";
			std::cout << std::endl;
		}
	}

    /// Constructors
	realMatrix() : gsl_base<gsl_matrix>(0), M_pStart(0), M_pFinish(0), M_uRows(0), M_uCols(0) {}

	realMatrix( realMatrix::size_type r, realMatrix::size_type c ) throw( std::bad_alloc ) :
        gsl_base<gsl_matrix>( gsl_matrix_alloc(  r , c  ) )
    {
        if ( isNull() )
            throw std::bad_alloc();

        M_set_all_fields( this->const_ptr() );
    }

	~realMatrix(){}

private:

	pointer M_pStart;
	pointer M_pFinish;
	size_type M_uRows;
	size_type M_uCols;

	/// Set all private fields to correspond to a gsl_matrix structs values
	/// m must not be NULL
	///
	/// Will not throw
	inline void M_set_all_fields( const gsl_matrix* m )
	{
		assert( m != 0 );
		M_pStart = m->data;
		M_pFinish = M_pStart + m->size1 * m->size2;
		M_uRows = m->size1;
		M_uCols = m->size2;
	}
};

}

int main()
{
    const int iRows = 4;
    const int iCols = 5;

    gsl::realMatrix m_1( iRows, iCols );

    std::cout << "m_1 addresses = ";
    m_1.printAddresses();

    gsl::realMatrix mat_2( iRows, iCols );

    std::cout << "m_1 addresses = ";
    m_1.printAddresses();

    std::cout << "m_2 addresses = ";
    mat_2.printAddresses();

    std::cin.get();
    return 0;
}

Also, did you try std::cout << M_pGSLData->data + i*M_pGSLData->tda + j << " "; yet?
This can quickly tell you whether the problem is on the M_pStart pointer or somewhere deeper.

From the print out, it is clear that the stride is correct, so the use of M_pGSLData->tda is not going to solve the problem, although it is recommended to use that value and not assume another (like M_uCols). Clearly, the value of the pointer M_start is being set to the value 1, sometime between the first and second printout. First, I would not recommend that you duplicate the matrix information (like start-pointer, size1 and size2), this will inevitably require some assumption about how GSL deals with the data, or way too many updating of the state, and there is no tangible gain to doing this. Use the data in M_pGSLData directly as much as possible.

As for why your M_start pointer gets overwritten with the value 1, that's hard to tell without complete code. You must identity what objects are directly next to the m_1 object in memory and how they could write memory out-of-bound that would spill over to the m_1 object.

Also, remember that GSL is pretty notorious for having memory corruption problems because of lot of its code is hand-made translations from fortran by not-so-good programmers (good scientists, but not good programmers), which is very error-prone. We almost destroyed a robot once because of a memory corruption problem that was internal to some GSL code (guess why I don't ever "just use" GSL code anymore, I always carefully scrutinise it before).

Thanks for all the replies, I've been at work all day so I haven't had a chance to test anything out until now. Some replies:

What if you try std::cout << [B]M_pGSLData->data[/B] + i*[B]M_pGSLData->tda[/B] + j << " "; instead?

This doesn't make any difference, The problem is that the start position of the matrix data is being set to 0x1 , so the problem starts before I would even get to use the M_pGSLData->tda part.

Can you please check sizeof(real) and sizeof(double) is same or not...

As expected, they're the same (8)

Second is, can you use valgrind to see where the invalid memory access is happening, Valgrind should point you exactly where you are accessing the invalid memory.

I'll give that a go when I get a bit of time.

Why you are allocating memory in the intializatin list, its very dangerous and could produce potential memory leak

Do you think this is still the case here? I'm using gsl_matrix_alloc in the initialisation list, which leaves M_pGSLData equal to NULL if it's unsuccessful. In the body of the constructor, I'm immediately checking for a NULL pointer with the isNull() member of gsl_base and throwing if it wasn't.

Solution of this is smart pointers. if you store value in smart pointer, lets say you you store a in the smart pointer, even b fails, smart pointer will take care of memory deallocation of a

Yeah, a smart pointer is generally safer than using a raw pointer, but I didn't want to pay the (probably small) performance penalty here :o)

Last but not least, I would suggest not to define any Exception specification except throw(), it doesn't do any good.

I like to do it for documentation purposes mainly, so that I can see which exceptions are thrown just by looking at the declaration. Otherwise, I would end up just writing this information in a comment or something :o)

I tried to recreate the problem but I couldn't. Does this work ok for you?

Really, that's interesting. I'm starting to wonder if I need to re-install GCC or something?

First, I would not recommend that you duplicate the matrix information (like start-pointer, size1 and size2), this will inevitably require some assumption about how GSL deals with the data, or way too many updating of the state, and there is no tangible gain to doing this. Use the data in M_pGSLData directly as much as possible.

I did this originally, but I have half an eye on moving away from the gsl representation into using my own memory management at some point, so I'm trying to keep the details of the underlying gsl_matrix quite abstracted. Maybe it's not worth it at this stage!

As for why your M_start pointer gets overwritten with the value 1, that's hard to tell without complete code. You must identity what objects are directly next to the m_1 object in memory and how they could write memory out-of-bound that would spill over to the m_1 object.

To be honest, what I posted is pretty much the complete code. It's really an annoying problem, I have a vector version that is working without any issue at all (using the same underlying gsl_base base class). Maddening!

I was unaware of any memory corruption issues in gsl. I will have to be more careful in future!

OK, this problem was to do with having the constructor definition in a separate file to the declaration. If I define them in the class declaration (as in m4ster_r0shi's example) then it all works as expected. I think this has something to do with the templated inheritance trick. Any offers?

This article has been dead for over six months. Start a new discussion instead.