Hi On 07/22/2014 01:24 AM, Erik Schultheis wrote:

Dear all,

i am simulating particle trajectories, and depending on the data i want to generate i may or may not have to track additional data [boost.ublas c_matrix] along those trajectories. Right now i enable/disable that via a template parameter (so, using std::tuple(vec, vec) or tuple(vec, vec, c_matrix) and fusion_algebra), which has a few drawbacks. The algorithm is already templated to support varying dimensions via a template parameter, so the additional data practically doubles the number of template instantiations for system and all possible observers. This actually got so bad that gcc (mingw gcc 4.8) refuses to compile with debugging enabled. Therefore, i would like to make the switch between tracking the matrix and ignoring it a runtime choice instead of a template parameter.

I don't understand fully what you exactly want to do. Can you show us some more code? Can you put the tracking into an observer? I mean, do you need to solve an ODE to track the additional data?

However, the cost of tracking that matrix is quite high (that is why it is optional), but simply disableing the d matrix / dt calculation in the system does not speed up the calculation much (86 s -> 58 s vs 27 s without matrix). A lot of time is spent in adding and multiplying numbers that are not used anywhere, or actually even copying them (e.g std::swap(c_matrix, c_matrix) shows up quite high in gprof (even though i am not sure how reliable the -pg option is in conjunction with -O3)).

So what i am asking here is for suggestions how to best let odeint treat the c_matrix as uninitialized memory and never perform any operations on that depending on a runtime parameter. Options i thought of were: * deriving a class from c_matrix and providing my own versions of operator overloads that check a flag: as far as i know, boost.ublas matrices use some template magic to let operators only return expression templates to optimize evaluation if operations are chained, and i would not be sure about how to actually retain that behaviour in case the matrix is tracked. * creating my own algebra class: from what i could tell from a few quick glances into the algebra files, this would probably be conceptually easier than the other option, but would require a lot of repetetive coding.

Does anyone have suggestion for a more efficient implementation strategy?

Have you considered to use a fusion view which only points to tuple( vec , vec ). This could be used as the state type you use. Nevertheless, you need to instantiate the stepper with a real resizable state_type which you could be again a tuple( vec , vec ).

From the options you showed here, an own algebra looks most promising. This algebra you also be implemented in terms of the existing fusion_algebra and redirect all calls to fusion_algebra.