Eigen, How to handle element-wise division by zero-CodePudding

How to handle an element-wise division in Eigen if some divider could be 0 ?
I would like the result to be 0.

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

int main()
{
  ArrayXd a(3), b(3), c;
  a << 1, 2, 0;
  b << 1.2, 3.4, 5.6;
  c = b / a;
}

CodePudding user response：

You can handle division by zero by using a binaryExpr:

ArrayXd a(3), b(3);
a << 1, 2, 0;
b << 1.2, 3.4, 5.6;

Eigen::ArrayXd c = a.binaryExpr(b, [](auto x, auto y) { return y==0 ? 0 : x/y; });

CodePudding user response：

I have found this solution using binaryExpr but not sure it's efficient. Need to test the speed. I have tens of thousands of values to compute.

#include <iostream>
#include <Eigen/Dense>

using namespace Eigen;

template<typename Scalar> 
struct SecureDivisionOp {
  SecureDivisionOp(Scalar v_ = 0) { v = v_; }
  Scalar v;
  Scalar operator()(const Scalar& a, const Scalar& b) const { 
    return b==0 ? v : a/b; 
  }
};

int main()
{
  ArrayXd a(3), b(3), c;
  a << 1, 2, 0;
  b << 1.2, 3.4, 5.6;

  c = b.binaryExpr(a, SecureDivisionOp<double>(0));
  std::cout << c << std::endl;
}

CodePudding user response：

If your compiler is reasonably good with auto-vectorization, you can use the .select( ) mechanism of Eigen:

c = (a!=0).select(b/a, 0.0);

Clang 6 or newer (with -O3 -DNDEBUG -march=skylake) will compile the main-loop of that to something equivalent to:

        vmovupd ymm1, ymmword ptr [rax]
        vcmpneqpd       ymm2, ymm1, ymm0
        vmaskmovpd      ymm3, ymm2, ymmword ptr [rdx]
        vdivpd  ymm1, ymm3, ymm1
        vandpd  ymm1, ymm2, ymm1
        vmovupd ymmword ptr [rax], ymm1

The only slight performance improvement would be to avoid the vmaskmovpd -- but this is necessary from the compiler's viewpoint, since strictly speaking, b[i] is not read for every i where a[i]==0.