multilarge_nlinear_function_fdf.hpp

Go to the documentation of this file.

/*
 * $Id
 * Copyright (C) 2010, 2011, 2012, 2020 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_MULTILARGE_NLINEAR_FUNCTION_FDF_HPP
#define CCGSL_MULTILARGE_NLINEAR_FUNCTION_FDF_HPP
 
// For std::swap
#if __cplusplus <= 199711L
#include<algorithm>
#else
#include<utility>
#endif
#include<gsl/gsl_multilarge_nlinear.h>
#include"vector.hpp"
#include"matrix.hpp"
 
namespace gsl {
  namespace multilarge {
    namespace nlinear {
      class function_fdf;
      template<typename T> void function_constructor( function_fdf& f, T& t );
      class function_fdf : public gsl_multilarge_nlinear_fdf {
      public:
    struct concept_fdf {
      virtual int f( gsl::vector const& x, gsl::vector& f ) = 0;
      virtual int df( CBLAS_TRANSPOSE_t TransJ, gsl::vector const& x, gsl::vector const& u,
              gsl::vector& v,gsl::matrix& JTJ ) = 0;
      virtual size_t xSize() const = 0;
      virtual size_t fSize() const = 0;
    };
    struct concept_fvv {
      virtual int fvv( gsl::vector const& x, gsl::vector const& v, gsl::vector& fvv ) = 0;
    };
#ifndef DOXYGEN_SKIP
    struct address_of {
      typedef int (*f_t)(gsl_vector const*,void*,gsl_vector*);
      typedef int (*df_t)(CBLAS_TRANSPOSE_t,gsl_vector const*,
                  gsl_vector const*,void*,gsl_vector*,gsl_matrix*);
      typedef int (*fvv_t)(gsl_vector const*,gsl_vector const*,void*,gsl_vector*);
      // Maybe not needed: f must not be nullptr
      template<typename U, int (U::*)(gsl::vector const&, gsl::vector&)>
      struct fSfinae {};
      template<typename U, int (U::*)(gsl::vector const&, gsl::vector&) const>
      struct fSfinaeC {};
      template<typename U> static f_t f(fSfinae<U,&U::f>*){return &F<U>::fn;}
      template<typename U> static f_t f(fSfinaeC<U,&U::f>*){return &F<U>::fn;}
      template<typename U> static f_t f(...){return 0;}
      // address_of<X>::f = nonzero if f implemented
      template<typename U, int (U::*)(CBLAS_TRANSPOSE_t,gsl::vector const&,
                      gsl::vector const&,gsl_vector*,gsl::matrix&)>
      struct dfSfinae {};
      template<typename U, int (U::*)(CBLAS_TRANSPOSE_t,gsl::vector const&,
                      gsl::vector const&,gsl_vector*,gsl::matrix&) const>
      struct dfSfinaeC {};
      template<typename U> static df_t df(dfSfinae<U,&U::df>*){return &F<U>::dfn;}
      template<typename U> static df_t df(dfSfinaeC<U,&U::df>*){return &F<U>::dfn;}
      template<typename U> static df_t df(...){return 0;}
      // address_of<X>::df = nonzero if df implemented
      template<typename U, int (U::*)(gsl::vector const&, gsl::vector const&, gsl::vector&)>
      struct fvvSfinae {};
      template<typename U, int (U::*)(gsl::vector const&, gsl::vector const&, gsl::vector&)
           const>
      struct fvvSfinaeC {};
      template<typename U> static fvv_t fvv(fvvSfinae<U,&U::fvv>*){return &F<U>::fvvn;}
      template<typename U> static fvv_t fvv(fvvSfinaeC<U,&U::fvv>*){return &F<U>::fvvn;}
      template<typename U> static fvv_t fvv(...){return 0;}
      // address_of<X>::fvv = nonzero if fvv implemented
    };
      private:
    struct base_F {
      typedef int (*function_t)(gsl_vector const*,void*,gsl_vector*);
      typedef int (*dfunction_t)(CBLAS_TRANSPOSE_t,gsl_vector const*,
                     gsl_vector const*,void*,gsl_vector*,gsl_matrix*);
      typedef int (*fvvfunction_t)(gsl_vector const*,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 );
        uv.wrap_gsl_vector_without_ownership( 0 );
        vv.wrap_gsl_vector_without_ownership( 0 );
        fv.wrap_gsl_vector_without_ownership( 0 );
        JTJv.wrap_gsl_matrix_without_ownership( 0 );
      }
    public:
      static int fn( gsl_vector const* x, void* params, gsl_vector* f ){
        // 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 = f;
        return ft->t.f( ft->xv, ft->fv ); }
      static int dfn( CBLAS_TRANSPOSE_t TransJ, gsl_vector const* x, gsl_vector const* u,
              void* params, gsl_vector* v, gsl_matrix* JTJ ){
        // 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->uv.ccgsl_pointer = const_cast<gsl_vector*>( u );
        ft->vv.ccgsl_pointer = v;
        ft->JTJv.ccgsl_pointer = JTJ;
        return ft->t.df( TransJ, ft->xv, ft->uv, ft->vv, ft->JTJv ); }
      static int fvvn( gsl_vector const* x, gsl_vector const* v,
               void* params, gsl_vector* f ){
        // 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->vv.ccgsl_pointer = const_cast<gsl_vector*>( v );
        ft->fv.ccgsl_pointer = f;
        return ft->t.fvv( ft->xv, ft->vv, ft->fv ); }
    private:
      T& t;
      gsl::vector xv;
      gsl::vector uv;
      gsl::vector vv;
      gsl::vector fv;
      gsl::matrix JTJv;
    };
    class Fn : public base_F {
    public:
      Fn( int (*const f)(gsl::vector const&,gsl::vector&),
          int (*const df)(CBLAS_TRANSPOSE_t,gsl::vector const&,gsl::vector const&,
                  gsl::vector&,gsl::matrix&),
          int (*const fvv)(gsl::vector const&,gsl::vector&,gsl::vector&),
          size_t const xSize, size_t fSize )
        : f( f ), df( df ), fvv( fvv ), xSize( xSize ), fSize( fSize ){}
      function_t function(){ return &fn; }
      dfunction_t dfunction(){ return &dfn; }
      fvvfunction_t fvvfunction(){ return &fvvn; }
    private:
      int (*const f)(gsl::vector const&,gsl::vector&);
      int (*const df)(CBLAS_TRANSPOSE_t,gsl::vector const&,gsl::vector const&,
              gsl::vector&,gsl::matrix&);
      int (*const fvv)(gsl::vector const&,gsl::vector&,gsl::vector&);
    public:
      size_t const xSize;
      size_t const fSize;
      static int fn( gsl_vector const* x, void* params, gsl_vector* f ){
        // share gsl_vector* and gsl::vector contents.
        Fn* ft = reinterpret_cast<Fn*>( params );
        ft->fv.ccgsl_pointer = f;
        return ft->f( ft->xv, ft->fv ); }
      static int dfn( CBLAS_TRANSPOSE_t TransJ, gsl_vector const* x, gsl_vector const* u,
              void* params, gsl_vector* v, gsl_matrix* JTJ ){
        // share gsl_vector* and gsl::vector contents.
        Fn* ft = reinterpret_cast<Fn*>( params );
        ft->xv.ccgsl_pointer = const_cast<gsl_vector*>( x );
        ft->uv.ccgsl_pointer = const_cast<gsl_vector*>( u );
        ft->vv.ccgsl_pointer = v;
        ft->JTJv.ccgsl_pointer = JTJ;
        return ft->df( TransJ, ft->xv, ft->uv, ft->vv, ft->JTJv ); }
      static int fvvn( gsl_vector const* x, gsl_vector const* v,
               void* params, gsl_vector* f ){
        // share gsl_vector* and gsl::vector contents.
        Fn* ft = reinterpret_cast<Fn*>( params );
        ft->xv.ccgsl_pointer = const_cast<gsl_vector*>( x );
        ft->vv.ccgsl_pointer = const_cast<gsl_vector*>( v );
        ft->fv.ccgsl_pointer = f;
        return ft->fvv( ft->xv, ft->vv, ft->fv ); }
    private:
      gsl::vector xv;
      gsl::vector vv;
      gsl::vector uv;
      gsl::vector fv;
      gsl::matrix JTJv;
    };
#endif
    base_F* func;
    size_t* count;
      public:
    function_fdf(){
      func = 0;
      f = 0;
      df = 0;
      fvv = 0;
      p = 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_multilarge_nlinear_fdf& v ){
      // explicitly set function = 0 to indicate nothing to be deleted.
      func = 0;
      // copy
      f = v.f;
      df = v.df;
      fvv = v.fvv;
      p = v.p;
      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)(CBLAS_TRANSPOSE_t,gsl::vector const&, gsl::vector const&,
                       gsl::vector&, gsl::matrix&),
              int (* const fvv)(gsl::vector const&, gsl::vector&, gsl::vector&),
              size_t const xSize, size_t const fSize ){
      func = 0;
      count = 0;
      func = new function_fdf::Fn( f, df, fvv, xSize, fSize );
      this->f = (0 == f) ? 0 : Fn::fn;
      this->df = (0 == df) ? 0 : Fn::dfn;
      this->fvv = (0 == fvv) ? 0 : Fn::fvvn;
      this->n = fSize;
      this->p = xSize;
      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;
      fvv = v.fvv;
      n = v.n;
      p = v.p;
      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;
      fvv = v.fvv;
      n = v.n;
      p = v.p;
      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( fvv, v.fvv );
      std::swap( n, v.n );
      std::swap( p, v.p );
      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( fvv, v.fvv );
      std::swap( n, v.n );
      std::swap( p, v.p );
      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::address_of::f<T>(0);
    f.df = function_fdf::address_of::df<T>(0);
    f.fvv = function_fdf::address_of::fvv<T>(0);
    // Checks coplete
    f.n = t.fSize();
    f.p = t.xSize();
    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.fvv = v.fvv;
    f.n = v.n; f.p = v.p; f.params = v.params; if( f.count != 0 ) ++*f.count;
    /* function is now shared. */ }
      template<>
      inline void function_constructor( function_fdf& f, gsl_multilarge_nlinear_fdf& v ){
    f.f = v.f; f.df = v.df; f.fvv = v.fvv; f.n = v.n; f.p = v.p; f.params = v.params; }
#endif
      template<typename T>
      inline function_fdf make_function_fdf( T& t ){
    function_fdf fn( t );
    return fn;
      }
    }
  }
}
 
#endif