3
\$\begingroup\$

I have made a C++ class to manipulate 3D vectors in Cartesian coordinates.

However, the performance of my class is much slower (about 2.5x) than simply using something like double p [3] and then running for loops for additions, subtractions, etc.

My specific concerns are:

  • Is there something fundamentally wrong that I am doing, that is causing this slowdown?
  • Is there a way to achieve what I am trying, but without as much of a slowdown?
  • Is this just a fundamental limitation of classes and operator overloading that I have to live with, in return for the convenience?

Here is the class:

#ifndef COORD_H
#define COORD_H


#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstring>
#include <complex>
#include <vector>
#include <map>
#include <stdlib.h>
#include "math.h"


// Useful reference: http://courses.cms.caltech.edu/cs11/material/cpp/donnie/cpp-ops.html


/*! \brief Class to store and manipulate points or position vectors in 3D, in Cartesian coordinates.*/
class vect3D
{
public:

    double x = 0.0, y = 0.0, z = 0.0;
    double tol = 1.0e-15;


    // Constructors
    vect3D() {};
    vect3D(std::initializer_list<double> RHS) { x = *RHS.begin(); y = *(RHS.begin()+1); z = *(RHS.begin()+2); };
    vect3D(std::vector<double> &RHS) { x = RHS[0]; y = RHS[1]; z = RHS[2]; };
    vect3D(double *RHS) { x = RHS[0]; y = RHS[1]; z = RHS[2]; };


    // Assignment
    vect3D& operator=(const vect3D &RHS) { x = RHS.x; y = RHS.y; z = RHS.z; return *this; };


    // Addition and subtraction
    vect3D& operator+=(const vect3D &RHS) { x += RHS.x; y += RHS.y; z += RHS.z; return *this; };
    vect3D& operator-=(const vect3D &RHS) { x -= RHS.x; y -= RHS.y; z -= RHS.z; return *this; };

    vect3D& operator+=(const double &RHS) { x += RHS; y += RHS; z += RHS; return *this; };
    vect3D& operator-=(const double &RHS) { x -= RHS; y -= RHS; z -= RHS; return *this; };

    vect3D operator+(const vect3D &RHS) { return vect3D(*this) += RHS; };
    vect3D operator-(const vect3D &RHS) { return vect3D(*this) -= RHS; };

    vect3D operator+(const double &RHS) { return vect3D(*this) += RHS; };
    vect3D operator-(const double &RHS) { return vect3D(*this) -= RHS; };


    // Scalar product and division
    vect3D& operator*=(const double &RHS) { x *= RHS; y *= RHS; z *= RHS; return *this; };
    vect3D& operator/=(const double &RHS) { x /= RHS; y /= RHS; z /= RHS; return *this; };

    vect3D operator*(const double &RHS) { return vect3D(*this) *= RHS; };
    vect3D operator/(const double &RHS) { return vect3D(*this) /= RHS; };

    friend vect3D operator*(double c, vect3D &vec) { return vec*c; };
    friend vect3D operator/(double c, vect3D &vec) { return vec/c; };


    // Comparisons
    bool operator==(const vect3D &RHS) { return ((x - RHS.x < x*tol) && (y - RHS.y < y*tol) && (z - RHS.z < z*tol)); };
    bool operator!=(const vect3D &RHS) { return !(*this == RHS); };

    bool operator>=(const vect3D &RHS) { return ((x >= RHS.x) && (y >= RHS.y) && (z >= RHS.z)); };
    bool operator<=(const vect3D &RHS) { return ((x <= RHS.x) && (y <= RHS.y) && (z <= RHS.z)); };
    bool operator>(const vect3D &RHS) { return !(*this <= RHS); };
    bool operator<(const vect3D &RHS) { return !(*this >= RHS); };


    // Euclidean norm
    double norm2() { return std::sqrt(std::pow(x, 2) + std::pow(y, 2) + std::pow(z, 2)); };
    friend double norm2(vect3D const &a) { return std::sqrt(std::pow(a.x, 2) + std::pow(a.y, 2) + std::pow(a.z, 2)); };


    // Dot product
    friend double dot(vect3D const &a, vect3D const &b) { return a.x*b.x + a.y*b.y + a.z*b.z; };

    // Cross product
    friend vect3D cross(vect3D const &a, vect3D const &b) { return {(a.y*b.z - a.z*b.y), (a.z*b.x - a.x*b.z), (a.x*b.y - a.y*b.x)}; };


    // Print to stream
    friend std::ostream& operator<<(std::ostream &stream, vect3D const &p) { return stream << "(" << p.x << ", " << p.y << ", " << p.z << ")" << std::flush; };


    // Function to explicitly return coordinates as an array of doubles, if needed
    void DoubleArray(double *v) { v[0] = x; v[1] = y; v[2] = z; return; };


    // To access coordinates using square bracket notation, for convenience
    std::vector<double *> p = std::vector<double *> {&x, &y, &z};
    double operator [] (int ii) const { return *(p[ii]); };

};

#endif

Thanks!

