matrix_ulong.hpp

Go to the documentation of this file.

/*
 * Copyright (C) 2010, 2011, 2012, 2019, 2020, 2021, 2024 John D Lamb
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */
 
#ifndef CCGSL_MATRIX_ULONG_HPP
#define CCGSL_MATRIX_ULONG_HPP
 
#include<gsl/gsl_matrix.h>
#include<new>
#include<gsl/gsl_permute_matrix_ulong.h>
 
#include"exception.hpp"
#include"vector_ulong.hpp"
#include"permutation.hpp"
 
// This file is autogenerated
 
namespace gsl {
  class matrix_ulong {
  public:
    matrix_ulong(){
      owns_data = false;
      ccgsl_pointer = 0;
      count = 0; // initially nullptr will do
    }
    // Refines random access container
    // Refines assignable
    explicit matrix_ulong( size_t const n1, size_t const n2 ) : owns_data(true){
      if( n1 > 0 and n2 > 0 ) ccgsl_pointer = gsl_matrix_ulong_alloc( n1, n2 );
      else { ccgsl_pointer = new gsl_matrix_ulong; ccgsl_pointer->size1 = n1;
    ccgsl_pointer->size2 = n2; ccgsl_pointer->data = 0; }
      // just plausibly we could allocate matrix_ulong but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
    // try to tidy up before rethrowing
    if( n1 > 0 and n2 > 0 ) gsl_matrix_ulong_free( ccgsl_pointer );
    else delete ccgsl_pointer;
    throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
    }
    explicit matrix_ulong( gsl_matrix_ulong* v ) : owns_data(true){
      ccgsl_pointer = v;
      // just plausibly we could fail to allocate count: no further action needed.
      count = new size_t;
      *count = 1; // initially there is just one reference to ccgsl_pointer
    }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
    matrix_ulong( std::initializer_list<std::initializer_list<unsigned long> > initializer_list ) : owns_data(true){
      size_t const n1 = initializer_list.size();
      size_t const n2 = initializer_list.begin()->size();
      for( auto l : initializer_list ){
    if( l.size() != n2 ){
      gsl::exception e( "matrix_ulong rows have unequal sizes", __FILE__, __LINE__,
                exception::GSL_EBADLEN );
      throw( e );
    }
      }
      ccgsl_pointer = gsl_matrix_ulong_alloc( n1, n2 );
      // just plausibly we could allocate matrix_ulong but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
    // try to tidy up before rethrowing
    if( n1 > 0 and n2 > 0 ) gsl_matrix_ulong_free( ccgsl_pointer );
    else delete ccgsl_pointer;
    throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
      // now copy
      int r = 0;
      for( auto row : initializer_list ){
    size_t c = 0;
    for( auto x : row ){ set( r, c, x ); ++c; }
    ++r;
      }
    }
#endif
    // copy constructor
    matrix_ulong( matrix_ulong const& v ) : owns_data(true),
                        ccgsl_pointer( v.ccgsl_pointer ), count( v.count ){
      if( count != 0 ) ++*count; // matrix_ulong is now shared.
    }
    // assignment operator
    matrix_ulong& operator=( matrix_ulong const& v ){
      // first, possibly delete anything pointed to by this
      if( count == 0 or --*count == 0 ){
    if( ccgsl_pointer != 0 and ccgsl_pointer->size1 > 0 and ccgsl_pointer->size2 > 0 ) gsl_matrix_ulong_free( ccgsl_pointer );
    else delete ccgsl_pointer;
    delete count;
      }
      // Then copy
      owns_data = v.owns_data;
      ccgsl_pointer = v.ccgsl_pointer;
      count = v.count;
      if( count != 0 ) ++*count; // block_ulong is now shared.
      return *this;
    }
    // clone()
    matrix_ulong clone() const {
      matrix_ulong copy( get()->size1, get()->size2 );
      // Now copy
      gsl_matrix_ulong_memcpy( copy.get(), get() );
      // return new object
      return copy;
    }
    // destructor
    ~matrix_ulong(){
      if( count != 0 and --*count == 0 ){
    // could have allocated null pointer
    if( ccgsl_pointer != 0 and owns_data ){
      if( ccgsl_pointer->size1 > 0 and ccgsl_pointer->size2 > 0 )
        gsl_matrix_ulong_free( ccgsl_pointer );
      else delete ccgsl_pointer; }
    delete count;
      }
    }
 
