function_fdf.hpp

Go to the documentation of this file.

/*
 * $Id: function_fdf.hpp 51 2012-01-07 13:55:48Z jdl3 $
 * Copyright (C) 2010, 2011, 2012, 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_FUNCTION_FDF_HPP
#define CCGSL_FUNCTION_FDF_HPP
 
// For std::pair
#include<utility>
// For std::swap
#if __cplusplus <= 199711L
#include<algorithm>
#endif
#include"function_scl.hpp"
 
namespace gsl {
  class function_fdf;
  template<typename T>
  void function_fdf_constructor( function_fdf& fdf, T& t );
  class function_fdf : public gsl_function_fdf {
  public:
    struct Concept {
      virtual double f( double x ) = 0;
      virtual double df( double x ) = 0;
      virtual std::pair<double,double> fdf( double x ) = 0;
    };
#ifndef DOXYGEN_SKIP
  private:
    struct base_F {
      typedef double (*function_t)(double,void*);
      typedef void (*function_fdf_t)(double,void*,double*,double*);
    };
    template<typename T>
    struct FDF : public base_F {
    public:
      FDF( T& t ) : t( t ){}
    public:
      static double fn( double x, void* params ){
        return reinterpret_cast<FDF<T>*>( params )->t.f( x ); }
      static double dfn( double x, void* params ){
        return reinterpret_cast<FDF<T>*>( params )->t.df( x ); }
      static void fdfn( double x, void* params, double* d, double* df ){
    std::pair<double,double> pair = reinterpret_cast<FDF<T>*>( params )->t.fdf( x );
    *d = pair.first; *df = pair.second; }
    private:
      T& t;
    };
    struct function_scl_adapter : public base_F, public Concept {
    public:
      function_scl_adapter( function_scl& function, function_scl&  dfunction )
    : scl_function( function ), scl_dfunction( dfunction ){}
      double f( double x ){ return fn_eval( scl_function, x ); }
      double df( double x ){ return fn_eval( scl_dfunction, x ); }
      std::pair<double,double> fdf( double x ){
    return std::make_pair( fn_eval( scl_function, x ), fn_eval( scl_dfunction, x ) ); }
    public:
      static double fn( double x, void* params ){
        return reinterpret_cast<function_scl_adapter*>( params )->f( x ); }
      static double dfn( double x, void* params ){
        return reinterpret_cast<function_scl_adapter*>( params )->df( x ); }
      static void fdfn( double x, void* params, double* d, double* df ){
    std::pair<double,double> pair = reinterpret_cast<function_scl_adapter*>( params )->fdf( x );
    *d = pair.first; *df = pair.second; }
    private:
      function_scl& scl_function;
      function_scl& scl_dfunction;
    };
#endif
    base_F* function;
    size_t* count;
  public:
    function_fdf(){
      function = 0;
      f = 0;
      df = 0;
      fdf = 0;
      params = 0;
      count = 0; // initially nullptr will do
    }
    template<typename T>
    friend void function_fdf_constructor( function_fdf& fdf, T& t );
    explicit function_fdf( gsl_function_fdf& v ){
      // explicitly set function = 0 to indicate nothing to be deleted.
      function = 0;
      // copy
      f = v.f;
      df = v.df;
      fdf = v.fdf;
      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_fdf( T& t ){
      function_fdf_constructor<T>( *this, t );
    }
    function_fdf( function_scl& function, function_scl&  dfunction ){
      this->function = 0;
      count = 0;
      // try constructing from function_scl_adapter
      this->function = new function_scl_adapter( function, dfunction );
      f = &function_scl_adapter::fn;
      df = &function_scl_adapter::dfn;
      fdf = &function_scl_adapter::fdfn;
      params = this->function;
      // 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_fdf( function_fdf const& v ) : function( v.function ), count( v.count ){
      f = v.f;
      df = v.df;
      fdf = v.fdf;
      params = v.params;
      if( count != 0 ) ++*count; // function_fdf is now shared.
    }
    // assignment operator
    function_fdf& operator=( function_fdf const& v ){
      // first, possibly delete anything pointed to by this
      if( count == 0 or --*count == 0 ){
    delete function;
    delete count;
      }
      // Then copy
      function = v.function;
      f = v.f;
      df = v.df;
      fdf = v.fdf;
      params = v.params;
      // count = v.count;
      if( count != 0 ) ++*count; // function_fdf is now shared.
      return *this;
    }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
    function_fdf( function_fdf&& v ) : function( v.function ), count( nullptr ){
      std::swap( f, v.f );
      std::swap( df, v.df );
      std::swap( fdf, v.fdf );
      std::swap( params, v.params );
      std::swap( count, v.count );
      v.function = nullptr;
    }
    function_fdf& operator=( function_fdf&& v ){
      std::swap( function, v.function );
      std::swap( f, v.f );
      std::swap( df, v.df );
      std::swap( fdf, v.fdf );
      std::swap( params, v.params );
      std::swap( count, v.count );
      return *this;
    }
#endif
    ~function_fdf(){ 
      // first, possibly delete anything pointed to by f
      if( count == 0 or --*count == 0 ){
    delete function;
    delete count;
      }
    }
  };
#ifndef DOXYGEN_SKIP
  template<typename T>
  inline void function_fdf_constructor( function_fdf& fdf, T& t ){
    fdf.function = 0;
    fdf.count = 0;
    // try constructing from T
    fdf.function = new function_fdf::FDF<T>( t );
    fdf.f = &function_fdf::FDF<T>::fn;
    fdf.df = &function_fdf::FDF<T>::dfn;
    fdf.fdf = &function_fdf::FDF<T>::fdfn;
    fdf.params = fdf.function;
    // just plausibly we could fail to allocate count: no further action needed.
    fdf.count = new size_t;
    *fdf.count = 1; // initially there is just one reference to function
  }
  template<>
  inline void function_fdf_constructor( function_fdf& fdf, function_fdf& v ){
    fdf.function = v.function; fdf.count = v.count; fdf.f = v.f; fdf.df = v.df;
    fdf.fdf = v.fdf; fdf.params = v.params; if( fdf.count != 0 ) ++*fdf.count;
    /* function_fdf is now shared. */ }
  template<>
  inline void function_fdf_constructor( function_fdf& fdf, gsl_function_fdf& v ){
    fdf.f = v.f; fdf.df = v.df; fdf.fdf = v.fdf; fdf.params = v.params; }
#endif
  
  template<typename T>
  inline function_fdf make_function_fdf( T& t ){
    function_fdf fn( t );
    return fn;
  }
 
  inline double fn_fdf_eval_f( gsl_function_fdf const* const F, double const x ){
    return (F->f)( x, F->params );
  }
  
  inline double fn_fdf_eval_f( gsl_function_fdf const& F, double const x ){
    return (F.f)( x, F.params );
  }
 
  inline double fn_fdf_eval_df( gsl_function_fdf const* const F, double const x ){
    return (F->df)( x, F->params );
  }
  
  inline double fn_fdf_eval_df( gsl_function_fdf const& F, double const x ){
    return (F.df)( x, F.params );
  }
 
  inline void fn_fdf_eval_df( gsl_function_fdf const* const F, double const x,
                  double* f, double* df ){
    return (F->fdf)( x, F->params, f, df );
  }
  
  inline void fn_fdf_eval_df( gsl_function_fdf const& F, double const x,
                  double& f, double& df ){
    return (F.fdf)( x, F.params, &f, &df );
  }
}
 
#endif