vector_complex.hpp

Go to the documentation of this file.

/*
 * $Id: vector_complex.hpp 314 2014-02-22 14:31:16Z jdl3 $
 * Copyright (C) 2010, 2011, 2012, 2019, 2020, 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_VECTOR_COMPLEX_HPP
#define CCGSL_VECTOR_COMPLEX_HPP
 
#include<gsl/gsl_vector_complex_double.h>
#include<new>
#include<iterator>
#ifdef __GXX_EXPERIMENTAL_CXX0X__
#include<complex>
#endif
 
#include"exception.hpp"
#include"block_complex.hpp"
#include"vector.hpp"
 
// This file is a template
#define CCGSL_MTY 2
 
namespace gsl {
  // declare matrix_complex class
  class matrix_complex;
  class vector_complex {
  public:
    vector_complex(){
      owns_data = false;
      ccgsl_pointer = 0;
      count = 0; // initially nullptr will do
    }
    // Refines random access container
    // Refines assignable
    explicit vector_complex( size_t const n )
      : owns_data(true) {
      if( n > 0 ) ccgsl_pointer = gsl_vector_complex_alloc( n );
      else { ccgsl_pointer = new gsl_vector_complex; ccgsl_pointer->size = 0; ccgsl_pointer->data = 0; }
      // just plausibly we could allocate vector_complex but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
    // try to tidy up before rethrowing
    if( n > 0 ) gsl_vector_complex_free( ccgsl_pointer );
    else delete ccgsl_pointer;
    throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
    }
#ifndef DOXYGEN_SKIP
    explicit vector_complex( int const n )
      : owns_data(true) {
      if( n > 0 ) ccgsl_pointer = gsl_vector_complex_alloc( static_cast<size_t>( n ) );
      else if(n<0)
    gsl_error("failed tring to make a vector of negative length",
          __FILE__, __LINE__, exception::GSL_EDOM );
      else { ccgsl_pointer = new gsl_vector_complex; ccgsl_pointer->size = 0; ccgsl_pointer->data = 0; }
      // just plausibly we could allocate vector_complex but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
    // try to tidy up before rethrowing
    if( n > 0 ) gsl_vector_complex_free( ccgsl_pointer );
    else delete ccgsl_pointer;
    throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
    }
#endif // DOXYGEN_SKIP
    explicit vector_complex( gsl_vector_complex* 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__
    vector_complex( std::initializer_list<std::complex<double> > initializer_list )
      : owns_data{true}{
      size_t const n = initializer_list.size();
      ccgsl_pointer = gsl_vector_complex_alloc( n );
      // just plausibly we could allocate vector but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
    // try to tidy up before rethrowing
    if( n > 0 ) gsl_vector_complex_free( ccgsl_pointer );
    else delete ccgsl_pointer;
    throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
      // now copy
      auto p = begin();
      for( auto x : initializer_list ){
    *p = complex( std::move( x ) ); ++p;
      }
    }
#endif
    // copy constructor
    vector_complex( vector_complex const& v )
      : owns_data(v.owns_data),
    ccgsl_pointer( v.ccgsl_pointer ), count( v.count ){
      if( count != 0 ) ++*count; // vector_complex is now shared.
    }
    // assignment operator
    vector_complex& operator=( vector_complex const& v ){
      // first, possibly delete anything pointed to by @c this
      if( count == 0 or --*count == 0 ){
    if( ccgsl_pointer != 0 ){
      if( ccgsl_pointer->size > 0 ) gsl_vector_complex_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; // vector_complex is now shared.
      return *this;
    }
    // clone()
    vector_complex clone() const {
      vector_complex copy( get()->size );
      // Now copy
      gsl_vector_complex_memcpy( copy.get(), get() );
      // return new object
      return copy;
    }
    // destructor
    ~vector_complex(){
      if( count != 0 and --*count == 0 ){
    // could have allocated null pointer
    if( ccgsl_pointer != 0 and owns_data ){
      if( ccgsl_pointer->size > 0 ) gsl_vector_complex_free( ccgsl_pointer );
      else delete ccgsl_pointer; }
    delete count;
      }
    }
    // Allow possibility of assigning from gsl_vector_complex without sharing
    void wrap_gsl_vector_complex_without_ownership( gsl_vector_complex* v ){
      if( count != 0 and --*count == 0 ){
    // could have allocated null pointer
    if( ccgsl_pointer != 0 ){
      if( ccgsl_pointer->size != 0 ) gsl_vector_complex_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
    }
    // Refines equality comparable
    // == operator
    bool operator==( vector_complex const& v ) const {
      // trivially equal if gsl_vector_complex*s are identical
      if( ccgsl_pointer == v.ccgsl_pointer ) return true;
      // trivially not equal if one is zero: != should be same as xor here
      if( (ccgsl_pointer == 0) != (v.ccgsl_pointer == 0 ) ) return false;
      // trivially not equal if sizes are different
      if( ccgsl_pointer->size != v.ccgsl_pointer->size ) return false;
      // check elementwise for equality
      for( size_t i = 0; i < CCGSL_MTY * ccgsl_pointer->size; ++i )
    if( *(ccgsl_pointer->data + i) != *(v.ccgsl_pointer->data + i) )
      return false; 
      return true;
    }
    void reset(){ vector_complex().swap( *this ); }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
    vector_complex( vector_complex&& v )
      : owns_data (v.owns_data),
    ccgsl_pointer( v.ccgsl_pointer ), count( nullptr ){
      std::swap( count, v.count );
      v.ccgsl_pointer = nullptr;
    }
    vector_complex& operator=( vector_complex&& v ){
      vector_complex( std::move( v ) ).swap( *this );
      return *this;
    }
#endif
    // != operator
    bool operator!=( vector_complex const& v ) const { return not operator==( v ); }
    // Refines forward container
    // Refines less than comparable
    // operator<
    bool operator<( vector_complex const& v ) const {
      // null vector_complex comes first
      if( ccgsl_pointer == 0 ) return v.ccgsl_pointer != 0;
      // if v is null then this > v
      if( v.ccgsl_pointer == 0 ) return false;
      // Now compare elementwise
      size_t const  size = ccgsl_pointer->size;
      size_t const  v_size = v.ccgsl_pointer->size;
      size_t const min = size > v_size ? size : v_size;
      for( size_t i = 0; i < min; ++i ){
    complex const t = gsl_vector_complex_get( ccgsl_pointer, i );
    complex const u =gsl_vector_complex_get( v.ccgsl_pointer, i );
    if( t < u ) return true;
    if( u < t ) return false;
      }
      // elements match.
      return size < v_size;
    }
    // operator>
    bool operator>( vector_complex const& v ) const {
      // null vector_complex comes first
      if( ccgsl_pointer == 0 ) return false;
      // if v is null then this > v
      if( v.ccgsl_pointer == 0 ) return true;
      // Now compare elementwise
      size_t const  size = ccgsl_pointer->size;
      size_t const  v_size = v.ccgsl_pointer->size;
      size_t const min = size > v_size ? size : v_size;
      for( size_t i = 0; i < min; ++i ){
    complex const t = gsl_vector_complex_get( ccgsl_pointer, i );
    complex const u =gsl_vector_complex_get( v.ccgsl_pointer, i );
    if( t > u ) return true;
    if( u > t ) return false;
      }
      // elements match.
      return size > v_size;
    }
    // operator<=
    bool operator<=( vector_complex const& v ) const {
      return operator<( v ) or operator==( v );
    }
    // operator>=
    bool operator>=( vector_complex const& v ) const {
      return operator>( v ) or operator==( v );
    }
    // Refines container
    // type value_type
    typedef complex value_type;
    // type reference
    typedef  complex_ref reference;
    // type const_reference
    typedef  reference const const_reference;
    // type pointer
    typedef  complex_ptr pointer;
    // type const_pointer
    typedef  complex_ptr const const_pointer;
    // type iterator
  private:
    template<typename container, typename content,bool reverse_t>
    class iterator_base {
      friend class vector_complex;
    public:
      typedef complex value_type;
      typedef complex_ref reference;
      typedef complex_ptr pointer;
      typedef std::random_access_iterator_tag iterator_category;
      // // type iterator_traits<vector_complex>::difference_type
      typedef ptrdiff_t difference_type;
    public:
      // // operator*
      reference operator*() const {
    // First check that iterator is initialised.
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return complex_ref( 0 );
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return complex_ref( 0 );
    }
    // Check that position make sense
    if( position >= static_cast<difference_type>( v->size() ) ){
      gsl_error( "trying to dereference beyond rbegin()", __FILE__, __LINE__, exception::GSL_EFAILED );
      return complex_ref( 0 );
    }
    if( position <= -1 ){
      gsl_error( "trying to dereference beyond begin()", __FILE__, __LINE__, exception::GSL_EFAILED );
      return complex_ref( 0 );
    }
    // position and v are valid: return data
    return complex_ref( v->ccgsl_pointer->data + CCGSL_MTY * position * v->ccgsl_pointer->stride );
      }
      // // operator->
      pointer operator->() const {
    // First check that iterator is initialised.
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return pointer( 0 );
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return pointer( 0 );
    }
    // Check that position make sense
    if( position >= static_cast<difference_type>( v->size() ) ){
      gsl_error( "trying to dereference end()", __FILE__, __LINE__, exception::GSL_EFAILED );
      return pointer( 0 );
    }
    if( position <= -1 ){
      gsl_error( "trying to dereference rend()", __FILE__, __LINE__, exception::GSL_EFAILED );
      return pointer( 0 );
    }
    // position and v are valid: return data
    return pointer( v->ccgsl_pointer->data + CCGSL_MTY * position * v->ccgsl_pointer->stride );
      }
      // // operator[]
      reference operator[]( difference_type const n ) const {
    // First check that iterator is initialised.
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return reference( 0 );
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return reference( 0 );
    }
    // Check that position make sense
    difference_type const p = reverse_t ? position - n : position + n;
    if( p >= static_cast<difference_type>( v->size() ) ){
      gsl_error( "trying to dereference beyond rbegin()", __FILE__, __LINE__, exception::GSL_EFAILED );
      return reference( 0 );
    }
    if( p <= -1 ){
      gsl_error( "trying to dereference beyond begin()", __FILE__, __LINE__, exception::GSL_EFAILED );
      return reference( 0 );
    }
    // p is a valid position
    return reference( v->ccgsl_pointer->data + CCGSL_MTY * p * v->ccgsl_pointer->stride );
      }
      // // operator-: find distance between two iterators
      difference_type operator-( iterator_base<container,content,reverse_t> const& i ) const {
    // Warn if either iterator undefined
    if( v == 0 or i.v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return 0;
    } else if( v->ccgsl_pointer == 0 or i.v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return 0;
    }
    // Warn if iterators do not point to same vector_complex
    if( v->ccgsl_pointer != i.v->ccgsl_pointer ){
      gsl_error( "trying to take difference of iterators for different vector_complex objects", __FILE__, __LINE__,
             exception::GSL_EFAILED );
      return 0;
    }
    return reverse_t ? i.position - position : position - i.position;
      }
      // // operator!=
      // // operator<
      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 );
      }
      bool operator<( iterator_base<container,content,reverse_t> const& i ) const {
    // Warn if either iterator undefined
    if( v == 0 or i.v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return false;
    }
    // Warn if iterators do not point to same vector_complex
    if( v->ccgsl_pointer != i.v->ccgsl_pointer ){
      gsl_error( "trying to take difference of iterators for different vector_complex objects", __FILE__, __LINE__,
             exception::GSL_EFAILED );
      return false;
    }
    return reverse_t ? i.position < position : position < i.position;
      }
    protected:
      void increment(){
    // Only makes sense if v points to a vector_complex
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    }
    // increment position and check against size
    if( reverse_t ){ if( position >= 0 ) --position; }
    else { if( position < static_cast<difference_type>( v->size() ) ) ++position; }
      }
      void decrement(){
    // Only makes sense if v points to a vector_complex
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    }
    // decrement position and check against size
    if( reverse_t ){ if( position < static_cast<difference_type>( v->size() ) ) ++position; }
    else { if( position >= 0 ) --position; }
      }
      void shift( difference_type const n ){
    // Warn if iterator undefined
    if( v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    } else if( v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return;
    }
    position += reverse_t ? -n : n;
      }
      iterator_base(){ v = 0; }
      iterator_base( container* v, difference_type position )
    : v( v ), position( position ){}
      container* v;
      difference_type 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<vector_complex,double,reverse_t>{
    public:
      // // Refines output iterator
      // // operator=
      iterator_t<reverse_t>& operator=( iterator_t<reverse_t> const& i ){
    iterator_base<vector_complex,double,reverse_t>::v = i.v;
    iterator_base<vector_complex,double,reverse_t>::position = i.position;
    return *this;
      }
      // // Refines forward iterator
      // // operator++ (both)
      iterator_t<reverse_t>& operator++(){
    iterator_base<vector_complex,double,reverse_t>::increment();
    return *this;
      }
      iterator_t<reverse_t> operator++( int ){
    // return value
    iterator_t<reverse_t> result( *this );
    iterator_base<vector_complex,double,reverse_t>::increment();
    return result;
      }
      // // Refines bidirectional iterator
      // // operator-- (both)
      iterator_t<reverse_t>& operator--(){
    iterator_base<vector_complex,double,reverse_t>::decrement();
    return *this;
      }
      iterator_t<reverse_t> operator--( int ){
    // return value
    iterator_t<reverse_t> result( *this );
    iterator_base<vector_complex,double,reverse_t>::decrement();
    return result;
      }
      typedef typename iterator_base<vector_complex,double,reverse_t>::difference_type
    difference_type;
      // // operator+=
      iterator_t<reverse_t>& operator+=( difference_type const n ){
    this->shift( n );
    return *this;
      }
      // // operator-=
      iterator_t<reverse_t>& operator-=( difference_type const n ){
    this->shift( -n );
    return *this;
      }
      // // operator+ (n+i)(i+n)
      iterator_t<reverse_t> operator+( difference_type const n ) const {
    iterator_t<reverse_t> result( *this );
    result.shift( n );
    return result;
      }
      // // operator- (n-i)(i-n)(i-j)
      iterator_t<reverse_t> operator-( difference_type const n ) const {
    iterator_t<reverse_t> result( *this );
    result.shift( -n );
    return result;
      }
      difference_type operator-( iterator_t<reverse_t> const& i ) const {
    return iterator_base<vector_complex,double,reverse_t>::operator-( i );
      }
      difference_type operator-( const_iterator_t<reverse_t> const& i ) const {
    // Warn if either iterator undefined
    if( this->v == 0 or i.v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return 0;
    } else if( this->v->ccgsl_pointer == 0 or i.v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return 0;
    }
    // Warn if iterators do not point to same vector_complex
    if( this->v->ccgsl_pointer != i.v->ccgsl_pointer ){
      gsl_error( "trying to take difference of iterators for different vector_complex objects", __FILE__, __LINE__,
             exception::GSL_EFAILED );
      return 0;
    }
    return reverse_t ? i.position - this->position : this->position - i.position;
      }
      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 );
      }
      bool operator<( const_iterator_t<reverse_t> const& i ) const {
    // Warn if either iterator undefined
    if( this->v == 0 or i.v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return false;
    }
    // Warn if iterators do not point to same vector_complex
    if( this->v->ccgsl_pointer != i.v->ccgsl_pointer ){
      gsl_error( "trying to take difference of iterators for different vector_complex objects", __FILE__, __LINE__,
             exception::GSL_EFAILED );
      return false;
    }
    return reverse_t ? i.position < this->position : this->position < i.position;
      }
      iterator_t() : iterator_base<vector_complex,double,reverse_t>(){}
    protected:
      friend class vector_complex;
      // We need a constructor for vector_complex
      iterator_t( vector_complex* v, difference_type position )
      : iterator_base<vector_complex,double,reverse_t>( v, position ){}
    };
    template<bool reverse_t> class const_iterator_t
      : public iterator_base<vector_complex const,double,reverse_t>{
    public:
      // // Refines output iterator
      // // operator=
      const_iterator_t<reverse_t>& operator=( const_iterator_t<reverse_t> const& i ){
    iterator_base<vector_complex const,double,reverse_t>::v = i.v;
    iterator_base<vector_complex const,double,reverse_t>::position = i.position;
    return *this;
      }
      // // Refines forward iterator
      // // operator++ (both)
      const_iterator_t<reverse_t>& operator++(){
    iterator_base<vector_complex const,double,reverse_t>::increment();
    return *this;
      }
      const_iterator_t<reverse_t> operator++( int ){
    // return value
    const_iterator_t<reverse_t> result( *this );
    iterator_base<vector_complex const,double,reverse_t>::increment();
    return result;
      }
      // // Refines bidirectional iterator
      // // operator-- (both)
      const_iterator_t<reverse_t>& operator--(){
    iterator_base<vector_complex const,double,reverse_t>::decrement();
    return *this;
      }
      const_iterator_t<reverse_t> operator--( int ){
    // return value
    const_iterator_t<reverse_t> result( *this );
    iterator_base<vector_complex const,double,reverse_t>::decrement();
    return result;
      }
      typedef typename iterator_base<vector_complex const,double,reverse_t>::difference_type
    difference_type;
      // // operator+=
      const_iterator_t<reverse_t>& operator+=( difference_type const n ){
    this->shift( n );
    return *this;
      }
      // // operator-=
      const_iterator_t<reverse_t>& operator-=( difference_type const n ){
    this->shift( -n );
    return *this;
      }
      // // operator+ (n+i)(i+n)
      const_iterator_t<reverse_t> operator+( difference_type const n ) const {
    const_iterator_t<reverse_t> result( *this );
    result += n;
    return result;
      }
      // // operator- (n-i)(i-n)(i-j)
      const_iterator_t<reverse_t> operator-( difference_type const n ) const {
    const_iterator_t<reverse_t> result( *this );
    result -= n;
    return result;
      }
      difference_type operator-( const_iterator_t<reverse_t> const& i ) const {
    return iterator_base<vector_complex const,double,reverse_t>::operator-( i );
      }
      difference_type operator-( iterator_t<reverse_t> const& i ) const {
    // Warn if either iterator undefined
    if( this->v == 0 or i.v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return 0;
    } else if( this->v->ccgsl_pointer == 0 or i.v->ccgsl_pointer == 0 ){
      gsl_error( "vector_complex not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return 0;
    }
    // Warn if iterators do not point to same vector_complex
    if( this->v->ccgsl_pointer != i.v->ccgsl_pointer ){
      gsl_error( "trying to take difference of iterators for different vector_complex objects", __FILE__, __LINE__,
             exception::GSL_EFAILED );
      return 0;
    }
    return reverse_t ? i.position - this->position : this->position - i.position;
      }
      const_iterator_t() : iterator_base<vector_complex,double,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<( iterator_t<reverse_t> const& i ) const {
    // Warn if either iterator undefined
    if( this->v == 0 or i.v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return false;
    }
    // Warn if iterators do not point to same vector_complex
    if( this->v->ccgsl_pointer != i.v->ccgsl_pointer ){
      gsl_error( "trying to take difference of iterators for different vector_complex objects", __FILE__, __LINE__,
             exception::GSL_EFAILED );
      return false;
    }
    return reverse_t ? i.position < this->position : this->position < i.position;
      }
      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 );
      }
      bool operator<( const_iterator_t<reverse_t> const& i ) const {
    // Warn if either iterator undefined
    if( this->v == 0 or i.v == 0 ){
      gsl_error( "iterator not initialised", __FILE__, __LINE__, exception::GSL_EFAILED );
      return false;
    }
    // Warn if iterators do not point to same vector_complex
    if( this->v->ccgsl_pointer != i.v->ccgsl_pointer ){
      gsl_error( "trying to take difference of iterators for different vector_complex objects", __FILE__, __LINE__,
             exception::GSL_EFAILED );
      return false;
    }
    return reverse_t ? i.position < this->position : this->position < i.position;
      }
    protected:
      // We need a constructor for vector_complex
      friend class vector_complex;
      const_iterator_t( vector_complex const* v, difference_type position )
      : iterator_base<vector_complex const,double,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;
    // type difference_type
    typedef  const_iterator::difference_type difference_type;
    // type size_type
    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, size() );
    }
    const_iterator end() const {
      if( ccgsl_pointer == 0 ) return const_iterator( this, 0 );
      return const_iterator( this, size() );
    }
    // size()
    size_type size() const { return ccgsl_pointer == 0 ? 0 : ccgsl_pointer->size; }
    double* data() { 
      if( ccgsl_pointer == 0 ) gsl_error( "null vector", __FILE__, __LINE__, GSL_EFAULT );
#ifndef GSL_RANGE_CHECK_OFF
      if( ccgsl_pointer->stride != 1 )
    gsl_error( "vector does not have stride of size 1", __FILE__, __LINE__, GSL_EBADLEN );
#endif
      return ccgsl_pointer->data; }
    double const* data() const {
      if( ccgsl_pointer == 0 ) gsl_error( "null vector", __FILE__, __LINE__, GSL_EFAULT );
#ifndef GSL_RANGE_CHECK_OFF
      if( ccgsl_pointer->stride != 1 )
    gsl_error( "vector does not have stride of size 1", __FILE__, __LINE__, GSL_EBADLEN );
#endif
      return ccgsl_pointer->data; }
    // max_size()
    size_type max_size() const { return ccgsl_pointer == 0 ? 0 : ccgsl_pointer->size; }
    // empty()
    bool empty() const { return ccgsl_pointer == 0 or ccgsl_pointer->size == 0; }
    // swap() --- should work even if sizes don't match
    void swap( vector_complex& v ){
      gsl_vector_complex* tmp = ccgsl_pointer; ccgsl_pointer = v.ccgsl_pointer; v.ccgsl_pointer = tmp;
      size_t* tmp2 = count; count = v.count; v.count = tmp2;
    }
    // Refines reversible container
    // rbegin()
    reverse_iterator rbegin(){
      if( ccgsl_pointer ==0 ) return reverse_iterator( this, 0 );
      return reverse_iterator( this, size() - 1 );
    }
    const_reverse_iterator rbegin() const {
      if( ccgsl_pointer ==0 ) return const_reverse_iterator( this, 0 );
      return const_reverse_iterator( this, size() - 1 );
    }
    // rend()
    reverse_iterator rend(){
      return reverse_iterator( this, -1 );
    }
    const_reverse_iterator rend() const {
      return const_reverse_iterator( this, -1 );
    }
    // operator[]
    complex_ref operator[]( size_t const n ){
      // First check that iterator is initialised.
      if( ccgsl_pointer == 0 ){
    gsl_error( "vector_complex is null", __FILE__, __LINE__, exception::GSL_EFAILED );
    return complex_ref( 0 );
      }
#ifndef GSL_RANGE_CHECK_OFF
      // Check that position make sense
      if( n >= size() ){
    gsl_error( "trying to read beyond end of vector_complex", __FILE__, __LINE__, exception::GSL_EFAILED );
    return complex_ref( 0 );
      }
      // n is a valid position
#endif
      return complex_ref( ccgsl_pointer->data + CCGSL_MTY * n * ccgsl_pointer->stride );
    }
    complex_ref const  operator[]( size_t const n ) const {
      // First check that iterator is initialised.
      if( ccgsl_pointer == 0 ){
    gsl_error( "vector_complex is null", __FILE__, __LINE__, exception::GSL_EFAILED );
    return complex_ref( 0 );
      }
#ifndef GSL_RANGE_CHECK_OFF
      // Check that position make sense
      if( n >= size() ){
    gsl_error( "trying to read beyond end of vector_complex", __FILE__, __LINE__, exception::GSL_EFAILED );
    return complex_ref( 0 );
      }
      // n is a valid position
#endif
      return complex_ref( ccgsl_pointer->data + n * ccgsl_pointer->stride );
    }
  private:
    bool owns_data;
    gsl_vector_complex* ccgsl_pointer;
    size_t* count;
  public:
    // shared reference functions
    gsl_vector_complex* get() { return ccgsl_pointer; }
    gsl_vector_complex 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 vector_complex calloc( size_t const n ){ return vector_complex( gsl_vector_complex_calloc( n ) ); }
    void set_zero(){ gsl_vector_complex_set_zero( get() ); }
    void set_all( complex x ){ gsl_vector_complex_set_all( get(), x ); }
    int set_basis( size_t i ){ return gsl_vector_complex_set_basis( get(), i ); }
    int memcpy( vector_complex const& src ){ return gsl_vector_complex_memcpy( get(), src.get() ); }
    int reverse(){ return gsl_vector_complex_reverse( get() ); }
    int swap_elements( size_t const i, size_t const j ){
      return gsl_vector_complex_swap_elements( get(), i, j ); }
    int add( vector_complex const& b ){ return gsl_vector_complex_add( get(), b.get() ); }
    int sub( vector_complex const& b ){ return gsl_vector_complex_sub( get(), b.get() ); }
    int mul( vector_complex const& b ){ return gsl_vector_complex_mul( get(), b.get() ); }
    int div( vector_complex const& b ){ return gsl_vector_complex_div( get(), b.get() ); }
    int scale( complex const x ){ return gsl_vector_complex_scale( get(), x ); }
    int add_constant( complex const x ){ return gsl_vector_complex_add_constant( get(), x ); }
    int axpby( complex const alpha, vector_complex const& x,
           complex const beta ){
      return gsl_vector_complex_axpby( alpha, x.get(), beta, get() );
    }
    int isnull() const { return gsl_vector_complex_isnull( get() ); }
    int ispos() const { return gsl_vector_complex_ispos( get() ); }
    int isneg() const { return gsl_vector_complex_isneg( get() ); }
    int isnonneg() const { return gsl_vector_complex_isnonneg( get() ); }
    complex get( size_t const i ) const { return gsl_vector_complex_get( get(), i ); }
    void set( size_t const i, complex x ){ gsl_vector_complex_set( get(), i, x ); }
    complex_ptr ptr( size_t const i ){
      if( i >= ccgsl_pointer->size )
    gsl_error( "Index out of range", __FILE__, __LINE__, exception::GSL_EINVAL );
      return complex_ptr( ccgsl_pointer->data + CCGSL_MTY * i ); }
    complex_ptr const const_ptr( size_t const i ){
      if( i >= ccgsl_pointer->size )
    gsl_error( "Index out of range", __FILE__, __LINE__, exception::GSL_EINVAL );
      return complex_ptr( ccgsl_pointer->data + CCGSL_MTY * i ); }
    int fread( FILE* stream ){ return gsl_vector_complex_fread( stream, get() ); }
    int fwrite( FILE* stream ) const { return gsl_vector_complex_fwrite( stream, get() ); }
    int fscanf( FILE* stream ){ return gsl_vector_complex_fscanf( stream, get() ); }
    int fprintf( FILE* stream, char const* format ) const {
      return gsl_vector_complex_fprintf( stream, get(), format ); }
    vector_complex( block_complex& b, size_t const offset, size_t const n, size_t const stride = 1 ){
      ccgsl_pointer = gsl_vector_complex_alloc_from_block( b.get(), offset, n, stride );
      // just plausibly we could allocate vector_complex but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
    // try to tidy up before rethrowing
    gsl_vector_complex_free( ccgsl_pointer );
    throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
    }
    vector_complex( vector_complex& v, size_t const offset, size_t const n, size_t const stride = 1 ){
      ccgsl_pointer = gsl_vector_complex_alloc_from_vector( v.get(), offset, n, stride );
      // just plausibly we could allocate vector_complex but not count
      try { count = new size_t; } catch( std::bad_alloc& e ){
    // try to tidy up before rethrowing
    gsl_vector_complex_free( ccgsl_pointer );
    throw e;
      }
      *count = 1; // initially there is just one reference to ccgsl_pointer
    }
    gsl::vector real(){
      return gsl::vector::view_array_with_stride( data(), CCGSL_MTY, size()  ); }
    gsl::vector const const_real() const {
      return gsl::vector::const_view_array_with_stride( data(), CCGSL_MTY, size()  ); }
    gsl::vector const real() const {
      return gsl::vector::const_view_array_with_stride( data(), CCGSL_MTY, size()  ); }
    gsl::vector imag(){
      return gsl::vector::view_array_with_stride( data() + CCGSL_MTY / 2, CCGSL_MTY, size()  ); }
    gsl::vector const const_imag() const {
      return gsl::vector::const_view_array_with_stride( data() + CCGSL_MTY / 2, CCGSL_MTY, size()  ); }
    gsl::vector const imag() const {
      return gsl::vector::const_view_array_with_stride( data() + CCGSL_MTY / 2, CCGSL_MTY, size()  ); }
    static vector_complex view_array( double* v, size_t n ){
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_view_array( v, n ).vector;
      return vector_complex( w );
    }
 
 
    static vector_complex view_array_with_stride( double* base, size_t stride, size_t n ){
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_view_array_with_stride( base, stride, n ).vector;
      return vector_complex( w );
    }
 
 
    static vector_complex const const_view_array( double const* v, size_t n ){
      gsl_vector_complex* w = static_cast<gsl_vector_complex *>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_view_array( v, n ).vector;
      return vector_complex( w );
    }
 
 
    static vector_complex const const_view_array_with_stride( double const* base, size_t stride, size_t n ){
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_view_array_with_stride( base, stride, n ).vector;
      return vector_complex( w );
    }
 
 
#ifndef DOXYGEN_SKIP
    static vector_complex const view_array( double const* v, size_t n ){
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_view_array( v, n ).vector;
      return vector_complex( w );
    }
#endif //DOXYGEN_SKIP
 
 
#ifndef DOXYGEN_SKIP
    static vector_complex const view_array_with_stride( double const* base, size_t stride, size_t n ){
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_view_array_with_stride( base, stride, n ).vector;
      return vector_complex( w );
    }  
#endif //DOXYGEN_SKIP
 
 
    vector_complex subvector( size_t i, size_t n ){
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_subvector( get(), i, n ).vector;
      return vector_complex( w );
    }
 
 
    vector_complex subvector_with_stride( size_t i, size_t stride, size_t n ){
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_subvector_with_stride( get(), i, stride, n ).vector;
      return vector_complex( w );
    }
 
 
    vector_complex const const_subvector( size_t i, size_t n ) const {
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_subvector( get(), i, n ).vector;
      return vector_complex( w );
    }
 
 
    vector_complex const const_subvector_with_stride( size_t i, size_t stride, size_t n ) const {
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_subvector_with_stride( get(), i, stride, n ).vector;
      return vector_complex( w );
    }
 
 
#ifndef DOXYGEN_SKIP
    vector_complex const subvector( size_t i, size_t n ) const {
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_subvector( get(), i, n ).vector;
      return vector_complex( w );
    }
#endif //DOXYGEN_SKIP
 
 
#ifndef DOXYGEN_SKIP
    vector_complex const subvector_with_stride( size_t i, size_t stride, size_t n ) const {
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_subvector_with_stride( get(), i, stride, n ).vector;
      return vector_complex( w );
    }
#endif //DOXYGEN_SKIP
 
 
    template<typename ARRAY>
    static vector_complex view_array( ARRAY& v, size_t n = 0 ){
      if(n > v.size())
        gsl_error( "n is too big", __FILE__, __LINE__, exception::GSL_EBADLEN );
      if(0 == n)
        n = v.size();
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_view_array( v.data(), n ).vector;
      return vector_complex( w );
    }
 
 
    template<typename ARRAY>
    static vector_complex view_array_with_stride( ARRAY& base, size_t stride, size_t n=0 ){
      if(0 == n)
        n = base.size()/stride;
      if((n-1)*stride > base.size())
        gsl_error( "n is too big", __FILE__, __LINE__, exception::GSL_EBADLEN );
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_view_array_with_stride( base.data(), stride, n ).vector;
      return vector_complex( w );
    }
 
 
    template<typename ARRAY>
    static vector_complex const const_view_array( ARRAY const& v, size_t n=0 ){
      if(n > v.size())
        gsl_error( "n is too big", __FILE__, __LINE__, exception::GSL_EBADLEN );
      if(0 == n)
        n = v.size();
      gsl_vector_complex* w = static_cast<gsl_vector_complex *>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_view_array( v.data(), n ).vector;
      return vector_complex( w );
    }
 
 
    template<typename ARRAY>
    static vector_complex const const_view_array_with_stride( ARRAY const& base, size_t stride, size_t n=0 ){
      if(0 == n)
        n = base.size()/stride;
      if((n-1)*stride > base.size())
        gsl_error( "n is too big", __FILE__, __LINE__, exception::GSL_EBADLEN );
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_view_array_with_stride( base.data(), stride, n ).vector;
      return vector_complex( w );
    }
 
 
#ifndef DOXYGEN_SKIP
    template<typename ARRAY>
    static vector_complex const view_array( ARRAY const& v, size_t n=0 ){
      if(n > v.size())
        gsl_error( "n is too big", __FILE__, __LINE__, exception::GSL_EBADLEN );
      if(0 == n)
        n = v.size();
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_view_array( v.data(), n ).vector;
      return vector_complex( w );
    }
#endif //DOXYGEN_SKIP
 
 
#ifndef DOXYGEN_SKIP
    template<typename ARRAY>
    static vector_complex const view_array_with_stride( ARRAY const& base, size_t stride, size_t n ){
      if(0 == n)
        n = base.size()/stride;
      if((n-1)*stride > base.size())
        gsl_error( "n is too big", __FILE__, __LINE__, exception::GSL_EBADLEN );
      gsl_vector_complex* w = static_cast<gsl_vector_complex*>( malloc( sizeof( gsl_vector_complex ) ) );
      *w = gsl_vector_complex_const_view_array_with_stride( base.data(), stride, n ).vector;
      return vector_complex( w );
    }  
#endif //DOXYGEN_SKIP
 
    // Extra allocators from matrix_complex objects. The definition must come after matrix_complex.
    static vector_complex alloc_row_from_matrix( matrix_complex& m, size_t const i );
    static vector_complex alloc_col_from_matrix( matrix_complex& m, size_t const j );
  };
 
  inline vector_complex::iterator operator+
  ( vector_complex::iterator::difference_type const n, vector_complex::iterator const& i ){ return i + n; }
  
  inline vector_complex::const_iterator operator+
  ( vector_complex::const_iterator::difference_type const n, vector_complex::const_iterator const& i ){ return i + n; }
  
  inline vector_complex::reverse_iterator operator+
  ( vector_complex::reverse_iterator::difference_type const n, vector_complex::reverse_iterator const& i ){
    return i + n; }
  
  inline vector_complex::const_reverse_iterator operator+
  ( vector_complex::const_reverse_iterator::difference_type const n, vector_complex::const_reverse_iterator const& i ){
    return i + n; }
  
}
#undef CCGSL_MTY
#endif