fp_gsl_matrix.cc

Go to the documentation of this file.

#include "fp_gsl_matrix.h"
#include "linear_algebra.h"
#include <cstdlib>

using namespace FullPhysics;

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

void GslMatrix::reset(gsl_matrix* M, bool Owned)
{
  if(gsl_matrix_ && block_owned_)
    free(gsl_matrix_->block);
  if(gsl_matrix_ && owned_)
    gsl_matrix_free(gsl_matrix_);
  gsl_matrix_ = M;
  blitz::Array<double, 2> a(M->data, blitz::shape(M->size1, M->size2), 
                            blitz::shape(M->tda, 1), blitz::neverDeleteData);
  blitz_array_.reference(a);
  owned_ = Owned;
  block_owned_ = false;
}

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

void GslMatrix::reset(blitz::Array<double, 2>& M)
{
  if(gsl_matrix_ && block_owned_)
    free(gsl_matrix_->block);
  if(gsl_matrix_ && owned_)
    gsl_matrix_free(gsl_matrix_);
  blitz_array_.reference(to_c_order(M));
  gsl_matrix_ = (gsl_matrix*) malloc(sizeof(gsl_matrix));
  gsl_matrix_->size1 = blitz_array_.rows();
  gsl_matrix_->size2 = blitz_array_.cols();
  gsl_matrix_->tda = blitz_array_.cols();
  gsl_matrix_->data = blitz_array_.data();
  gsl_matrix_->owner = 0;
  gsl_matrix_->block = (gsl_block*) malloc(sizeof(gsl_block));
  gsl_matrix_->block->size = blitz_array_.size();
  gsl_matrix_->block->data = blitz_array_.data();
  owned_ = true;
  block_owned_ = true;
}

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

GslMatrix::~GslMatrix()
{
  if(gsl_matrix_ && block_owned_)
    free(gsl_matrix_->block);
  if(gsl_matrix_ && owned_)
    gsl_matrix_free(gsl_matrix_);
}

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

void GslVector::reset(gsl_vector* M, bool Owned)
{
  if(gsl_vector_ && block_owned_)
    free(gsl_vector_->block);
  if(gsl_vector_ && owned_)
    gsl_vector_free(gsl_vector_);
  gsl_vector_ = M;
  blitz::Array<double, 1> a(M->data, blitz::shape(M->size), 
                            blitz::shape(M->stride), blitz::neverDeleteData);
  blitz_array_.reference(a);
  owned_ = Owned;
  block_owned_ = false;
}

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

void GslVector::reset(blitz::Array<double, 1>& M)
{
  if(gsl_vector_ && block_owned_)
    free(gsl_vector_->block);
  if(gsl_vector_ && owned_)
    gsl_vector_free(gsl_vector_);
  blitz_array_.reference(to_c_order(M));
  gsl_vector_ = (gsl_vector*) malloc(sizeof(gsl_vector));
  gsl_vector_->size = blitz_array_.rows();
  gsl_vector_->stride = 1;
  gsl_vector_->data = blitz_array_.data();
  gsl_vector_->owner = 0;
  gsl_vector_->block = (gsl_block*) malloc(sizeof(gsl_block));
  gsl_vector_->block->size = blitz_array_.size();
  gsl_vector_->block->data = blitz_array_.data();
  owned_ = true;
  block_owned_ = true;
}

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

GslVector::~GslVector()
{
  if(gsl_vector_ && block_owned_)
    free(gsl_vector_->block);
  if(gsl_vector_ && owned_) 
    gsl_vector_free(gsl_vector_);
}