8

I am writing some code to integrate ODE's. This question is as much a request for coding advice as for a solution, so if you have an alternative suggestion to the one I am about to offer please let me know!

The "objects" to be integrated by the ODE integrator come in "blocks" of 6... The reason for this is that I have a std::vector of doubles, and they are arranged in the following way:

The first 3 doubles are position coordinates; x, y and z. The next 3 doubles are velocity coordinates; x, y and z.

So, now you know that, I have a function which takes pairs of "position" """vectors""" as arguments and returns some sort of result... See where I'm going with this?

Currently the function expects 2 lots of position coordinates in the following manner:

std::vector<double> magic_gravity_formula(const std::vector<double> &r1,
          const std::vector<double> &r2, const double m1, const double m2)

I don't want to copy all the data in group of 3 to new vectors - that's a crazy (and very slow) way to program something.

I could use pointers to the raw data instead... and just pass a pointer to the x coordinate (first item in a block of 3 doubles) - this seems kind of okay to me but perhaps there's a better method? Kind of something like Python or Matlab array slicing? Can I do something like that?

I kind of want to pass a new vector (or some sort of wrapper class?), created from the data already stored in the array... Something kind of like

std::vector<double> sub_section_of_data = data[0..2] // Obviously pseudocode!

Okay so the above is nonsensical, because presumably a language which implemented that syntax would still make a copy operation, which is likely to be slow - exactly what I am trying to avoid...

So yeah I am not sure exactly the best way to proceed here - can anyone suggest a "good" solution? (In the non-subjective way!)

EDIT: To make this really clear - the problem is I don't want to do something like:

std::vector<double> r1_temp;
r1_temp.push_back(data[0]); // Copy ! Bad !
r1_temp.push_back(data[1]);
r1_temp.push_back(data[2]);

... same for an r2 ...

std::vector<double> force = magic_gravity_formula(r1, r2, m1, m2);

EDIT 2: Consider compiler options - would the compiler optimize my code for me by doing something like changing the function to accept arguments in the following way:

std::vector<double> super_gravity_formula(double x1, double y1, double z1, double x2, double y2, double z2, double m1, double m2)

In which case, perhaps this question isn't important? (Other than form a "make your code look nice to read" point of view.)

Edit 3: Therefore it is still important.

Drax
  • 12,682
  • 7
  • 45
  • 85
FreelanceConsultant
  • 13,167
  • 27
  • 115
  • 225
  • Have you looked into using move semantics? – shuttle87 Jun 08 '15 at 20:01
  • I have not, what are these things of which you speak? – FreelanceConsultant Jun 08 '15 at 20:05
  • 1
    *I could use pointers to the raw data instead.* - That sounds like a plan. – Captain Giraffe Jun 08 '15 at 20:06
  • I think I misunderstood the use case here, just to clarify: you have one `std::vector` and you want to use use the data from that to call the function with the signature `magic_gravity_formula(const std::vector &r1, const std::vector &r2, const double m1, const double m2)`? If that's the case then I think Yakk's answer does exactly what you want. – shuttle87 Jun 08 '15 at 20:12
  • 1
    Essentially the function needs to use the data from the "big vector with all the data in it" - but these things must be passed in blocks of three objects, it's no good, for example, to pass the *entire array as a reference* - that's not going to work! – FreelanceConsultant Jun 08 '15 at 20:13
  • Slicing 1D array can cause significant cache performance problem. Be aware. Also take a look at std::valarray – user3528438 Jun 08 '15 at 20:14
  • @user3728501 Whats wrong with a begin and an end pointer? – Captain Giraffe Jun 08 '15 at 20:15
  • 2
    I think you should profile before deciding what's faster. If your vectors are only 6 elements in size copying may prove to be faster than potentially invalidating parts of the `CPU` cache by following pointers. – Galik Jun 08 '15 at 20:17

4 Answers4

8

You want a view into an existing vector.

std::experimental::array_view<T>, or roll your own.

An array view is a pair of T*, and accessors that let you treat it sort of like an array operator[] .begin() .back() size() empty() etc.

Here is one of many such implementations I've rolled. That one is a bit heavy weight, with a range<Iterator> view that adaptively handles random access iterators, and an array_view<T> that inherits from range<T*>.

