bpmnr/gsl_matrix.c

Go to the documentation of this file.

#include <bpm/bpm_messages.h>
#include <bpm/bpm_nr.h>

int gsl_matrix_swap_columns(gsl_matrix * m, const size_t i, const size_t j) {
  
  const size_t size1 = m->size1;
  const size_t size2 = m->size2;

  if (i >= size2)
    {
      bpm_error("first column index is out of range in gsl_matrix_swap_columns(...)",
                __FILE__, __LINE__ );
      return BPM_FAILURE;
    }

  if (j >= size2)
    {
      bpm_error("second column index is out of range in gsl_matrix_swap_columns(...)",
                __FILE__, __LINE__ );
      return BPM_FAILURE;
    }

  if (i != j)
    {
      double *col1 = m->data + i;
      double *col2 = m->data + j;
      
      size_t p;
      
      for (p = 0; p < size1; p++)
        {
          size_t k;
          size_t n = p * m->tda;
 
          for (k = 0; k < 1; k++)
            {
              double tmp = col1[n+k] ;
              col1[n+k] = col2[n+k] ;
              col2[n+k] = tmp ;
            }
        }
    }

  return BPM_SUCCESS;
}


// -------------------------------------------------------------------

_gsl_vector_view gsl_matrix_column (gsl_matrix * m, const size_t j) {
 
  _gsl_vector_view view = NULL_VECTOR_VIEW;
  
  if (j >= m->size2)
    {
      bpm_error("column index is out of range in gsl_matrix_column", 
                __FILE__, __LINE__ );
      return view;
    }

  {
    gsl_vector v = NULL_VECTOR;
    
    v.data = m->data + j;
    v.size = m->size1;
    v.stride = m->tda;
    v.block = m->block;
    v.owner = 0;

    view.vector = v;
    return view;
  }
}

double gsl_matrix_get(const gsl_matrix * m, const size_t i, const size_t j) {

  return *(double*) (m->data + (i * m->tda + j));
}




void gsl_matrix_set(gsl_matrix * m, const size_t i, const size_t j, const double x) { 

  *(double*) (m->data + (i * m->tda + j)) = x;

}



_gsl_matrix_view gsl_matrix_submatrix (gsl_matrix * m, const size_t i, const size_t j, 
                                       const size_t n1, const size_t n2) {

  _gsl_matrix_view view = NULL_MATRIX_VIEW; 

  if (i >= m->size1)
    {
      bpm_error("row index is out of range in gsl_matrix_submatrix(...)",
                __FILE__, __LINE__ );
      return view;
    }
  else if (j >= m->size2)
    {
      bpm_error("column index is out of range in gsl_matrix_submatrix(...)",
                __FILE__, __LINE__ );
      return view;
    }
  else if (n1 == 0)
    {
      bpm_error("first dimension must be non-zero in gsl_matrix_submatrix(...)",
                __FILE__, __LINE__ );
      return view;
    }
  else if (n2 == 0)
    {
      bpm_error("second dimension must be non-zero in gsl_matrix_submatrix(...)",
                __FILE__, __LINE__ );
      return view;
    }
  else if (i + n1 > m->size1)
    {
      bpm_error("first dimension overflows matrix in gsl_matrix_submatrix(...)",
                __FILE__, __LINE__ );
      return view;
    }
  else if (j + n2 > m->size2)
    {
      bpm_error("second dimension overflows matrix in gsl_matrix_submatrix(...)",
                __FILE__, __LINE__ );
      return view;
    }

  {
     gsl_matrix s = NULL_MATRIX;

     s.data = m->data + (i * m->tda + j);
     s.size1 = n1;
     s.size2 = n2;
     s.tda = m->tda;
     s.block = m->block;
     s.owner = 0;
     
     view.matrix = s;     
     return view;
  }
}

