multiroot_function.hpp

Go to the documentation of this file.

/*
 * $Id: multiroot_function.hpp 235 2012-08-08 19:46:06Z jdl3 $
 * Copyright (C) 2010, 2011, 2012, 2016, 2020, 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_MULTIROOT_FUNCTION_HPP
#define CCGSL_MULTIROOT_FUNCTION_HPP
 
#include<gsl/gsl_multiroots.h>
// For std::swap
#if __cplusplus <= 199711L
#include<algorithm>
#else
#include<utility>
#endif
 
namespace gsl {
  namespace multiroot {
    class function;
    template<typename T>
    void function_constructor( function&, T& );
    class function : public gsl_multiroot_function {
    public:
      struct Concept {
    virtual int f( gsl::vector const& x, gsl::vector& f ) = 0;
    virtual size_t size() const = 0;
      };
#ifndef DOXYGEN_SKIP
    private:
      struct base_F {
    typedef int (*function_t)(gsl_vector const*,void*,gsl_vector*);
      };
      template<typename T>
      struct F : public base_F {
      public:
    F( T& t ) : t( t ){
      xv.wrap_gsl_vector_without_ownership( 0 );
      fv.wrap_gsl_vector_without_ownership( 0 );
    }
      public:
    static int fn( gsl_vector const* x, void* params, gsl_vector* fx ){
      // share gsl_vector* and gsl::vector contents.
      F<T>* ft = reinterpret_cast<F<T>*>( params );
      ft->xv.ccgsl_pointer = const_cast<gsl_vector*>( x );
      ft->fv.ccgsl_pointer = fx;
      return ft->t.f( ft->xv, ft->fv ); }
      private:
    T& t;
    gsl::vector xv;
    gsl::vector fv;
      };
      class Fn : public base_F {
      public:
    Fn( int (*const f)(gsl::vector const&,gsl::vector&), size_t n ) : f( f ), n( n ){
      xv.wrap_gsl_vector_without_ownership( 0 );
      fv.wrap_gsl_vector_without_ownership( 0 );
    }
    function_t function(){ return &fn; }
      private:
    int (*const f)(gsl::vector const&,gsl::vector&);
      public:
    size_t const n;
    static int fn( gsl_vector const* x, void* params, gsl_vector* fx ){
      // share gsl_vector* and gsl::vector contents.
      Fn* ft = reinterpret_cast<Fn*>( params );
      ft->xv.ccgsl_pointer = const_cast<gsl_vector*>( x );
      ft->fv.ccgsl_pointer = fx;
      return ft->f( ft->xv, ft->fv ); }
      private:
    gsl::vector xv;
    gsl::vector fv;
      };
      // /**
      //  * Used internally to point to function objects.
      //  */
      // template<typename T>
      // struct gsl_F : public base_F {
      // public:
      //    /**
      //     * The constructor requires an object that implements Concept.
      //     * @param A function object that implements Concept.
      //     */
      //    gsl_F( T& t ) : t( t ){}
      // public:
      //    /**
      //     * A static function that gsl::function_f can use.
      //     * @param x The arguments of the functions.
      //     * @param params Always a pointer to an object of type gsl_F.
      //     * @param f The values of the functions.
      //     * @return The function result.
      //     */
      //    static int fn( gsl_vector const* x,  void* params, gsl::vector* f ){
      //      return reinterpret_cast<F<T>*>( params )->t.f( x, f ); }
      // private:
      //    /**
      //     * Reference to function object from which this was constructed.
      //     */
      //    T& t;
      // };
      // /**
      //  * Subclass that allows us to use arbitrary functions.
      //  */
      // class gsl_Fn : public base_F {
      // public:
      //    /**
      //     * The constructor requires a function.
      //     * @param f A function.
      //     * @param n The number of functions and arguments.
      //     */
      //    gsl_Fn( int (*const f)(gsl_vector const*,gsl_vector*), size_t n ) : f( f ), n( n ){}
      //    /**
      //     * @return A function for gsl_multiroot_function.
      //     */
      //    function_t function(){ return &fn; }
      // private:
      //    /**
      //     * The original function.
      //     */
      //    int (*const f)(gsl_vector const*,gsl_vector*);
      // public:
      //    /**
      //     * The number of functions and arguments.
      //     */
      //    size_t const n;
      //    /**
      //     * A static function that gsl::function_scl can use.
      //     * @param x The function argument.
      //     * @param params Always a pointer to an object of type gsl_Fn
      //     * @return The function result.
      //     */
      //    static int fn( gsl_vector const* x, void* params, gsl_vector* f ){
      //      return reinterpret_cast<gsl_Fn*>( params )->f( x, f ); }
      // };
#endif
      base_F* func;
      size_t* count;
    public:
      function(){
    func = 0;
    f = 0;
    n = 0;
    params = 0;
    count = 0; // initially nullptr will do
      }
      template<typename T>
      friend void function_constructor( function&, T& );
      explicit function( gsl_multiroot_function& v ){
    // explicitly set function = 0 to indicate nothing to be deleted.
    func = 0;
    // copy
    f = v.f;
    n = v.n;
    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>
      function( T& t ){ function_constructor<T>( *this, t ); }
      function( int (* const f)(gsl::vector const&, gsl::vector&), size_t n ){
    func = 0;
    count = 0;
    func = new function::Fn( f, n );
    this->f = Fn::fn;
    this->n = n;
    params = func;
    // 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
      }
       
      // Copy and move constructors
      // copy constructor
      function( function const& v ) : func( v.func ), count( v.count ){
    f = v.f;
    n = v.n;
    params = v.params;
    if( count != 0 ) ++*count; // function is now shared.
      }
      // assignment operator
      function& operator=( function const& v ){
    // first, possibly delete anything pointed to by this
    if( count == 0 or --*count == 0 ){
      delete func;
      delete count;
    }
    // Then copy
    func = v.func;
    f = v.f;
    n = v.n;
    params = v.params;
    // count = v.count;
    if( count != 0 ) ++*count; // function is now shared.
    return *this;
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      function( function&& v ) : func( v.func ), count( nullptr ){
    std::swap( f, v.f );
    std::swap( n, v.n );
    std::swap( params, v.params );
    std::swap( count, v.count );
    v.func = nullptr;
      }
      function& operator=( function&& v ){
    std::swap( func, v.func );
    std::swap( f, v.f );
    std::swap( n, v.n );
    std::swap( params, v.params );
    std::swap( count, v.count );
    return *this;
      }
#endif
      ~function(){ 
    // first, possibly delete anything pointed to by f
    if( count == 0 or --*count == 0 ){
      delete func;
      delete count;
    }
      }
    };
#ifndef DOXYGEN_SKIP
    template<typename T>
    inline void function_constructor( function& f, T& t ){
      f.func = 0;
      f.count = 0;
      // try constructing from T
      f.func = new function::F<T>( t );
      f.f = &function::F<T>::fn;
      f.n = t.size();
      f.params = f.func;
      // just plausibly we could fail to allocate count: no further action needed.
      f.count = new size_t;
      *f.count = 1; // initially there is just one reference to function
    }
    template<>
    inline void function_constructor( function& f, function& v ){
      f.func = v.func; f.count = v.count; f.f = v.f; f.n = v.n;
      f.params = v.params; if( f.count != 0 ) ++*f.count;
      /* function is now shared. */ }
    template<>
    inline void function_constructor( function& f, gsl_multiroot_function& v ){
      f.f = v.f; f.n = v.n; f.params = v.params; }
#endif
 
    template<typename T>
    inline function make_function( T& t ){
      function fn( t );
      return fn;
    }
  }
}
 
#endif