odeiv2_system.hpp

Go to the documentation of this file.

/*
 * $Id: odeiv2_system.hpp 289 2012-12-04 15:32:17Z jdl3 $
 * Copyright (C) 2012, 2016, 2021 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_ODEIV2_SYSTEM_HPP
#define CCGSL_ODEIV2_SYSTEM_HPP
 
#include"exception.hpp"
#include"vector.hpp"
#include"matrix.hpp"
#include<gsl/gsl_odeiv2.h>
 
namespace gsl {
  namespace odeiv2 {
    class system;
    template<typename T>
    void system_constructor( system& sys, T& t );
    class system : public gsl_odeiv2_system {
    public:
      struct Concept {
    template<typename T>
    int function( double const t, T const& y, T& dydt );
    template<typename T,typename U>
    int jacobian( double const t, T const&, U& dfdy, T& dfdt ){
      return exception::GSL_SUCCESS; }
    virtual size_t size() const = 0;
      };
#ifndef DOXYGEN_SKIP
    private:
      struct base_F {
    typedef int (*function_t)(double const,double const*,double*,void*);
    typedef int (*jacobian_t)(double const,double const*,double*,double*,void*);
    virtual ~base_F(){};
    virtual bool has_jacobian() const = 0;
      };
      template<typename T>
      struct Function : public base_F {
      public:
    Function( T& t ) : t( t ){
      size_t const size = t.size();
      vy.size = size;
      vy.stride = 1;
      vy.data = 0;
      vy.block = 0;
      vy.owner = 0;
      yv.wrap_gsl_vector_without_ownership( &vy );
      vdydt.size = size;
      vdydt.stride = 1;
      vdydt.data = 0;
      vdydt.block = 0;
      vdydt.owner = 0;
      dydtv.wrap_gsl_vector_without_ownership( &vdydt );
      vdfdy.size = size * size;
      vdfdy.stride = 1;
      vdfdy.data = 0;
      vdfdy.block = 0;
      vdfdy.owner = 0;
      dfdyv.wrap_gsl_vector_without_ownership( &vdfdy );
      mdfdy.size1 = size;
      mdfdy.size2 = size;
      mdfdy.tda = size;
      mdfdy.data = 0;
      mdfdy.block = 0;
      mdfdy.owner = 0;
      dfdym.wrap_gsl_matrix_without_ownership( &mdfdy );
      vdfdt.size = size;
      vdfdt.stride = 1;
      vdfdt.data = 0;
      vdfdt.block = 0;
      vdfdt.owner = 0;
      dfdtv.wrap_gsl_vector_without_ownership( &vdfdt );
    }
      public:
    static int function( double const t, double const* y, double* dydt, void* params ){
      Function<T>* this_t = reinterpret_cast<Function<T>*>( params );
      return this_t->function_detail( t, y, dydt );
    }
    static int jacobian( double const t, double const* y, double* dfdy, double* dfdt,
                 void* params ){
      Function<T>* this_t = reinterpret_cast<Function<T>*>( params );
      return this_t->jacobian_detail( t, y, dfdy, dfdt );
    }
      private:
    template<bool N> struct bool_to_type {};
    template<typename V>
    struct Fvector {
      template<typename W,int (V::*)(double const,gsl::vector const&,gsl::vector&)> struct S{};
      template<typename W> static char test( S<W,&W::function>* );
      template<typename W,int (V::*)(double const,gsl::vector&,gsl::vector&)> struct Su{};
      template<typename W> static char test( Su<W,&W::function>* );
      template<typename W> static int test( ... );
      static const bool vector = sizeof(test<V>( 0 )) == sizeof(char);
    };
    enum { ptr_ptr = 0, ptr_vec = 1, ptr_mat = 2, vec_ptr = 4, vec_vec = 5, vec_mat = 6,
           undefined = 8 };
    template<typename V>
    struct Jacobian {
      template<typename W,int (V::*)(double const,double const*,double*,double const*)>
      struct PtrPtr{};
      template<typename W> static char test_ptr_ptr( PtrPtr<W,&W::jacobian>* );
      template<typename W,int (V::*)(double const,double*,double*,double*)>
      struct UPtrPtr{};
      template<typename W> static char test_ptr_ptr( UPtrPtr<W,&W::jacobian>* );
      template<typename W> static int test_ptr_ptr( ... );
      template<typename W,int (V::*)(double const,double const*,gsl::vector&,double*)>
      struct PtrVec{};
      template<typename W> static char test_ptr_vec( PtrVec<W,&W::jacobian>* );
      template<typename W,int (V::*)(double const,double*,gsl::vector&,double*)>
      struct UPtrVec{};
      template<typename W> static char test_ptr_vec( UPtrVec<W,&W::jacobian>* );
      template<typename W> static int test_ptr_vec( ... );
      template<typename W,int (V::*)(double const,double const*,gsl::matrix&,double*)>
      struct PtrMat{};
      template<typename W> static char test_ptr_mat( PtrMat<W,&W::jacobian>* );
      template<typename W,int (V::*)(double const,double*,gsl::matrix&,double*)>
      struct UPtrMat{};
      template<typename W> static char test_ptr_mat( UPtrMat<W,&W::jacobian>* );
      template<typename W> static int test_ptr_mat( ... );
      template<typename W,int (V::*)(double const,gsl::vector const&,double*,gsl::vector&)>
      struct VecPtr{};
      template<typename W> static char test_vec_ptr( VecPtr<W,&W::jacobian>* );
      template<typename W,int (V::*)(double const,gsl::vector&,double*,gsl::vector&)>
      struct UVecPtr{};
      template<typename W> static char test_vec_ptr( UVecPtr<W,&W::jacobian>* );
      template<typename W> static int test_vec_ptr( ... );
      template<typename W,int (V::*)(double const,gsl::vector const&,gsl::vector&,gsl::vector&)>
      struct VecVec{};
      template<typename W> static char test_vec_vec( VecVec<W,&W::jacobian>* );
      template<typename W,int (V::*)(double const,gsl::vector&,gsl::vector&,gsl::vector&)>
      struct UVecVec{};
      template<typename W> static char test_vec_vec( UVecVec<W,&W::jacobian>* );
      template<typename W> static int test_vec_vec( ... );
      template<typename W,int (V::*)(double const,gsl::vector const&,gsl::matrix,gsl::vector&)>
      struct VecMat{};
      template<typename W> static char test_vec_mat( VecMat<W,&W::jacobian>* );
      template<typename W,int (V::*)(double const,gsl::vector&,gsl::matrix,gsl::vector&)>
      struct UVecMat{};
      template<typename W> static char test_vec_mat( UVecMat<W,&W::jacobian>* );
      template<typename W> static int test_vec_mat( ... );
      static const int type
      = sizeof(test_ptr_ptr<V>( 0 )) == sizeof(char) ? ptr_ptr :
        sizeof(test_ptr_vec<V>( 0 )) == sizeof(char) ? ptr_vec :
        sizeof(test_ptr_mat<V>( 0 )) == sizeof(char) ? ptr_mat :
        sizeof(test_vec_ptr<V>( 0 )) == sizeof(char) ? vec_ptr :
        sizeof(test_vec_vec<V>( 0 )) == sizeof(char) ? vec_vec :
        sizeof(test_vec_mat<V>( 0 )) == sizeof(char) ? vec_mat :
        undefined;
    };
    bool has_jacobian() const { return Jacobian<T>::type != undefined; }
    template<int N> struct int_to_type {};
    template<typename U>
    int function( double const t, double const* y, double* dydt, U& u, bool_to_type<true> ){
      vy.data = const_cast<double*>( y );
      vdydt.data = dydt;
      return u.function( t, yv, dydtv ); 
    }
    template<typename U>
    int function( double const t, double const* y, double* dydt, U& u, bool_to_type<false> ){
      return u.function( t, y, dydt ); 
    }
    int function_detail( double const t, double const* y, double* dydt ){
      return function( t, y, dydt, this->t, bool_to_type<Fvector<T>::vector>() );
    }
    template<typename U>
    int jacobian( double const t, double const* y, double* dfdy, double* dfdt,
              U& u, int_to_type<undefined> ){
      return 0;
    }
    template<typename U>
    int jacobian( double const t, double const* y, double* dfdy, double* dfdt,
              U& u, int_to_type<ptr_ptr> ){
      return u.jacobian( t, y, dfdy, dfdt );
    }
    template<typename U>
    int jacobian( double const t, double const* y, gsl::vector& dfdy, double* dfdt,
              U& u, int_to_type<ptr_vec> ){
      vdfdy.data = dfdy; 
      return u.jacobian( t, y, dfdyv, dfdt );
    }
    template<typename U>
    int jacobian( double const t, double const* y, gsl::vector& dfdy, double* dfdt,
              U& u, int_to_type<ptr_mat> ){
      mdfdy.data = dfdy; 
      return u.jacobian( t, y, dfdym, dfdt );
    }
    template<typename U>
    int jacobian( double const t, gsl::vector const& y, double* dfdy, gsl::vector& dfdt,
              U& u, int_to_type<vec_ptr> ){
      vdfdt.data = dfdt;
      vy.data = y;
      return u.jacobian( t, yv, dfdy, dfdtv );
    }
    template<typename U>
    int jacobian( double const t, gsl::vector const& y, gsl::vector& dfdy, gsl::vector& dfdt,
              U& u, int_to_type<vec_vec> ){
      vdfdt.data = dfdt; 
      vy.data = y;
      vdfdy.data = dfdy; 
      return u.jacobian( t, yv, dfdyv, dfdtv );
    }
    template<typename U>
    int jacobian( double const t, gsl::vector const& y, gsl::matrix& dfdy, gsl::vector& dfdt,
              U& u, int_to_type<vec_mat> ){
      vdfdt.data = dfdt; 
      vy.data = y;
      mdfdy.data = dfdy; 
      return u.jacobian( t, yv, dfdym, dfdtv );
    }
    int jacobian_detail( double const t, double const* y, double* dfdy, double* dfdt ){
      return jacobian( t, y, dfdy, dfdt, this->t, int_to_type<Jacobian<T>::type>() );
    }
    gsl_vector vy;
    gsl::vector yv;
    gsl_vector vdydt;
    gsl::vector dydtv;
    gsl_vector vdfdy;
    gsl::vector dfdyv;
    gsl_matrix mdfdy;
    gsl::matrix dfdym;
    gsl_vector vdfdt;
    gsl::vector dfdtv;
    T& t;
      };
#endif
      base_F* f;
      size_t* count;
    public:
      system(){
    f = 0;
    function = 0;
    jacobian = 0;
    dimension = 0;
    params = 0;
    count = 0; // initially nullptr will do
      }
      template<typename T>
      friend void system_constructor( system& sys, T& t );
      explicit system( gsl_odeiv2_system& v ){
    // explicitly set f = 0 to indicate nothing to be deleted.
    f = 0;
    // copy
    function = v.function;
    jacobian = v.jacobian;
    dimension = v.dimension;
    params = v.params;
    // 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 function
      }
      template<typename T>
      system( T& t ){ system_constructor<T>( *this, t ); }
       
      // Copy and move constructors
      // copy constructor
      system( system const& v ) : f( v.f ), count( v.count ){
    function = v.function;
    jacobian = v.jacobian;
    dimension = v.dimension;
    params = v.params;
    if( count != 0 ) ++*count; // system is now shared.
      }
      // assignment operator
      system& operator=( system const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      delete f;
      delete count;
    }
    // Then copy
    f = v.f;
    function = v.function;
    jacobian = v.jacobian;
    dimension = v.dimension;
    params = v.params;
    // count = v.count;
    if( count != 0 ) ++*count; // system is now shared.
    return *this;
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      system( system&& v ) : f( v.f ), count( nullptr ){
    std::swap( function, v.function );
    std::swap( jacobian, v.jacobian );
    std::swap( dimension, v.dimension );
    std::swap( params, v.params );
    std::swap( count, v.count );
    v.function = nullptr;
      }
      system& operator=( system&& v ){
    std::swap( f, v.f );
    std::swap( function, v.function );
    std::swap( jacobian, v.jacobian );
    std::swap( dimension, v.dimension );
    std::swap( params, v.params );
    std::swap( count, v.count );
    return *this;
      }
#endif
      ~system(){ 
    // first, possibly delete anything pointed to by f
    if( count == 0 or --*count == 0 ){
      delete f;
      delete count;
    }
      }
    };
#ifndef DOXYGEN_SKIP
    template<typename T>
    inline void system_constructor( system& sys, T& t ){
      sys.function = 0;
      sys.count = 0;
      // try constructing from T
      sys.f = new system::Function<T>( t );
      sys.function = &system::Function<T>::function;
      if( sys.f->has_jacobian() )
    sys.jacobian = 0;
      else
    sys.jacobian = &system::Function<T>::jacobian;
      sys.dimension = t.size();
      sys.params = sys.f;
      // just plausibly we could fail to allocate count: no further action needed.
      sys.count = new size_t;
      *sys.count = 1; // initially there is just one reference to function
    }
    template<>
    inline void system_constructor( system& sys, system& v ){
      sys.f = v.f; sys.function = v.function; sys.jacobian = v.jacobian;
      sys.dimension = v.dimension; sys.params = v.params;
      if( sys.count != 0 ) ++*sys.count; // system is now shared.
      /* sys is now shared. */ }
    template<>
    inline void system_constructor( system& sys, gsl_odeiv2_system& v ){
      sys.function = v.function; sys.jacobian = v.jacobian;
      sys.dimension = v.dimension; sys.params = v.params; }
#endif
  
    template<typename T>
    inline system make_system( T& t ){
      system fn( t );
      return fn;
    }
      
  }
}
#endif