We can improve this code.

1: First of all, I urge you tochange your underlying data structure.double *** G; hurts my eyes. Consider using for examplestd::vector<std::vector<std::vector<double> > > instead. (If you are using C++11, you can drop the spaces. Yay!) Usingstd::vector has a plethora of advantages. The two most important are that it will be easier to get the code to work correctly, the other is readability.

If for some masochistic reason you don't want to usestd::vector or a similar type, at least use arrays instead of raw pointers.

2:L,M,N andG are horrific variable names. Consider using for examplewidth,length,depth andmatrix -- or something along those lines -- instead. The same goes for your getter functions.

3: Your output stream operator (operator<<) doesn't use any private members ofGrid3d (which is good!), so it should not be declared as afriend of the class.

4: You can safely usestd::swap() directly on your object.

5: Your class is usingdestructive assignment. Thiscould be okay if it's a well thought-out design decision. In almost all cases, however, it is better to make a (deep) copy of your object.

6: Your class isnot exception-safe. Much of the reason for this is that it fiddles with raw pointers. In C++, we try to avoid pointers to owned resources whenever we can. Putting yourdoubles into astd::vector will take care of a lot of this.

7: Use more whitespace.for (i=0;i<L;i++){ is not very easy on the eye.for (i = 0; i < L; i++) { is much better.

8:You probably want touseoperator[]instead ofoperator()for indexing operations. When in doubt, do as the built-ins do, and arrays usefoo[42] for indexing. Update: For multiple indices, it's perfectly okay to useoperator().

9: Use as few#includes as possible in your header file. Move what you can into your implementation file. After removing thefriend declaration (as discussed in point 3), you can moveall your#includes to your implementation file.

10: Don't ever haveusing statements orusing directives in a header file. Any file#includeing your header will have its global namespace polluted. You don't want that. Use full namespace specifiers in header files. You are already doing this, so theusing directive can safely be removed.

There are some other issues, but these are the most important ones. They should give you something to work on. I suggest you create a new thread with your improved code after incorporating these changes, for further feedback.

• 1
  \$\begingroup\$Agreed, except for8. It's fairly common practice to overloadoperator() when you need to index something with multiple indices.\$\endgroup\$

  Yuushi

  – Yuushi

  2013-06-20 01:03:01 +00:00

  CommentedJun 20, 2013 at 1:03
• \$\begingroup\$@Yuushi Oh, right, multiple indices. You are completely correct, I'll edit my answer.\$\endgroup\$

  Lstor

  – Lstor

  2013-06-20 01:05:03 +00:00

  CommentedJun 20, 2013 at 1:05
• \$\begingroup\$There can be efficiency considerations with a contiguous memory block. So potentially std::vector<double> to do the memory management and then do the appropriate calculations inoperator() oroperator[]\$\endgroup\$

  Loki Astari

  – Loki Astari

  2013-06-22 18:19:11 +00:00

  CommentedJun 22, 2013 at 18:19

I have a few things to add to Lstor's review (and a few things to repeat).

double*** makes me cry. That's just waiting for some weird exception safety problem or a memory leak. As Lstor said, use a vector (thoug technicallystd::vector< std::vector< std::vector<double> > >, notstd::vector<std::vector<double> > as he said).

Or, if you really want to use raw allocation, at least usestd::unique_ptr so you have exception safety.

Fun bonus: if you use vectors, you actually won't need your swap (which you don't need now), destructor oroperator= as the default ones will work.

Your dimensions should bestd::size_ts rather thanints.

It doesn't apply quite as much here since meaningful names are rather hard to come up with, but in general, you should always use named parameters in your class declarations.

Grid3d(int,int,int); doesn't exactly tell me much about what I'm passing to the constructor.Grid3d(int L, int M, int N); tells me more.Grid3d(int width, int height, int depth); tells me a lot more (though it brings up a rather difficult problem of what you would name a fourth dimension if you were to have one)

clean is a bad name. I wouldn't expect clean to mean set everything toNAN. I would name itreset ornullify or something. Or, better yet, have a method calledsetAll(double val) that sets the entire matrix. That way the ambiguity of what you're setting it to is gone (and you don't have to duplicate your constructor loop and your clear loops).

Declare variables as tightly scoped as possible.

void Grid3d::clean(){    for (int i=0; i < L; i++){        for (int j=0; j < M; j++){            for (int k=0; k < N; k++){                G[i][j][k] = NAN;            }           }    }}

(And once again, those ints should bestd::size_t.)

You should have a const getter too. Not having one destroys all possibilities of const correctness with usage of your class. For example, it would be nice ifoperator<< took a constGrid3d instead of a mutable one. Or what aboutcurl_r? If I saw the declaration of that, I would wonder why in the world calculating the curl changes the vectors.

T v = ...;return v;

Should usually just bereturn ...;

But in this case, your statement should probably be broken into more digestable pieces:

const double part1 = ...; //With a better name than part1, ideallyconst double part2 = ...;const double part2 = ...;return (part1 * part2) + part3;

(This in no way resembles your formula, but I have depressingly little memory of what the different parts mean, so I have no idea how to meaningfully break it apart).

I like to declare variables that I don't plan on changing asconst. That way, 4 lines later, when I decide I want to reusewidth (or whatever), I get an error instead of then having the wrong value 2 more lines down when I think I'm using the original width again. Then again, there is definitely a potential to over do it :).

This one is mostly personal preference, but I find spaced out arithmetic operators a million times easier to read.4+3*2 vs4 + 3 * 2 (though I would also parenthesize that).

(This one is 80% personal preference.)

I find middle-reference style awkward. I like eitherT& t orT &t (personally, I hugely preferT& t as I see the reference being part of the type, not part of the variable).

std::ostream & operator << (std::ostream & os , Grid3d & g) Is weird to read for me. I would use:std::ostream& operator<<(std::ostream& os , Grid3d& g).

(And I would actually usestd::ostream& operator<<(std::ostream& os , const Grid3d& g).

Being a bit pedantic,clean could use the class' getter rather than direct access to G. That way if you change something about G, you only have to alter it in one place.

(If you're worried about performance,get will likely be inlined anyway meaning it will have no overhead. You would of course want to profile though if performance is an issue.)

• 3
  \$\begingroup\$inline is basically a compiler suggestion. The compiler is free to do its own cost/benefit analysis and completely ignore it. For something as simple as a getter/setter, however, on any kind of reasonable optimization level, it's almost 100% guaranteed to be inlined automatically anyway.\$\endgroup\$

  Yuushi

  – Yuushi

  2013-06-20 10:48:40 +00:00

  CommentedJun 20, 2013 at 10:48
• 1
  \$\begingroup\$+1. Good points! Also, like @Yuushi points out,inline is just a hint to the compiler, which is free to do what it wants. Defining a function in the class definition is the same as declaring itinline.\$\endgroup\$

  Lstor

  – Lstor

  2013-06-20 18:05:56 +00:00

  CommentedJun 20, 2013 at 18:05
• \$\begingroup\$Ah, seems I need to do some reading on inline. I've been procrastinating it for years now :).\$\endgroup\$

  Corbin

  – Corbin

  2013-06-20 19:24:53 +00:00

  CommentedJun 20, 2013 at 19:24

Quickie style issues

Several of these were mentioned in the other answers and some of them weren't.

• You should generally not write code in the global namespace. There are a number of good reasons for this, but this margin is too small to contain them.The programmers stackexchange has some of them.

• You usually shouldn't callstd::swap directly, though here it is harmless. You should get in the habit of swapping like this:using std::swap; swap(foo, bar); It won't cause a differentswap to be used here, but you shouldn't try to keep in your head which namespaceswap comes from for each type you use (seethis stackoverflow question for more discussion onswap and argument-dependent lookup). I discussswap a bit more later.

• Dimensions should generally be passed assize_t, notint.

• You should not putusing statements in a header unless inside a function (or sometimes class) definition. You shouldespecially not putusing statements in the global namespace in a header (and usually not in a source file -- putusing statements in a namespace to prevent symbol clashes and ODR violations during linking).

• You should name the arguments in your header file. Names are the primary form of API documentation and you've got none of that at all in this header file. As LStor mentioned, the names you do have are not great, especiallyL,get_L, etc.

• What isclean supposed to be doing? Are you trying to pollute your surrounding data with NANs in order to detect use of an uninitialized grid? Considerclear instead with 0s instead of NANs, or maybesetAll as Corbin suggests.

Underlying data structure

I (mostly) disagree with the other answers aboutdouble***. This is exactly the sort of class that should use a contiguous array as its underlying data structure. You get exception safety by deleting the memory in your destructor (that, after all, is all that smart pointers do!). However, you're losing any efficiency gain of using the raw data structure by creating an array for each row. Instead, use adouble* as your underlying memory like so:

class Grid3d {  private:    size_t L, M, N;    double* buffer;  public:    Grid3d(size_t L, size_t M, size_t N)      : L(L), M(M), N(N), buffer(new double[L * M * N]) {}    ~Grid3d() { delete [] buffer; }    /* ... */};

This has the following advantages over yourdouble*** approach and over thevector-based approaches suggested in the other answers.

• The buffer is allocated with a single allocation, which:

  • reduces heap fragmentation
  • reduces the number of system calls (so is faster)
  • eliminates the possibility of partial allocation, which genuinely was an exception-safety issue with the previous implementation.

• Element accesses will require only one memory lookup, instead of 3.

This approach is in line withHerb Sutter's excellent guideline:

Don’t use explicit new, delete, and owning * pointers, except in rare cases encapsulated inside the implementation of a low-level data structure.

This is one of those low-level data structures.

Allow your class to be used correctly when const.

It is a good general guideline to declare variables const whenever possible. If you didn't wish to make your class assignable and swappable, I'd suggest making its dimensions const. As it is, your users have no way to access the data in the grid when theGrid3d is passed to functions by constant reference, for example (as theGrid3d arguments should be passed toSystem::curl_r in your example). Here's a solution:

Add a private function to share the index calculation logic:

size_t compute_index(size_t i, size_t j, size_t k) const {    return i * (M * N) + j * N + k;  // untested}

Replaceget(int, int, int) byget(int, int, int) const:

double Grid3d::get(size_t i, size_t j, size_t k) const {    /* add bounds checking */    return buffer[compute_index(i, j, k)];}

Addconst and non-const versions ofoperator():

double& Grid3d::operator()(size_t i, size_t j, size_t k) {    return buffer[compute_index(i, j, k)];}double Grid3d::operator()(size_t i, size_t j, size_t k) const {    return buffer[compute_index(i, j, k)];}

Swap correctness

std::swap actually correctly swaps your type (as LStor points out) because its member variables are all built-in types (this is sufficient but not necessary). You needn't provide aswap method at all, but if you do you can write it as

void Grid3d::swap(Grid3d& other) {  using std::swap;  swap(*this, other);}

Using theusing std::swap; construction as I have is unnecessary because you know that you actually will be usingstd::swap, not doing ADL. However, as mentioned above, it's good to get in the habit of doing this because it's absolutely necessary when the arguments toswap are of unknown type, as when writing a template (begin andend are the other two functions which merit this treatment).

The assignment operator

You'll want to provide two assignment operators:

Grid3d& Grid3d::operator=(Grid3d&& other) {  swap(other);  return *this;}Grid3d& Grid3d::operator=(Grid3d other) {  swap(other);  return *this;}

The move-assignment operator prevents an unnecessary copy in expressions likegrid1 = std::move(grid2).

There are, of course, other areas of improvement, but this should be enough to get you started.

• \$\begingroup\$Many good points. However, except if writing an industrial-strength library for numeric computations, your suggestion about adouble* is simply premature optimization.\$\endgroup\$

  Lstor

  – Lstor

  2013-06-20 18:26:29 +00:00

  CommentedJun 20, 2013 at 18:26
• 2
  \$\begingroup\$You're right that adouble* would have much better performance than adouble***. One of us probably should have said that sooner. 3 levels of indirection is going to murder the cache, not to mention all the load/store. It could still easily be a std::vector<double> instead of double* though. And, this isnot a low level data structure. Using a single dimensional vector is going to offer the same performance as a single dimensional array on all but the worst compilers.\$\endgroup\$

  Corbin

  – Corbin

  2013-06-20 19:28:46 +00:00

  CommentedJun 20, 2013 at 19:28
• \$\begingroup\$His assignment operator is fine by the way. It's just the classic copy and swap idiom (which is typically done as a simple method of getting a strong exception guarantee -- something your operator= does not offer).\$\endgroup\$

  Corbin

  – Corbin

  2013-06-20 19:37:40 +00:00

  CommentedJun 20, 2013 at 19:37
• \$\begingroup\$@Corbin With a double*** that gets deleted by the constructor, it is absolutely not fine. It goes like this. Start with the code Grid3d a(10, 10, 10); Grid3d b; b = a; When b = a is computed, a is copied (and in particular, a.G is copied). Now the copy and b are swapped. b.G == a.G. When a and b go out of scope, they'll both delete G. However, if Grid3d were implemented with a vector, you'd be right; copy and swap would be fine.\$\endgroup\$

  ruds

  – ruds

  2013-06-20 19:57:19 +00:00

  CommentedJun 20, 2013 at 19:57
• \$\begingroup\$@Lstor It's really not a premature optimization. You should not sacrifice code clarity in return for unnecessary optimization, that is true. But a vector<vector<vector<double>>> is actually far less clear than a double*, and insisting on it is premature pessimization. Whether you use a vector or a double* here is a matter of taste, I suppose, but it's not completely clear either way. After all, this buffer is not logically a vector.\$\endgroup\$

  ruds

  – ruds

  2013-06-20 20:01:14 +00:00

  CommentedJun 20, 2013 at 20:01

A few points mostly regarding efficiency:

• As @ruds said, use a contiguous array. I would use avector<double> instead of thedouble * he suggests, as this simplifies many member functions, at a trivial extra overhead.
• Try to decouple thesize of the leading dimension of your array (i.e. the one whose elements you store contiguously) from itsallocation. This allows you to e.g. store 35*10*10 array in a structure with 36*10*10 elements. If you want to make use of the vector processing units (SSE3, Velocity Engine, etc.) that many modern processors have, having structures aligned to 2- or 4-element boundaries can bring tremendous performance benefits. For this reason, numeric libraries like BLAS and LAPACK tend to have separate parameters for the size and allocation of the leading dimension.
• If your problem is suitable forfloat instead ofdouble, that is also a tremendous performance win.

• c++
• matrix