If it doesn't work, search for another post on SO by "yakk" with the word "array_view" or debug it yourself. I've written it at least a half dozen times, with more or less debugging, with different const correctness rules.

Once you have it, { vec.data()+index, vec.data()+index+3 } will construct an array_view<double> with next to zero overhead.

Community
  • 1
  • 1
Yakk - Adam Nevraumont
  • 262,606
  • 27
  • 330
  • 524
  • Do you know whether the {first,last} initializer convention agrees with proposed `std::string_view` and friends? (`boost::string_ref` seems to have gone with `{begin, length}`) – sehe Jun 08 '15 at 20:05
  • Hmm that's a complicated set of code you've got there - I will try to understand it – FreelanceConsultant Jun 08 '15 at 20:18
  • @sehe The `array_view` I use in my day job supports `{first, length}` *and* `{first, last}`, because I find both useful. `first, length` is awesome when working with legacy code, and `first, last` lets me do sub-views much easier. – Yakk - Adam Nevraumont Jun 08 '15 at 20:36
  • @T.C. Sure, but what is the practical difference between official experimental and unofficial? A very slight degree of increased version stability, and possible you-already-have-it-on-your-system support for voted in stuff? Anyhow, here is a proposal that includes array_view, I don't know if it is the most recent one: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2015/n4512.html – Yakk - Adam Nevraumont Jun 08 '15 at 20:39
  • @user as noted, that is an overly complex one. Really, you just need `T* b=0; T* e=0;` then write some easy accessors that you need. I like being able to implicitly convert from array-like constructs as well, but that is a lot of boilerplate. – Yakk - Adam Nevraumont Jun 08 '15 at 20:41
  • The user has vectors, math sense, of 6 elements. Rather than optimizing the copy of a few bytes I would remove the `std::vector` altogether (dynamic allocation) and would provide proper types to make the code more resilient to mistakes (provide a view into a mix of position and velocity coordinates?)... just saying... I find the solution far more complex than the problem. – David Rodríguez - dribeas Jun 08 '15 at 21:44
1

If what you want is vectors, in the math sense, of 6 elements of type double, std::vector<double> might not be the most efficient representation. I would just go ahead and define a 3D point and then the position + velocity as:

struct point {
   double x, y, z;
};
struct position_and_velocity {
   point pos;
   point vel;
};

Operating on this types will probably be more efficient than operating on vectors (no dynamic allocations), and you could just operate on this as values. A whole copy of this structure is probably cheaper than inserting the first element into the std::vector, and a bit slower but not much than providing two pointers into an array to get a slice.

Additionally it will be less error prone, as in the case of slicing an array you could by mistake use the wrong values and end up with things like (pseudo code) slice[1..3] that is two of the position coordinates with one of the velocity coordinates (probably non-sensical).

The function signature is part of the binary interface of your program, a compiler will not change the function signature unless it inlines the code. Furthermore, I would very much doubt that a compiler would remove the dynamic allocation inherent to vectors in any of the transformations (unless it can remove the whole variable).

David Rodríguez - dribeas
  • 204,818
  • 23
  • 294
  • 489
  • It would be a good solution, but unfortunately I have ODE's for other things which operate on single variables and vectors in spaces with more or less dimension than 3 (or 6)... so this method could be used, but would require me to "fill" unused data in these structs with junk... which might not be safe and would make the code confusing for sure! – FreelanceConsultant Jun 09 '15 at 10:47
0

Since the code requires const std::vector<double> & you cannot be smarter.

Otherwise (if the library functions were more generic) you could easily pass a range (see e.g. Boost Range, but you could easily one-off a vector subrange much like boost::string_ref or std::string_view)

sehe
  • 374,641
  • 47
  • 450
  • 633
  • "const double * const" then it is? – FreelanceConsultant Jun 08 '15 at 20:06
  • 1
    Wait. Can you change that interface? (I assumed `magic_gravity_formula` was ... magic and you couldn't change it) – sehe Jun 08 '15 at 20:07
  • 2
    And to answer the question: no, never drop to raw pointers voluntarily.See @Yakk's answer. Use tiny wrappers to encode the semantics (which can then be seen and checked). This will improve your life and the compiler optimizes it away – sehe Jun 08 '15 at 20:09
0

'I could use pointers to the raw data instead... and just pass a pointer to the x coordinate'.

You answered your question, yourself. std::algorithms's do nothing, else.