gsl_matrix * gsl_matrix_alloc (const size_t n1, const size_t n2)
{

  gsl_block * block;
  gsl_matrix * m;

  if (n1 == 0)
    {
      bpm_error("matrix dimension n1 must be positive integer in gsl_matrix_alloc(...)",
                __FILE__, __LINE__);
      return NULL;
    }
  else if (n2 == 0)
    {
      bpm_error("matrix dimension n2 must be positive integer in gsl_matrix_alloc(...)",
                __FILE__, __LINE__);
      return NULL;
    }

  m = (gsl_matrix *) malloc (sizeof(gsl_matrix));

  if (m == 0)
    {
      bpm_error("failed to allocate space for matrix struct in gsl_matrix_alloc(...)",
                __FILE__, __LINE__);
      return NULL;
    }

  /* FIXME: n1*n2 could overflow for large dimensions */

  block = gsl_block_alloc (n1 * n2) ;

  if (block == 0)
    {
      bpm_error("failed to allocate space for block in gsl_matrix_alloc(...)",
                __FILE__, __LINE__);
      return NULL;
    }

  m->data = block->data;
  m->size1 = n1;
  m->size2 = n2;
  m->tda = n2; 
  m->block = block;
  m->owner = 1;

  return m;
}

gsl_matrix * gsl_matrix_calloc (const size_t n1, const size_t n2)
{
  size_t i;

  gsl_matrix * m = gsl_matrix_alloc( n1, n2);

  if (m == 0)
    return 0;

  /* initialize matrix to zero */

  for (i = 0; i < n1 * n2; i++)
    {
      m->data[i] = 0;
    }

  return m;
}

_gsl_vector_const_view gsl_matrix_const_row (const gsl_matrix * m, const size_t i)
{

  _gsl_vector_const_view view = NULL_VECTOR_VIEW;
  
  if (i >= m->size1)
    {
      bpm_error("row index is out of range in gsl_matrix_const_row(...)",
                __FILE__, __LINE__);
      return view;
    }
  
  {
    gsl_vector v = NULL_VECTOR;
    
    v.data = m->data + i * m->tda;
    v.size = m->size2;
    v.stride = 1;
    v.block = m->block;
    v.owner = 0;
    
    view.vector = v;
    return view;
  }
}

_gsl_vector_view gsl_matrix_row (gsl_matrix * m, const size_t i)
{

  _gsl_vector_view view = NULL_VECTOR_VIEW;
  
  if (i >= m->size1)
    {
      bpm_error("row index is out of range in gsl_matrix_row(...)",
                __FILE__, __LINE__);
      return view;
    }
  
  {
    gsl_vector v = NULL_VECTOR;
    
    v.data = m->data + i * m->tda;
    v.size = m->size2;
    v.stride = 1;
    v.block = m->block;
    v.owner = 0;
    
    view.vector = v;
    return view;
  }
}

_gsl_vector_const_view gsl_matrix_const_column (const gsl_matrix * m, const size_t i)
{

  _gsl_vector_const_view view = NULL_VECTOR_VIEW;
  
  if (i >= m->size1)
    {
      bpm_error("row index is out of range", __FILE__, __LINE__ );
      return view;
    }
  
  {
    gsl_vector v = NULL_VECTOR;
    
    v.data = m->data + i * m->tda;
    v.size = m->size2;
    v.stride = 1;
    v.block = m->block;
    v.owner = 0;
    
    view.vector = v;
    return view;
  }
}

void gsl_matrix_set_identity (gsl_matrix * m)
{

  size_t i, j;
  double * const data = m->data;
  const size_t p = m->size1 ;
  const size_t q = m->size2 ;
  const size_t tda = m->tda ;

  const double zero = 0.;
  const double one = 1.;

  for (i = 0; i < p; i++)
    {
      for (j = 0; j < q; j++)
        {
          *(double *) (data + (i * tda + j)) = ((i == j) ? 1. : 0.);
        }
    }
}

---