Using std::array as state variable

[vchuravy]

I have been using rodent on a project for a while now and it is a wonderful lightweight tool. The reason that I started using it was that I needed a replacement of odeint (I didn't want to pull in all the boost dependencies) and so far it has exceeded my expectations.

So I am solving a multidimensional PDE and I am storing the state in a std::vector. But size the size is constant and known I was wondering if using std::array would be more appropriate.

struct model {
  void operator()( double /* t */, const state &x , state &dxdt) const;
  int size() const { return DX*DX*NM; }
};

But when I change:

typedef std::vector< double > state;

to

typedef std::array< double, DX*DX*NM> state;

I run into the problem that rodent expect an constructor of the state with state(int) and array doesn't provide such a constructor since it should be know at compile time.

In file included from /home/wallnuss/work/anisotropic/Anisotropic1/main.cpp:17:0:
/home/wallnuss/work/anisotropic/Anisotropic1/lib/rodent/rodent/explicitrk.hpp: In instantiation of ‘rodent::AdaptiveRKCashKarp<T_Func, vecT, vecT_traits>::AdaptiveRKCashKarp(T_Func&) [with T_Func = anisotropic_model; vecT = std::array<double, 5120ul>; vecT_traits = rodent::vec_traits<std::array<double, 5120ul> >]’:
/home/wallnuss/work/anisotropic/Anisotropic1/main.cpp:112:73:   required from here
/home/wallnuss/work/anisotropic/Anisotropic1/lib/rodent/rodent/explicitrk.hpp:614:14: error: no matching function for call to ‘std::array<double, 5120ul>::array(const int&)’
       func(_f)
              ^
/home/wallnuss/work/anisotropic/Anisotropic1/lib/rodent/rodent/explicitrk.hpp:614:14: note: candidates are:
In file included from /home/wallnuss/work/anisotropic/Anisotropic1/model.hpp:3:0,
                 from /home/wallnuss/work/anisotropic/Anisotropic1/main.hpp:2,
                 from /home/wallnuss/work/anisotropic/Anisotropic1/main.cpp:14:
/usr/include/c++/4.9.2/array:81:12: note: std::array<double, 5120ul>::array()
     struct array
            ^
/usr/include/c++/4.9.2/array:81:12: note:   candidate expects 0 arguments, 1 provided
/usr/include/c++/4.9.2/array:81:12: note: constexpr std::array<double, 5120ul>::array(const std::array<double, 5120ul>&)
/usr/include/c++/4.9.2/array:81:12: note:   no known conversion for argument 1 from ‘const int’ to ‘const std::array<double, 5120ul>&’
/usr/include/c++/4.9.2/array:81:12: note: constexpr std::array<double, 5120ul>::array(std::array<double, 5120ul>&&)
/usr/include/c++/4.9.2/array:81:12: note:   no known conversion for argument 1 from ‘const int’ to ‘std::array<double, 5120ul>&&’

[jeanluct]

Hi Valentin,

Glad you're finding rodent useful, that's good to hear. What you ask is a reasonable feature. Probably it could be implemented using traits. Right now there is traits.hpp and jlt_traits.hpp, which provide generic traits for the STL and for my own vector/matrix library. The trick would be to add a constructor to the traits somehow, which would default to ignoring the size argument for fixed-size. Maybe... I'm not 100% sure this could be done elegantly (constructors are not things that can be easily aliased, some trickiness might be required). At one point there was support for Blitz++ but they changed the interface so much I couldn't get it to work anymore (there is a devel branch for that).

Definitely it is the kind of thing that you should fork then make a branch for, if you want to try it. It will likely completely break things for a while. There are quite a few things about rodent that I'd like to improve, given its age. Right now I'm mostly maintaining it so it will compile with modern compilers, or to fix bugs. But I'd be happy to help if I can, I'm just not sure I have the time to implement this now.

Best,

Jean-Luc

[vchuravy]

Hi Jean-Luc,

So what is the performance difference (if any) between using std::vector and your matrix library (jlt right?)?

Thanks for the pointers I will look into this if I find some time later this week.

Valentin

[jeanluct]

So what is the performance difference (if any) between using std::vector and your matrix library (jlt right?)?

None, really: jlt::mathvector is really to provide convenient math operations. It is derived from std::vector.

[vchuravy]

So this will apparently work

struct state : std::array<double, DX*DX*NM>
{
  state(){}
  state(const size_t n){}
};

Even thou I am not sure that is a valid abuse of the inheritance system...

[jeanluct]

Do you mean it works with rodent? Sounds perfectly valid to me!

[vchuravy]

Yeah this works with rodent :) doesn't lead to any speed ups though. I was hoping that the compiler would be a bit smarter about access patterns.

[jeanluct]

Sounds like these vectors are pretty long, so the advantage of fixed length is minimal. Worth a try though.

[vchuravy]

Yeah I am solving a PDE consisting of 5120 variables. Rodent itself is very efficent most of the time is spend in the model evaluation so I am trying to speed that up. Do you have any recommendation what kind of solver would be best to minimize the number of times the model is evaluated? Currently I am using AdaptiveRKCashKarp.

[jeanluct]

What kind of equation? It really depends on whether you're accuracy-limit or stiff.

On Mon, Mar 30, 2015 at 10:04 AM, Valentin Churavy <notifications@github.com

wrote:

Yeah I am solving a PDE consisting of 5120 variables. Rodent itself is very efficent most of the time is spend in the model evaluation so I am trying to speed that up. Do you have any recommendation what kind of solver would be best to minimize the number of times the model is evaluated? Currently I am using AdaptiveRKCashKarp.

— Reply to this email directly or view it on GitHub https://github.com/jeanluct/rodent/issues/5#issuecomment-87602677.

[vchuravy]

It is a reaction-diffusion model and it certainly has stiff properties.

[jeanluct]

Sorry about the long delay -- my only suggestion then is to try an implicit solver, so you can take larger time steps. But this only really helps if you can take the accuracy hit.

[vchuravy]

Thanks for the response :) we will see what is going to work.