    // Allow possibility of assigning from gsl_matrix_ulong without sharing
    void wrap_gsl_matrix_ulong_without_ownership( gsl_matrix_ulong* v ){
      if( count != 0 and --*count == 0 ){
    // could have allocated null pointer
    if( ccgsl_pointer != 0 ){
      if( ccgsl_pointer->size1 != 0 and ccgsl_pointer->size1 != 0 )
        gsl_matrix_ulong_free( ccgsl_pointer );
      else delete ccgsl_pointer; }
      }
      ccgsl_pointer = v;
      if(0 == count) count = new size_t;
      *count = 1;
      owns_data = false; // should never be able to delete ccgsl_pointer
    }
    void reset(){ matrix_ulong().swap( *this ); }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
    matrix_ulong( matrix_ulong&& v ) : owns_data {v.owns_data},
                       ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
      std::swap( count, v.count );
      v.ccgsl_pointer = nullptr;
    }
    matrix_ulong& operator=( matrix_ulong&& v ){
      matrix_ulong( std::move( v ) ).swap( *this );
      return *this;
    }
#endif
  private:
    struct vector_ulong_ptr : public vector_ulong {
      vector_ulong_ptr( vector_ulong const& v ) : vector_ulong( v ){}
      vector_ulong& operator*(){ return *this; }
      vector_ulong* operator->(){ return this; }
    };
    template<typename container,typename content,bool reverse_t> class iterator_base {
      friend class vector_ulong;
    public:
      typedef std::bidirectional_iterator_tag iterator_category;
      typedef vector_ulong value_type;
      typedef vector_ulong_ptr pointer;
      typedef value_type reference;
    public:
      // // operator*
      reference operator*() const {
    // Always try to return something
    static content something;
    // First check that iterator is initialised.
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return something;
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "matrix_ulong not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return something;
    }
    // Check that position make sense
    if( position >= v->size1() ){
      gsl_error( "trying to dereference end()", __FILE__, __LINE__, exception::GSL_EFAILED );
      return something;
    }
    // position and v are valid: return data
    return v->row( position );
      }
      // // operator->
      pointer operator->() const {
    // Always try to return something
    static content something_base;
    static vector_ulong_ptr something( something_base );
    // First check that iterator is initialised.
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return something;
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "matrix_ulong not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return something;
    }
    // Check that position make sense
    if( position >= v->size1() ){
      gsl_error( "trying to dereference end()", __FILE__, __LINE__, exception::GSL_EFAILED );
      return something;
    }
    // position and v are valid: return data
    vector_ulong_ptr ptr( v->row( position ) );
#ifdef __GXX_EXPERIMENTAL_CXX0X__
    return std::move( ptr );
#else
    return ptr;
#endif
      }
      bool operator==( iterator_base<container,content,reverse_t> const& i ) const {
    return this->v == i.v and this->position == i.position;
      }
      bool operator!=( iterator_base<container,content,reverse_t> const& i ) const {
    return not this->operator==( i );
      }
    protected:
      void increment(){
    // Only makes sense if v points to a matrix_ulong
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "matrix_ulong not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    }
    // increment position and check against size
    if( reverse_t ){ if( position >= 0 ) --position; }
    else { if( position < v->size1() ) ++position; }
      }
      void decrement(){
    // Only makes sense if v points to a matrix_ulong
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "matrix_ulong not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    }
    // decrement position and check against size
    if( reverse_t ){ if( position < v->size1() ) ++position; }
    else { if( position >= 0 ) --position; }
      }
      iterator_base(){ v = 0; }
      iterator_base( container* v, size_t position )
    : v( v ), position( position ){}
      container* v;
      size_t position;
    };
    // Need to know about const_iterator_t
    template<bool reverse_t> class const_iterator_t;
    template<bool reverse_t> class iterator_t : public iterator_base<matrix_ulong,vector_ulong,reverse_t>{
    public:
      // // Refines output iterator
      // // operator=
      iterator_t<reverse_t>& operator=( iterator_t<reverse_t> const& i ){
    iterator_base<matrix_ulong,vector_ulong,reverse_t>::v = i.v;
    iterator_base<matrix_ulong,vector_ulong,reverse_t>::position = i.position;
    return *this;
      }
      // // Refines forward iterator
      // // operator++ (both)
      iterator_t<reverse_t>& operator++(){
    iterator_base<matrix_ulong,vector_ulong,reverse_t>::increment();
    return *this;
      }
      iterator_t<reverse_t> operator++( int ){
    // return value
    iterator_t<reverse_t> result( *this );
    iterator_base<matrix_ulong,vector_ulong,reverse_t>::increment();
    return result;
      }
      // // Refines bidirectional iterator
      // // operator-- (both)
      iterator_t<reverse_t>& operator--(){
    iterator_base<matrix_ulong,vector_ulong,reverse_t>::decrement();
    return *this;
      }
      iterator_t<reverse_t> operator--( int ){
    // return value
    iterator_t<reverse_t> result( *this );
    iterator_base<matrix_ulong,vector_ulong,reverse_t>::decrement();
    return result;
      }
      bool operator==( const_iterator_t<reverse_t> const& i ) const {
        return this->v == i.v and this->position == i.position;
      }
      bool operator!=( const_iterator_t<reverse_t> const& i ) const {
        return not this->operator==( i );
      }
      iterator_t() : iterator_base<matrix_ulong,vector_ulong,reverse_t>(){}
    protected:
      friend class matrix_ulong;
      // We need a constructor for matrix_ulong
      iterator_t( matrix_ulong* v, size_t position )
      : iterator_base<matrix_ulong,vector_ulong,reverse_t>( v, position ){}
    };
    template<bool reverse_t> class const_iterator_t
      : public iterator_base<matrix_ulong const,vector_ulong const,reverse_t>{
    public:
      // // Refines output iterator
      // // operator=
      const_iterator_t<reverse_t>& operator=( const_iterator_t<reverse_t> const& i ){
    iterator_base<matrix_ulong const,vector_ulong const,reverse_t>::v = i.v;
    iterator_base<matrix_ulong const,vector_ulong const,reverse_t>::position = i.position;
    return *this;
      }
      // // Refines forward iterator
      // // operator++ (both)
      const_iterator_t<reverse_t>& operator++(){
    iterator_base<matrix_ulong const,vector_ulong const,reverse_t>::increment();
    return *this;
      }
      const_iterator_t<reverse_t> operator++( int ){
    // return value
    const_iterator_t<reverse_t> result( *this );
    iterator_base<matrix_ulong const,vector_ulong const,reverse_t>::increment();
    return result;
      }
      // // Refines bidirectional iterator
      // // operator-- (both)
      const_iterator_t<reverse_t>& operator--(){
    iterator_base<matrix_ulong const,vector_ulong const,reverse_t>::decrement();
    return *this;
      }
      const_iterator_t<reverse_t> operator--( int ){
    // return value
    const_iterator_t<reverse_t> result( *this );
    iterator_base<matrix_ulong const,vector_ulong const,reverse_t>::decrement();
    return result;
      }
      const_iterator_t() : iterator_base<matrix_ulong,vector_ulong,reverse_t>(){}
      const_iterator_t( iterator_t<reverse_t> const&i ){
    const_iterator_t<reverse_t>::v = i.v;
    const_iterator_t<reverse_t>::position = i.position;
      }
      bool operator==( iterator_t<reverse_t> const& i ) const {
        return this->v == i.v and this->position == i.position;
      }
      bool operator!=( iterator_t<reverse_t> const& i ) const {
        return not this->operator==( i );
      }
      bool operator==( const_iterator_t<reverse_t> const& i ) const {
        return this->v == i.v and this->position == i.position;
      }
      bool operator!=( const_iterator_t<reverse_t> const& i ) const {
        return not this->operator==( i );
      }
    protected:
      // We need a constructor for matrix_ulong
      friend class matrix_ulong;
      const_iterator_t( matrix_ulong const* v, size_t position )
      : iterator_base<matrix_ulong const,vector_ulong const,reverse_t>( v, position ){}
    };
  public:
    typedef const_iterator_t<false> const_iterator;
    typedef iterator_t<false> iterator;
    typedef const_iterator_t<true> const_reverse_iterator;
    typedef iterator_t<true> reverse_iterator;
    typedef size_t size_type;
    // begin()
    iterator begin(){
      return iterator( this, 0 );
    }
    const_iterator begin() const {
      return const_iterator( this, 0 );
    }
    // end()
    iterator end(){
      if( ccgsl_pointer == 0 ) return iterator( this, 0 );
      return iterator( this, size1() );
    }
    const_iterator end() const {
      if( ccgsl_pointer == 0 ) return const_iterator( this, 0 );
      return const_iterator( this, size1() );
    }
    // rbegin()
    reverse_iterator rbegin(){
      if( ccgsl_pointer ==0 ) return reverse_iterator( this, 0 );
      return reverse_iterator( this, size1() - 1 );
    }
    const_reverse_iterator rbegin() const {
      if( ccgsl_pointer == 0 ) return const_reverse_iterator( this, 0 );
      return const_reverse_iterator( this, size1() - 1 );
    }
    // rend()
    reverse_iterator rend(){
      return reverse_iterator( this, -1 );
    }
    const_reverse_iterator rend() const {
      return const_reverse_iterator( this, -1 );
    }
  public:
    // Sizes
    size_t size1() const { return ccgsl_pointer == 0 ? 0 : ccgsl_pointer->size1; }
    size_t size2() const { return ccgsl_pointer == 0 ? 0 : ccgsl_pointer->size2; }
    unsigned long* data() { 
      if( ccgsl_pointer == 0 ) gsl_error( "null vector_ulong", __FILE__, __LINE__, GSL_EFAULT );
#ifndef GSL_RANGE_CHECK_OFF
      if( ccgsl_pointer->size2 != ccgsl_pointer->tda )
    gsl_error( "matrix_ulong size2 and tda do not match", __FILE__, __LINE__, GSL_EBADLEN );
#endif
      return ccgsl_pointer->data; }
    unsigned long const* data() const {
      if( ccgsl_pointer == 0 ) gsl_error( "null vector_ulong", __FILE__, __LINE__, GSL_EFAULT );
#ifndef GSL_RANGE_CHECK_OFF
      if( ccgsl_pointer->size2 != ccgsl_pointer->tda )
    gsl_error( "matrix_ulong size2 and tda do not match", __FILE__, __LINE__, GSL_EBADLEN );
#endif
      return ccgsl_pointer->data; }
    void swap( matrix_ulong& m ){
      std::swap( ccgsl_pointer, m.ccgsl_pointer );
      std::swap( count, m.count );
    }
    void tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, matrix_ulong const& src){
      gsl_matrix_ulong_tricpy( Uplo, Diag, get(), src.get() );
    }
    void transpose_tricpy(CBLAS_UPLO_t Uplo, CBLAS_DIAG_t Diag, matrix_ulong const& src){
      gsl_matrix_ulong_transpose_tricpy( Uplo, Diag, get(), src.get() );
    }
    // view operations
    matrix_ulong submatrix( size_t const i, size_t const j, size_t const n1, size_t const n2 ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_submatrix( get(), i, j, n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    vector_ulong row( size_t const i ){
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_row( get(), i ).vector;
      return vector_ulong( w );
    }
    vector_ulong column( size_t const j ){
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_column( get(), j ).vector;
      return vector_ulong( w );
    }
    vector_ulong
    diagonal(){ gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_diagonal( get() ).vector;
      return vector_ulong( w );
    }
    vector_ulong subdiagonal( size_t const k ){
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_subdiagonal( get(), k ).vector;
      return vector_ulong( w );
    }
    vector_ulong superdiagonal( size_t const k ){
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_superdiagonal( get(), k ).vector;
      return vector_ulong( w );
    }
    vector_ulong subrow( size_t const i, size_t const offset, size_t const n ){
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_subrow( get(), i, offset, n ).vector;
      return vector_ulong( w );
    }
    vector_ulong subcolumn( size_t const j, size_t const offset, size_t const n ){
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_subcolumn( get(), j, offset, n ).vector;
      return vector_ulong( w );
    }
    static matrix_ulong view_array( unsigned long* base, size_t const n1, size_t const n2 ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_view_array( base, n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong view_array_with_tda( unsigned long* base, size_t const n1, size_t const n2, size_t const tda ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m =  gsl_matrix_ulong_view_array_with_tda( base, n1, n2, tda ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong view_vector( vector_ulong& v, size_t const n1, size_t const n2 ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m =  gsl_matrix_ulong_view_vector( v.get(), n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong view_vector_with_tda( vector_ulong& v, size_t const n1, size_t const n2, size_t const tda ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_view_vector_with_tda( v.get(), n1, n2, tda ).matrix;
      return matrix_ulong( m );
    }
    // const versions ...
    matrix_ulong const const_submatrix( size_t const i, size_t const j, size_t const n1, size_t const n2 ) const {
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_const_submatrix( get(), i, j, n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    vector_ulong const const_row( size_t const i ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_row( get(), i ).vector;
      return vector_ulong( w );
    }
    vector_ulong const const_column( size_t const j ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_column( get(), j ).vector;
      return vector_ulong( w );
    }
    vector_ulong const const_diagonal() const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_diagonal( get() ).vector;
      return vector_ulong( w );
    }
    vector_ulong const const_subdiagonal( size_t const k ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_subdiagonal( get(), k ).vector;
      return vector_ulong( w );
    }
    vector_ulong const const_superdiagonal( size_t const k ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_superdiagonal( get(), k ).vector;
      return vector_ulong( w );
    }
    vector_ulong const const_subrow( size_t const i, size_t const offset, size_t const n ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_subrow( get(), i, offset, n ).vector;
      return vector_ulong( w );
    }
    vector_ulong const const_subcolumn( size_t const j, size_t const offset, size_t const n ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_subcolumn( get(), j, offset, n ).vector;
      return vector_ulong( w );
    }
    matrix_ulong const submatrix( size_t const i, size_t const j, size_t const n1, size_t const n2 ) const {
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_const_submatrix( get(), i, j, n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    vector_ulong const row( size_t const i ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_row( get(), i ).vector;
      return vector_ulong( w );
    }
    vector_ulong const column( size_t const j ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_column( get(), j ).vector;
      return vector_ulong( w );
    }
    vector_ulong const diagonal() const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_diagonal( get() ).vector;
      return vector_ulong( w );
    }
    vector_ulong const subdiagonal( size_t const k ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_subdiagonal( get(), k ).vector;
      return vector_ulong( w );
    }
    vector_ulong const superdiagonal( size_t const k ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_superdiagonal( get(), k ).vector;
      return vector_ulong( w );
    }
    vector_ulong const subrow( size_t const i, size_t const offset, size_t const n ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_subrow( get(), i, offset, n ).vector;
      return vector_ulong( w );
    }
    vector_ulong const subcolumn( size_t const j, size_t const offset, size_t const n ) const {
      gsl_vector_ulong* w = static_cast<gsl_vector_ulong*>( malloc( sizeof( gsl_vector_ulong ) ) );
      *w = gsl_matrix_ulong_const_subcolumn( get(), j, offset, n ).vector;
      return vector_ulong( w );
    }
    static matrix_ulong const const_view_array( unsigned long const* base, size_t const n1, size_t const n2 ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_const_view_array( base, n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong const
    const_view_array_with_tda( unsigned long const* base, size_t const n1, size_t const n2, size_t const tda ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m =  gsl_matrix_ulong_const_view_array_with_tda( base, n1, n2, tda ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong const const_view_vector( vector_ulong const& v, size_t const n1, size_t const n2 ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m =  gsl_matrix_ulong_const_view_vector( v.get(), n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong const
    const_view_vector_with_tda( vector_ulong const& v, size_t const n1, size_t const n2, size_t const tda ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_const_view_vector_with_tda( v.get(), n1, n2, tda ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong const view_array( unsigned long const* base, size_t const n1, size_t const n2 ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_const_view_array( base, n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong const
    view_array_with_tda( unsigned long const* base, size_t const n1, size_t const n2, size_t const tda ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m =  gsl_matrix_ulong_const_view_array_with_tda( base, n1, n2, tda ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong const view_vector( vector_ulong const& v, size_t const n1, size_t const n2 ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m =  gsl_matrix_ulong_const_view_vector( v.get(), n1, n2 ).matrix;
      return matrix_ulong( m );
    }
    static matrix_ulong const
    view_vector_with_tda( vector_ulong const& v, size_t const n1, size_t const n2, size_t const tda ){
      gsl_matrix_ulong* m = static_cast<gsl_matrix_ulong*>( malloc( sizeof( gsl_matrix_ulong ) ) );
      *m = gsl_matrix_ulong_const_view_vector_with_tda( v.get(), n1, n2, tda ).matrix;
      return matrix_ulong( m );
    }
  private:
    bool owns_data;
    gsl_matrix_ulong* ccgsl_pointer;
    size_t* count;
  public:
    // shared reference functions
    gsl_matrix_ulong* get(){ return ccgsl_pointer; }
    gsl_matrix_ulong const* get() const { return ccgsl_pointer; }
    bool unique() const { return count != 0 and *count == 1; }
    size_t use_count() const { return count == 0 ? 0 : *count; }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
    explicit
#endif
    operator bool() const { return ccgsl_pointer != 0; }
    // GSL functions
    static matrix_ulong calloc( size_t const n1, size_t const n2 ){ return matrix_ulong( gsl_matrix_ulong_calloc( n1, n2 ) ); }
    void set_zero(){ gsl_matrix_ulong_set_zero( get() ); }
    void set_all( unsigned long x ){ gsl_matrix_ulong_set_all( get(), x ); }
    int memcpy( matrix_ulong const& src ){ return gsl_matrix_ulong_memcpy( get(), src.get() ); }
    unsigned long max() const { return gsl_matrix_ulong_max( get() ); }
    unsigned long min() const { return gsl_matrix_ulong_min( get() ); }
    void minmax( unsigned long* min_out, unsigned long* max_out ) const {
      gsl_matrix_ulong_minmax( get(), min_out, max_out ); }
    void minmax( unsigned long& min_out, unsigned long& max_out ) const {
      gsl_matrix_ulong_minmax( get(), &min_out, &max_out ); }
    int add( matrix_ulong const& b ){ return gsl_matrix_ulong_add( get(), b.get() ); }
    int sub( matrix_ulong const& b ){ return gsl_matrix_ulong_sub( get(), b.get() ); }
    int scale( unsigned long const x ){ return gsl_matrix_ulong_scale( get(), x ); }
    int add_constant( unsigned long const x ){ return gsl_matrix_ulong_add_constant( get(), x ); }
    int isnull() const { return gsl_matrix_ulong_isnull( get() ); }
    int ispos() const { return gsl_matrix_ulong_ispos( get() ); }
    int isneg() const { return gsl_matrix_ulong_isneg( get() ); }
    int isnonneg() const { return gsl_matrix_ulong_isnonneg( get() ); }
    unsigned long get( size_t const i, size_t const j ) const { return gsl_matrix_ulong_get( get(), i, j ); }
    void set( size_t const i, size_t const j, unsigned long x ){ gsl_matrix_ulong_set( get(), i, j, x ); }
    unsigned long* ptr( size_t const i, size_t const j ){ return gsl_matrix_ulong_ptr( get(), i, j ); }
    unsigned long const* const_ptr( size_t const i, size_t const j ) const {
      return gsl_matrix_ulong_const_ptr( get(), i, j ); }
    int fread( FILE* stream ){ return gsl_matrix_ulong_fread( stream, get() ); }
    int fwrite( FILE* stream ) const { return gsl_matrix_ulong_fwrite( stream, get() ); }
    int fscanf( FILE* stream ){ return gsl_matrix_ulong_fscanf( stream, get() ); }
    int fprintf( FILE* stream, char const* format ) const {
      return gsl_matrix_ulong_fprintf( stream, get(), format ); }
    matrix_ulong( block_ulong& b, size_t const offset, size_t const n1, size_t const n2, size_t const d2 ){
      ccgsl_pointer = gsl_matrix_ulong_alloc_from_block( b.get(), offset, n1, n2, d2 );
      // just plausibly we could allocate vector_ulong but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
        // try to tidy up before rethrowing
        gsl_matrix_ulong_free( ccgsl_pointer );
        throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
    }
    matrix_ulong( matrix_ulong& m, size_t const k1, size_t const k2, size_t const n1, size_t const n2 ){
      ccgsl_pointer = gsl_matrix_ulong_alloc_from_matrix( m.get(), k1, k2, n1, n2 );
      // just plausibly we could allocate matrix_ulong but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
        // try to tidy up before rethrowing
        gsl_matrix_ulong_free( ccgsl_pointer );
        throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
    }
    // More functions
    void set_identity(){ gsl_matrix_ulong_set_identity( get() ); }
    int swap_rows( size_t const i, size_t const j ){ return gsl_matrix_ulong_swap_rows( get(), i, j ); }
    int swap_columns( size_t const i, size_t const j ){
      return gsl_matrix_ulong_swap_columns( get(), i, j ); }
    int swap_rowcol( size_t const i, size_t const j ){ return gsl_matrix_ulong_swap_rowcol( get(), i, j ); }
    int transpose(){ return gsl_matrix_ulong_transpose( get() ); }
    int transpose_memcpy( matrix_ulong const& src ){
      return gsl_matrix_ulong_transpose_memcpy( get(), src.get() ); }
#ifndef DOXYGEN_SKIP
    void max_index( size_t* imax, size_t* jmax ) const {
      gsl_matrix_ulong_max_index( get(), imax, jmax ); }
    void min_index( size_t* imin, size_t* jmin ) const {
      gsl_matrix_ulong_min_index( get(), imin, jmin ); }
    void minmax_index( size_t* imin, size_t* jmin, size_t* imax, size_t* jmax ) const {
      gsl_matrix_ulong_minmax_index( get(), imin, jmin, imax, jmax ); }
#endif
    void max_index( size_t& imax, size_t& jmax ) const {
      gsl_matrix_ulong_max_index( get(), &imax, &jmax ); }
    void min_index( size_t& imin, size_t& jmin ) const {
      gsl_matrix_ulong_min_index( get(), &imin, &jmin ); }
    void minmax_index( size_t& imin, size_t& jmin, size_t& imax, size_t& jmax ) const {
      gsl_matrix_ulong_minmax_index( get(), &imin, &jmin, &imax, &jmax ); }
    int
    mul_elements( matrix_ulong const& b ){
      return gsl_matrix_ulong_mul_elements( get(), b.get() ); }
    int div_elements( matrix_ulong const& b ){
      return gsl_matrix_ulong_div_elements( get(), b.get() ); }
    unsigned long norm1() const {
      return gsl_matrix_ulong_norm1( get() ); }
    int scale_rows( vector_ulong const& x ){
      return gsl_matrix_ulong_scale_rows( get(), x.get() ); }
    int scale_columns( vector_ulong const& x ){
      return gsl_matrix_ulong_scale_columns( get(), x.get() ); }
    int add_diagonal( unsigned long const x ){
      return gsl_matrix_ulong_add_diagonal( get(), x ); }
    int get_row( vector_ulong& v, size_t const i ) const {
      return gsl_matrix_ulong_get_row( v.get(), get(), i ); }
    int get_col( vector_ulong& v, size_t const j ) const {
      return gsl_matrix_ulong_get_col( v.get(), get(), j ); }
    int set_row( size_t const i, vector_ulong const& v ){
      return gsl_matrix_ulong_set_row( get(), i, v.get() ); }
    int set_col( size_t const j, vector_ulong const& v ){
      return gsl_matrix_ulong_set_col( get(), j, v.get() ); }
    // Extra functions for []
    vector_ulong operator[]( size_t const i ){
      // First check that iterator is initialised.
      if( ccgsl_pointer == 0 ){
    gsl_error( "matrix_ulong is null", __FILE__, __LINE__, exception::GSL_EFAULT );
    return vector_ulong();
      }
#ifndef GSL_RANGE_CHECK_OFF
#ifndef __GXX_EXPERIMENTAL_CXX0X__
      static vector_ulong something;
#endif
      //Check that position make sense
      if( i >= size1() ){
    gsl_error( "trying to read beyond last row of matrix_ulong",
           __FILE__, __LINE__, exception::GSL_EINVAL );
#ifdef __GXX_EXPERIMENTAL_CXX0X__
    return vector_ulong();
#else
    return something;
#endif
      }
      // n is a valid position
#endif
      vector_ulong v(0);
      gsl_vector_ulong_view w = gsl_matrix_ulong_row( ccgsl_pointer, i );
      v.wrap_gsl_vector_ulong_without_ownership(&(w.vector));
      return v;
    }
    vector_ulong const operator[]( size_t const i ) const {
      // First check that iterator is initialised.
      if( ccgsl_pointer == 0 ){
    gsl_error( "matrix_ulong is null", __FILE__, __LINE__, exception::GSL_EFAULT );
    return vector_ulong();
      }
      vector_ulong v(0);
      gsl_vector_ulong_view w = gsl_matrix_ulong_row( ccgsl_pointer, i );
      v.wrap_gsl_vector_ulong_without_ownership(&(w.vector));
      return v;
    }
    inline int
    permute( permutation& p ){
      return gsl_permute_matrix_ulong( p.get(), get() ); } 
  };
 
  // Extra functions for vector_ulong allocation from matrix_ulong objects
  inline vector_ulong vector_ulong::alloc_row_from_matrix( matrix_ulong& m, size_t const i ){
    return vector_ulong ( gsl_vector_ulong_alloc_row_from_matrix( m.get(), i ) ); }
  inline vector_ulong vector_ulong::alloc_col_from_matrix( matrix_ulong& m, size_t const i ){
    return vector_ulong ( gsl_vector_ulong_alloc_col_from_matrix( m.get(), i ) ); }
}
#endif