title: Better memory management for GSL gsl_vector
body: <p>When I write linear algebra programs in C++, I use the Armadillo library. It is based on templates, which provide me a way to define vectors of any length that don't necessarily require additional memory allocation, since they are statically assigned an appropriate memory buffer at compile-time. When I use <code>arma::Col&lt;double&gt;::fixed&lt;3&gt;</code> the compiler creates a "new type" on the fly so that the vector contains a buffer of exactly 3 doubles.</p>
  
  <p>Now I'm working on a linear algebra program in C, and I'm using the GNU Scientific Library (GSL). In order to instantiate a 3D vector, I do: <code>gsl_vector_alloc(3)</code> that returns a <code>gsl_vector*</code>. The problem is that this operation causes the dynamic allocation of a small portions of memory, and this happens millions of times during program runtime. My program is wasting a lot of resources to perform tens of millions of <code>malloc</code>/<code>free</code> operations.</p>
  
  <p>I've inspected the internal structure of <code>gsl_vector</code>:</p>
  
  <pre><code>typedef struct
  {
     size_t size;
     size_t stride;
     double * data;
     gsl_block * block;
     int owner;
  } gsl_vector;
  </code></pre>
  
  <p>For the library to work correctly, <code>data</code> should point to the first element of the vector, usually inside a <code>gsl_block</code> structure like this:</p>
  
  <pre><code>typedef struct
  {
     size_t size;
     double * data;
  } gsl_block;
  </code></pre>
  
  <p>which contains another <code>data</code> pointer. So, for instantiating a simple 3D vector, this sequence of <code>malloc</code>s happen:</p>
  
  <ol>
  <li>A <code>gsl_vector</code> structure is <code>malloc</code>'d (something around 40 bytes on x86_64).</li>
  <li>A <code>gsl_block</code> structure is <code>malloc</code>'d (16 bytes) and the <code>block</code> pointer of the gsl_vector is set to the memory address of the gsl_block just allocated</li>
  <li>An array of 3 doubles is <code>malloc</code>'d and its memory address is assigned to both <code>data</code> pointers (the one in <code>gsl_block</code> and the one in <code>gsl_vector</code>).</li>
  </ol>
  
  <p>I obtained a 40% performance gain by removing two <code>malloc</code>s. I created my custom <code>gsl_vector</code> creation routine, that allocates an array of 3 doubles and sets the <code>data</code> pointer of the <code>gsl_vector</code> to the address of this array. Then I return a <code>gsl_vector</code> (not a pointer).<br>
  But doing so, I still get millions of <code>malloc(3 * sizeof(double))</code> operations.</p>
  
  <p>I didn't manage to "embed" the array of 3 doubles inside the <code>gsl_vector</code> struct, since if the <code>data</code> pointer points to something which is <em>inside</em> the struct itself (hacky!), then the pointer is no longer valid when the vector is copied elsewhere!</p>
  
  <p>Do you have any ideas (apart switching to C++ or rolling my own linear algebra library)? I am open to any suggestion.</p>