\$\endgroup\$
  • 2
    \$\begingroup\$ Hey, welcome to Code Review! If your code works but you want it to be better (in terms of speed or maintainability), you have come to the right place. Stack Overflow is a better place if you are looking for why something does not work or how to do something in the first place. \$\endgroup\$ – Graipher Nov 29 '18 at 14:33
  • 1
    \$\begingroup\$ Have you thought about inverting the access? E.g. store std::array<double, 3> and provide x(), y(), z() functions returning a reference to their corresponding elements. \$\endgroup\$ – Incomputable Nov 29 '18 at 20:03
  • \$\begingroup\$ @Incomputable Now that you mention it, I refer to the x, y, z members almost as frequently through [ ] as through .x, .y etc., so I'll try out what you suggest and assess the impact on performance. \$\endgroup\$ – EM_IE Nov 30 '18 at 18:32
4
\$\begingroup\$

This is a general review, not a specific performance review.


Many of the operators ought to be const-qualified. Relational operators are of questionable value, given there isn't a natural ordering for coordinate vectors.


The constructors are fragile because they don't enforce the requirement that the caller provide at least three values. The initializer_list constructor is pointless really - just provide an ordinary 3-argument constructor that will be properly checked at the call site, and prefer initialization to assignment (this last may be a performance issue):

vect3D(double x, double y, double z)
  : x{x}, y{y}, z{z}
{}

I don't think there's any need for the assignment operator, unless you really don't want to copy tol as the compiler-generated one does - in which case, I'd recommend an explanatory comment, because it just looks like you forgot that member.


You seem to missing an include of <cmath> for std::pow() and std::sqrt(). However, you don't want to be using those functions, given that <cmath> provides a (much better behaved and possibly more efficient) std::hypot() function:

// Euclidean norm
double norm2() const { return std::hypot(x, y, z); };
friend double norm2(vect3D const &a) { return a.norm2(); };

(Note: in C++14 and earlier, std::hypot() only took two arguments, so you had to write std::hypot(std::hypot(x, y), z) or similar.)


The extra std::vector of pointers that each vect3D carries around with it is likely to be a performance drag. If you really want that help for indexing, consider a static constexpr array of member-pointers instead (which won't take space in the object, nor need to allocate when it's copied):

    // To access coordinates using square bracket notation, for convenience
    double operator[](int i) const { return this->*p[i]; };
    double& operator[](int i) { return this->*p[i]; };

private:
    static constexpr double vect3D::*p[] = { &vect3D::x, &vect3D::y, &vect3D::z };

Don't stream std::flush when writing output. If callers want output to be flushed, they will do this themselves; if not, they certainly don't want it to be imposed upon them.


Header rationalisation - in conjunction with the other improvements I recommend, we only need

#include <cmath>
#include <ostream>
\$\endgroup\$
  • \$\begingroup\$ These are extremely helpful; thank you. In case it is helpful to others, I want to note that implementing your replacement of std::vector with the static constexpr improved performance by ~2x, so now the class is almost on par with using raw pointers. \$\endgroup\$ – EM_IE Nov 29 '18 at 15:48
  • 1
    \$\begingroup\$ I'm not sure whether you ever give different coordinates different tol values; if not, then making that member a static const would reduce the object size by 25%, probably with a corresponding speed boost. \$\endgroup\$ – Toby Speight Nov 29 '18 at 15:55
  • 1
    \$\begingroup\$ What do you think about this: a free standing template function component, which has component’s name as first argument and vector as second, input to function argument? That would allow crazy extensibility. Here is an example component<vec3d::x>(line3d). One could also write delegate like x_component(line3d) to make it less verbose. \$\endgroup\$ – Incomputable Nov 29 '18 at 20:26
3
\$\begingroup\$

Is this just a fundamental limitation of classes and operator overloading that I have to live with, in return for the convenience?

Not at all! C++ is all about abstractions with zero overhead. The design behind C++ is that it shouldn't leave any place for another language closer to the machine, while providing powerful abstractions to leverage native speed.

Further down a semi-complete review of your code (@TobySpeight already covered much ground).

#ifndef COORD_H
#define COORD_H


#include <iostream>
#include <fstream>
#include <iomanip>

That seems a bit heavy on I/O headers

#include <cstring>
#include <complex>
#include <vector>
#include <map>
#include <stdlib.h>

There's a bit of clean up to do here

//#include "math.h"


// Useful reference: http://courses.cms.caltech.edu/cs11/material/cpp/donnie/cpp-ops.html


/*! \brief Class to store and manipulate points or position vectors in 3D, in Cartesian coordinates.*/
class vect3D

Why not a struct? Since your data is public anyway

