[Help-gsl] Re: gsl ode solver - const a problem?

[Heiko Bauke]

Dear Brian,

On Mon, 21 Sep 2009 09:19:49 +0100
Brian Gough <address@hidden> wrote:

> At Fri, 18 Sep 2009 10:42:13 +0200,
> Heiko Bauke wrote:
> > There is no (obvious) reason why func should change the parameters. 
> 
> It is mainly for generality.  A simple example used in the GSL test
> programs is counting function evaluations.

I understand your desire to provide a library that is as general as
possible and does not limit its user by its interface. However, I think
counting function evaluations is a rather artificial example. One might
call this example a miss-use of the param-parameter. Quoting the GSL
manual: 

void * params
    This is a pointer to the arbitrary parameters of the system. 
 
> > The current design forces me regularly to apply a const_cast (in C+
> > +) when solving ODEs with GSL routines.
> 
> I'm not familiar with this situation, could you give an example -  is
> it specific to C++ or is there an equivalent problem in C?  Thank you.

No it's not really a C++ specific issue. However, const correctness is
a larger concern in C++, see
http://www.parashift.com/c++-faq-lite/const-correctness.html .

Imagine the gsl-routines are used in a more complex setting than in the
examples in the gsl manual. Furthermore, imagine that the parameters of
the ODE system are passes through a hierarchy of function calls via
constant pointers (or constant references in the case of C++), then
one has to cast the constness explicitly away. 

Here small program that might give you an idea what I am talking about.

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
#include <iostream>
#include <vector>

class vector {
  std::vector<double> v;
public:
  vector(int N) : v(N, 0) {
  }
  double & operator()(int i) {
    return v[i];
  }
  const double & operator()(int i) const {
    return v[i];
  }
  int size() const {
    return v.size();
  }
};

class matrix {
  int N;
  std::vector<double> A;
public:
  matrix(int N) : N(N), A(N*N, 0) {
  }
  double & operator()(int i, int j) {
    return A[i+j*N];
  }
  const double & operator()(int i, int j) const {
    return A[i+j*N];
  }
  int size() const {
    return N;
  }
};

class stepper {
  gsl_odeiv_step *s;
  gsl_odeiv_control *c;
  gsl_odeiv_evolve *e;
  gsl_odeiv_system sys;
  static int func(double t, const double y[], double f[], void *params)
{ matrix &A=*reinterpret_cast<matrix *>(params);
    for (int i=0, i_end=A.size(); i<i_end; ++i) {
      f[i]=0;
      for (int j=0, j_end=A.size(); j<j_end; ++j)
        f[i]+=A(i, j)*y[j];
    }
    return GSL_SUCCESS;
  }
public:
  // use a const reference to ensure that stepper does not change
matrix A stepper(const matrix &A) {
    const gsl_odeiv_step_type *T=gsl_odeiv_step_rkf45;
    s=gsl_odeiv_step_alloc(T, A.size());
    c=gsl_odeiv_control_y_new(1e-6, 0.0);
    e=gsl_odeiv_evolve_alloc(A.size());
    sys.function=func;
    sys.jacobian=0;
    sys.dimension=A.size();
    // must cast constness away because of the strong C++ type system
    // otherwise code will not compile
    sys.params=const_cast<matrix *>(&A);
  }
  ~stepper() {
    gsl_odeiv_evolve_free(e);
    gsl_odeiv_control_free(c);
    gsl_odeiv_step_free(s);
  }
  int operator()(double dt, double &t, double t_end, vector &x) {
    return gsl_odeiv_evolve_apply(e, c, s,
                                  &sys, 
                                  &t, t_end, &dt, &x(0));
  }
};
     
int main() {
  // system size
  int N=10;
  // fill coefficient matrix
  matrix A(N);
  for (int i=0; i<N; ++i)
    for (int j=0; j<N; ++j) {
      A(i, j)=0;
      if (i==j)
        A(i, j)=1;
    }
  // initial condition
  vector x(N);
  for (int i=0; i<N; ++i)
    x(i)=1;
  // solve system
  stepper step(A);
  double t=0.0, t_end=2.0, dt=1e-4;
  while (t<t_end) {
    int status=step(dt, t, t_end, x);
    if (status!=GSL_SUCCESS)
      break;
    std::cout << t;
    for (int i=0; i<N; ++i)
      std::cout << '\t' << x(i);
    std::cout << '\n';
  }
  return 0;
}


        Heiko 

-- 
-- Es gibt keine reine Wahrheit, aber ebensowenig einen reinen Irrtum.
-- (Friedrich Hebbel, dt. Dichter, 1813-1863)
-- Cluster Computing @ http://www.clustercomputing.de
--       Heiko Bauke @ http://www.mpi-hd.mpg.de/personalhomes/bauke
-- 
-- Der Pragmatiker entscheidet Fälle nicht nach Grundsätzen, sondern
-- fallweise. (Ron Kritzfeld)
-- Cluster Computing @ http://www.clustercomputing.de
--       Heiko Bauke @ http://www.mpi-hd.mpg.de/personalhomes/bauke