Vector and matrix views (and bug report)

[Fabrice Rossi <rossi at ufrmd dot dauphine dot fr>]

Hi.

I think I've found a bug in gsl_vector_alloc_from_block (init_source.c in
vector directory). If I'm correct, the line 94, which contains v->block =
block; should be rewritten v->block = 0; If the block is not set to zero, it
will be deallocated by gsl_vector_free, which should not happen (the block
must be manualy deallocated with gsl_block_free).

I've started to implement a very rough neural network (a two layer perceptron)
with gls. It's much easier to do this is the network parameters can be
observed sometimes as a big vector (in order to represent the modeling error
done by the network as a multidimensionnal function) and sometimes as set of
matrices and vectors. This is easy to do internaly thanks to alloc_from_block. 

The only trouble is when I'm going to use the network with the
multidimensionnal minimization algorithms. I will receive a parameter vector
externally allocated, and I cannot spend time in allocation for the internal
representation. What I need is an extension to
gsl_vector_view_row_from_matrix, that do not perform memory allocation (it
would be even better if range_checking in this method can be removed with the
GSL_RANGE_CHECK_OFF trick). Basically, I want to be able to view part of a
vector as a matrix and part of a vector as another vector. Here is a proposed
implementation (adapted from alloc_from_block and from row_col):

int
gsl_matrix_view_from_vector(gsl_matrix * m,gsl_vector * base,size_t offset, 
			    size_t n1, size_t n2, size_t d2)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (n1 == 0)
    {
      GSL_ERROR_RETURN ("matrix dimension n1 must be positive integer",
			GSL_EDOM, 0);
    }
  else if (n2 == 0)
    {
      GSL_ERROR_RETURN ("matrix dimension n2 must be positive integer",
			GSL_EDOM, 0);
    }
  else if (d2 < n2)
    {
      GSL_ERROR_RETURN ("matrix dimension d2 must be greater than n2",
			GSL_EDOM, 0);
    }
  else if (base->size < offset + n1 * d2)
    {
      GSL_ERROR_RETURN ("matrix size exceeds available vector size",
			GSL_EDOM, 0);
    }
#endif
  if (m->block != 0)
    {
      GSL_ERROR ("matrix already has memory allocated to it", GSL_ENOMEM);
    }

  m->data = base->data + offset;
  m->size1 = n1;
  m->size2 = n2;
  m->dim2 = d2;
  
  return GSL_SUCCESS;
}

int
gsl_vector_view_from_vector(gsl_vector * v,gsl_vector * base,size_t offset, 
			    size_t n,size_t stride)
{
#ifndef GSL_RANGE_CHECK_OFF
  if (n == 0)
    {
      GSL_ERROR_RETURN ("vector length n must be positive integer",
			GSL_EDOM, 0);
    }

  if (stride == 0)
    {
      GSL_ERROR_RETURN ("stride must be positive integer", GSL_EDOM, 0);
    }

  if (base->size <= offset + (n - 1) * stride)
    {
      GSL_ERROR_RETURN ("vector would extend past end of vector", GSL_EDOM,
0);
    }
#endif

  if (v->block != 0)
    {
      GSL_ERROR ("vector already has memory allocated to it", GSL_ENOMEM);
    }

  v->data = base->data + offset ;
  v->size = n;
  v->stride = stride;

  return GSL_SUCCESS;
}

Fabrice Rossi