I am a complete newbie to C++ and programming in general. I need to write something for scientific purposes and as such, performance is crucial.

I introduced two types, matrices and vectors with intervals as entries, and now want to make them "able to work with", i.e. define basic operations. The type definitions:

#include <eigen3/Eigen/Dense>
#include <boost/numeric/interval.hpp>

namespace bn = boost::numeric;
namespace bi = bn::interval_lib;

// Interval typedefs
using Interval = bn::interval<
double,
bi::policies<
bi::save_state<bi::rounded_transc_std<double> >,
bi::checking_base<double>
>
>;

using Matrix = Eigen::Matrix<Interval, 3, 3>;
using Vector = Eigen::Matrix<Interval, 3, 1>;


Luckily, matrix products work without extra definitions (although I need to call a.lazyProduct(b), otherwise I receive an "operator unambiguous" error, which I don't understand), but inner products and multiplications with constants do not. It seems that I have to manually overload the multiplications, which I did as follows:

Vector operator* (const double& x, const Vector& y)
{
Vector res;
for(int i = 0; i<3;i++) {
res[i] = x*y[i];
}
return res;
}

Matrix operator* (const double& x, const Matrix& y)
{
Matrix res;
for(int i = 0; i<3; i++) {
for(int j = 0; j<3; j++){
res(i,j) = x* y(i,j);
}
}
return res;
}

Interval inner_prod(const Vector& x, const Vector& y){
Interval res(0.0);
for(int i = 0; i<3; i++){
res += x[i]*y[i];
}
return res;
}

Interval inner_prod(const std::vector<double>& x, const Vector& y){
Interval res(0.0);
for(int i = 0; i<3; i++){
res += x[i] * y[i];
}
return res;
}


sizes are hard-coded as I will not be working with any others. I would love to have some feedback and improvement suggestions that focus on performance. Also for the inner product, I was hoping there was a way to make use of std::inner_product for possible better performance, but I can't figure out how.

The boost/numeric/interval.hpp header gives a definition for multiplication of intervals with scalars, I must have made some mistake before that this didn't work. That just leaves my custom types.

• Since you are using a vector of length 3, just use std::array<double, 3> rather than std::vector. That'll save some bytes which store the size in std::vector. – user14717 Jan 28 '18 at 3:50
• I switched to Eigen::Vector3d for the normal double entry vectors (this gives me also the preimplemented methods such as the norm). Still not sure if the straightforward for-loop overloading is really efficient – bernhard_e Jan 29 '18 at 11:31
• It should be fine. I've looked at the asm generated by clang in this situation and it always gets vectorized as long as you use '-march=native -O3' compile flag. – user14717 Jan 29 '18 at 16:20
• #include <eigen3/Eigen/Dense> I think you are supposed to write #include <Eigen/Dense> and compile with a suitable -I flag. – Marc Glisse Feb 4 '18 at 21:52
• The fact that you have to call a.lazyProduct(b) instead of a*b is a sign of a problem, either in Eigen or in Boost (it works fine with the interval type from CGAL), you may want to report that to one of those projects (probably Eigen, since AFAIK Boost.Interval is unmaintained). I wouldn't try to use the combination until that is fixed... – Marc Glisse Feb 4 '18 at 22:04

I'll start by saying I haven't used boost::numeric, so what I say here is based on a quick skim of its documentation and source.

Instead of writing the constant 3 in so much of the code, we should be using rows() or cols() as appropriate (or perhaps even RowsAtCompileTime and ColsAtCompileTime).

Don't pass simple value types by reference - double values should be simply passed by value.

In fact, since Eigen::Matrix has a member that multiplies by a scalar, we can simplify these functions by simply calling that. All we have to do is swap the argument order:

Vector operator*(double x, Vector y)
{
return y *= x;
}

Matrix operator*(double x, const Matrix& y)
{
return y * x;
}


We can also simplify the dot-product function by using the standard inner_product algorithm. We need to know that a vector's start and end iterators are obtained by data() and data() + size() respectively:

#include <numeric>

Interval inner_prod(const Vector& x, const Vector& y)
{
return std::inner_product(x.data(), x.data() + x.size(),
y.data(), Interval{});
}

Interval inner_prod(const std::vector<double>& x, const Vector& y)
{
if (x.size() != y.size())
return {};              // or throw an exception, or something
return std::inner_product(x.begin(), x.end(),
y.data(), Interval{});
}