help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

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


From: Brian Gough
Subject: Re: [Help-gsl] Re: gsl ode solver - const a problem?
Date: Tue, 22 Sep 2009 14:49:00 +0100
User-agent: Wanderlust/2.14.0 (Africa) Emacs/22.2 Mule/5.0 (SAKAKI)

At Mon, 21 Sep 2009 14:16:30 +0200,
Heiko Bauke wrote:
> 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. 

The more realistic case is interfacing to scripting languages, where the
function to be evaluated is defined in the scripting language and passed
through the params argument.  These are generally non-const objects, for
example in guile we would need to do scm_call_0 (SCM proc)


> 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. 

Can you put the const parameter inside a non-const struct to avoid the
cast? (possibly there is some more elegant way to express that in C++).
Here is what I mean:

#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;
  }
};

typedef struct { const matrix * matrix_param; } matrix_param_t;

class stepper {
  gsl_odeiv_step *s;
  gsl_odeiv_control *c;
  gsl_odeiv_evolve *e;
  gsl_odeiv_system sys;
  matrix_param_t p;

  static int func(double t, const double y[], double f[], void *params)
  { const matrix &A= *((matrix_param_t *)(params))->matrix_param;
    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
    p.matrix_param = &A;
    sys.params=&p;
  };

  ~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;
}




reply via email to

[Prev in Thread] Current Thread [Next in Thread]