{
public:

    double x = 0.0, y = 0.0, z = 0.0;
    double tol = 1.0e-15;

What you're looking for is standardized here https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon

    // Constructors
    vect3D() {};

Why would you define such a useless constructor? Just default it: vect3D() = default

By the way, there's a useless superfluous semi-colon at the end of all your function definitions.

    vect3D(std::initializer_list<double> RHS) { x = *RHS.begin(); y = *(RHS.begin()+1); z = *(RHS.begin()+2); };

If you conceive your vect3D as an aggregate (no user provided constructor, no virtual functions, and some other restrictions) you can use aggregate initialization directly (vect3D v{1., 2., 3.};) without having to define a brace-list constructor

    vect3D(std::vector<double> &RHS) { x = RHS[0]; y = RHS[1]; z = RHS[2]; };

So what if I provide a vector with one or two elements? Or even with six elements: that's clearly a mistake that won't be called out by the compiler.

    vect3D(double *RHS) { x = RHS[0]; y = RHS[1]; z = RHS[2]; };

Same here. The advice is to provide an interface that is hard to misuse

By the way, did you wonder why you got a lot of warnings about a deprecated implicit copy constructor? do you see why? Think of what will happen when the default copy constructor will copy your vector member by member, your vector of pointers included

    // Assignment
    vect3D& operator=(const vect3D &RHS) { x = RHS.x; y = RHS.y; z = RHS.z; return *this; };

    // Addition and subtraction
    vect3D& operator+=(const vect3D &RHS) { x += RHS.x; y += RHS.y; z += RHS.z; return *this; };
    vect3D& operator-=(const vect3D &RHS) { x -= RHS.x; y -= RHS.y; z -= RHS.z; return *this; };

    vect3D& operator+=(const double &RHS) { x += RHS; y += RHS; z += RHS; return *this; };
    vect3D& operator-=(const double &RHS) { x -= RHS; y -= RHS; z -= RHS; return *this; };

    vect3D operator+(const vect3D &RHS) { return vect3D(*this) += RHS; };
    vect3D operator-(const vect3D &RHS) { return vect3D(*this) -= RHS; };

    vect3D operator+(const double &RHS) { return vect3D(*this) += RHS; };
    vect3D operator-(const double &RHS) { return vect3D(*this) -= RHS; };


    // Scalar product and division
    vect3D& operator*=(const double &RHS) { x *= RHS; y *= RHS; z *= RHS; return *this; };
    vect3D& operator/=(const double &RHS) { x /= RHS; y /= RHS; z /= RHS; return *this; };

    vect3D operator*(const double &RHS) { return vect3D(*this) *= RHS; };
    vect3D operator/(const double &RHS) { return vect3D(*this) /= RHS; };

    friend vect3D operator*(double c, vect3D &vec) { return vec*c; };
    friend vect3D operator/(double c, vect3D &vec) { return vec/c; };


    // Comparisons
    bool operator==(const vect3D &RHS) { return ((x - RHS.x < x*tol) && (y - RHS.y < y*tol) && (z - RHS.z < z*tol)); };
    bool operator!=(const vect3D &RHS) { return !(*this == RHS); };

    bool operator>=(const vect3D &RHS) { return ((x >= RHS.x) && (y >= RHS.y) && (z >= RHS.z)); };
    bool operator<=(const vect3D &RHS) { return ((x <= RHS.x) && (y <= RHS.y) && (z <= RHS.z)); };
    bool operator>(const vect3D &RHS) { return !(*this <= RHS); };
    bool operator<(const vect3D &RHS) { return !(*this >= RHS); };


    // Euclidean norm
    double norm2() { return std::sqrt(std::pow(x, 2) + std::pow(y, 2) + std::pow(z, 2)); };

    friend double norm2(vect3D const &a) { return std::sqrt(std::pow(a.x, 2) + std::pow(a.y, 2) + std::pow(a.z, 2)); };

You don't need all those friends since x, y and z are public. That only clutters your class interface.

    // Dot product
    friend double dot(vect3D const &a, vect3D const &b) { return a.x*b.x + a.y*b.y + a.z*b.z; };

    // Cross product
    friend vect3D cross(vect3D const &a, vect3D const &b) { return {(a.y*b.z - a.z*b.y), (a.z*b.x - a.x*b.z), (a.x*b.y - a.y*b.x)}; };


    // Print to stream
    friend std::ostream& operator<<(std::ostream &stream, vect3D const &p) { return stream << "(" << p.x << ", " << p.y << ", " << p.z << ")" << std::flush; };


    // Function to explicitly return coordinates as an array of doubles, if needed
    void DoubleArray(double *v) { v[0] = x; v[1] = y; v[2] = z; return; };

    // To access coordinates using square bracket notation, for convenience
    std::vector<double *> p = std::vector<double *> {&x, &y, &z};

I'm amazed that someone with a sane mind (at least I hope so) would have such an extravagant idea. I'm not sure that square bracket notation is a convenience here -it might as well persuade the user that vect3D behaves like an array or a pointer- but even if you want to provide it, creating a vector of pointers is crazy, without even mentioning what will happen when a user naively calls for the third dimension with vec[3] (and they will). If you really want the bracket operator, implement it with a switch: case 0: return x; ...; default: throw std::out_of_bonds();

    double operator [] (int ii) const { return *(p[ii]); };

};

#endif
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.