Округление чисел ниже определенного порога до нуля в Eigen - PullRequest
Новая
       25

Округление чисел ниже определенного порога до нуля в Eigen

0 голосов
/ 03 февраля 2019

Я широко использую Eigen в научном приложении, которое я разрабатывал в течение некоторого времени.Поскольку я реализую численный метод, число ниже определенного порога (например, 1e-15) не представляет интереса, и это замедляет вычисления и увеличивает частоту появления ошибок.

Следовательно, я хочу округлитьотключить номера ниже этого порога до 0.Я могу сделать это с помощью цикла for, но использование нескольких относительно больших матриц (2М ячеек и выше на матрицу) с помощью цикла for - if стоит дорого и замедляет меня, поскольку мне нужно делать это несколько раз.

Есть ли более эффективный способ сделать это с библиотекой Eigen?

Другими словами, я пытаюсь исключить числа ниже определенного порога в моем конвейере вычислений.

Ответы [ 2 ]

0 голосов
/ 03 февраля 2019

Самый короткий способ написать то, что вы хотите, это 

void foo(Eigen::VectorXf& inout, float threshold)
{
    inout = (threshold < inout.array().abs()).select(inout, 0.0f);
}

Однако ни сравнения, ни метод select не векторизованы Eigen ( на данный момент ).

Если скорость важна, вам нужно либо написать некоторый код SIMD вручную, либо написать собственный функтор, который поддерживает метод packet (он использует внутреннюю функциональность Eigen, поэтому он не гарантированно стабилен!):

template<typename Scalar> struct threshold_op {
  Scalar threshold;
  threshold_op(const Scalar& value) : threshold(value) {}
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const{
    return threshold < std::abs(a) ? a : Scalar(0); }
  template<typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const {
    using namespace Eigen::internal;
    return pand(pcmp_lt(pset1<Packet>(threshold),pabs(a)), a);
  }
};
namespace Eigen { namespace internal {

template<typename Scalar>
struct functor_traits<threshold_op<Scalar> >
{ enum {
    Cost = 3*NumTraits<Scalar>::AddCost,
    PacketAccess = packet_traits<Scalar>::HasAbs };
};

}}

Затем его можно передать в unaryExpr: 

inout = inout.unaryExpr(threshold_op<float>(threshold));

Godbolt-Demo ( должен работать с SSE / AVX / AVX512 / NEON /)...): https://godbolt.org/z/bslATI

Фактически единственной причиной замедления могут быть субнормальные числа.В этом случае простой 

_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);

должен добиться цели (ср: Почему изменение от 0,1f до 0 замедляет производительность в 10 раз? )

0 голосов
/ 03 февраля 2019

Eigen имеет метод с именем UnaryExpr, который применяет данный указатель функции к каждому коэффициенту в матрице (он также имеет разреженные и массивные варианты).

Будет проверять его производительность иобновите этот ответ соответственно.

...