multiroot_function_fdf.hpp

Go to the documentation of this file.

/*
 * $Id: multiroot_function_fdf.hpp 171 2012-06-19 19:15:15Z 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_MULTIROOT_FUNCTION_FDF_HPP
#define CCGSL_MULTIROOT_FUNCTION_FDF_HPP
 
// For std::swap
#if __cplusplus <= 199711L
#include<algorithm>
#else
#include<utility>
#endif
 
namespace gsl {
  namespace multiroot {
    class function_fdf : public gsl_multiroot_function_fdf {
    public:
      struct Concept {
    virtual int f( gsl::vector const& x, gsl::vector& f ) = 0;
    virtual int df( gsl::vector const& x, gsl::matrix& J ) = 0;
    virtual int fdf( gsl::vector const& x, gsl::vector& f, gsl::matrix& J ) = 0;
    virtual size_t size() const = 0;
      };
#ifndef DOXYGEN_SKIP
    private:
      struct base_F {
    typedef int (*function_t)(gsl_vector const*,void*,gsl_vector*);
    typedef int (*dfunction_t)(gsl_vector const*,void*,gsl_matrix*);
    typedef int (*fdfunction_t)(gsl_vector const*,void*,gsl_vector*,gsl_matrix*);
      };
      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 );
      dfv.wrap_gsl_matrix_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 ); }
    static int dfn( gsl_vector const* x, void* params, gsl_matrix* dfx ){
      // 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->dfv.ccgsl_pointer = dfx;
      return ft->t.df( ft->xv, ft->dfv ); }
    static int fdfn( gsl_vector const* x, void* params, gsl_vector* fx, gsl_matrix* dfx ){
      // 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;
      ft->dfv.ccgsl_pointer = dfx;
      return ft->t.fdf( ft->xv, ft->fv, ft->dfv ); }
      private:
    T& t;
    gsl::vector xv;
    gsl::vector fv;
    gsl::matrix dfv;
      };
      class Fn : public base_F {
      public:
    Fn( int (*const f)(gsl::vector const&,gsl::vector&),
        int (*const df)(gsl::vector const&,gsl::matrix&),
        int (*const fdf)(gsl::vector const&,gsl::vector&,gsl::matrix&),
        size_t n ) : f( f ), df( df ), fdf( fdf ), n( n ){}
    function_t function(){ return &fn; }
    dfunction_t dfunction(){ return &dfn; }
    fdfunction_t fdfunction(){ return &fdfn; }
      private:
    int (*const f)(gsl::vector const&,gsl::vector&);
    int (*const df)(gsl::vector const&,gsl::matrix&);
    int (*const fdf)(gsl::vector const&,gsl::vector&,gsl::matrix&);
      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->fv.ccgsl_pointer = fx;
      return ft->f( ft->xv, ft->fv ); }
    static int dfn( gsl_vector const* x, void* params, gsl_matrix* dfx ){
      // share gsl_vector* and gsl::vector contents.
      Fn* ft = reinterpret_cast<Fn*>( params );
      ft->xv.ccgsl_pointer = const_cast<gsl_vector*>( x );
      ft->dfv.ccgsl_pointer = dfx;
      return ft->df( ft->xv, ft->dfv ); }
    static int fdfn( gsl_vector const* x, void* params, gsl_vector* fx, gsl_matrix* dfx ){
      // 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;
      ft->dfv.ccgsl_pointer = dfx;
      return ft->fdf( ft->xv, ft->fv, ft->dfv ); }
      private:
    gsl::vector xv;
    gsl::vector fv;
    gsl::matrix dfv;
      };
      template<typename T>
      struct gsl_F : public base_F {
      public:
    gsl_F( T& t ) : t( t ){}
      public:
    static int fn( gsl_vector const* x, void* params, gsl_vector* f ){
      return reinterpret_cast<F<T>*>( params )->t.f( x, f ); }
    static int dfn( gsl_vector const* x, void* params, gsl_matrix* df ){
      return reinterpret_cast<F<T>*>( params )->t.df( x, df ); }
    static int fdfn( gsl_vector const* x, void* params, gsl_vector* f, gsl_matrix* df ){
      return reinterpret_cast<F<T>*>( params )->t.fdf( x, f, df ); }
      private:
    T& t;
      };
      class gsl_Fn : public base_F {
      public:
    gsl_Fn( int (*const f)(gsl_vector const*,gsl_vector*),
        int (*const df)(gsl_vector const*,gsl_matrix*),
        int (*const fdf)(gsl_vector const*,gsl_vector*,gsl_matrix*),
        size_t n ) : f( f ), df( df ), fdf( fdf ), n( n ){}
    function_t function(){ return &fn; }
    dfunction_t dfunction(){ return &dfn; }
    fdfunction_t fdfunction(){ return &fdfn; }
      private:
    int (*const f)(gsl_vector const*,gsl_vector*);
    int (*const df)(gsl_vector const*,gsl_matrix*);
    int (*const fdf)(gsl_vector const*,gsl_vector*,gsl_matrix*);
      public:
    size_t const n;
    static int fn( gsl_vector const* x, void* params, gsl_vector* f ){
      return reinterpret_cast<gsl_Fn*>( params )->f( x, f ); }
    static int dfn( gsl_vector const* x, void* params, gsl_matrix* df ){
      return reinterpret_cast<gsl_Fn*>( params )->df( x, df ); }
    static int fdfn( gsl_vector const* x, void* params, gsl_vector* f, gsl_matrix* df ){
      return reinterpret_cast<gsl_Fn*>( params )->fdf( x, f, df ); }
      };
#endif
      base_F* func;
      size_t* count;
    public:
      function_fdf(){
    func = 0;
    f = 0;
    df = 0;
    fdf = 0;
    n = 0;
    params = 0;
    count = 0; // initially nullptr will do
      }
      template<typename T>
      friend void function_constructor( function_fdf&, T& );
      explicit function_fdf( gsl_multiroot_function_fdf& v ){
    // explicitly set function = 0 to indicate nothing to be deleted.
    func = 0;
    // copy
    f = v.f;
    df = v.df;
    fdf = v.fdf;
    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_fdf( T& t ){ function_constructor<T>( *this, t ); }
      function_fdf( int (* const f)(gsl::vector const&, gsl::vector&),
            int (* const df)(gsl::vector const&, gsl::matrix&),
            int (* const fdf)(gsl::vector const&, gsl::vector&, gsl::matrix&),
            size_t n ){
    func = 0;
    count = 0;
    func = new function_fdf::Fn( f, df, fdf, n );
    this->f = Fn::fn;
    this->df = Fn::dfn;
    this->fdf = Fn::fdfn;
    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_fdf( function_fdf const& v ) : func( v.func ), count( v.count ){
    f = v.f;
    df = v.df;
    fdf = v.fdf;
    n = v.n;
    params = v.params;
    if( count != 0 ) ++*count; // function 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 func;
      delete count;
    }
    // Then copy
    func = v.func;
    f = v.f;
    df = v.df;
    fdf = v.fdf;
    n = v.n;
    params = v.params;
    // count = v.count;
    if( count != 0 ) ++*count; // function is now shared.
    return *this;
      }
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      function_fdf( function_fdf&& v ) : func( v.func ), count( nullptr ){
    std::swap( f, v.f );
    std::swap( df, v.df );
    std::swap( fdf, v.fdf );
    std::swap( n, v.n );
    std::swap( params, v.params );
    std::swap( count, v.count );
    v.func = nullptr;
      }
      function_fdf& operator=( function_fdf&& v ){
    std::swap( func, v.func );
    std::swap( f, v.f );
    std::swap( df, v.df );
    std::swap( fdf, v.fdf );
    std::swap( n, v.n );
    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 func;
      delete count;
    }
      }
    };
#ifndef DOXYGEN_SKIP
    template<typename T>
    inline void function_constructor( function_fdf& f, T& t ){
      f.func = 0;
      f.count = 0;
      // try constructing from T
      f.func = new function_fdf::F<T>( t );
      f.f = &function_fdf::F<T>::fn;
      f.df = &function_fdf::F<T>::dfn;
      f.fdf = &function_fdf::F<T>::fdfn;
      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_fdf& f, function_fdf& v ){
      f.func = v.func; f.count = v.count; f.f = v.f; f.df = v.df; f.fdf = v.fdf; f.n = v.n;
      f.params = v.params; if( f.count != 0 ) ++*f.count;
      /* function is now shared. */ }
    template<>
    inline void function_constructor( function_fdf& f, gsl_multiroot_function_fdf& v ){
      f.f = v.f; f.df = v.df; f.fdf = v.fdf; f.n = v.n; f.params = v.params; }
#endif
 
    template<typename T>
    inline function_fdf make_function_fdf( T& t ){
      function_fdf fn( t );
#ifdef __GXX_EXPERIMENTAL_CXX0X__
      return std::move( fn );
#else
      return fn;
#endif
    }
  }
}
 
#endif