std::accumulate()只对一个复数std::vector的实部进行累加

Question

我曾经对称（Hermitian）矩阵（该矩阵位于 std::vector 中）求和，这是一种巨大的浪费，因为虚部总是加起来等于零（我指的是边长为n>1000 的大矩阵，所以我认为这很重要。

现在我只添加矩阵的上三角部分。但我想进一步优化此过程并避免添加复杂部分，因为我不需要它。

我目前使用的是：

std::real(std::accumulate(myMatrix.begin(), myMatrix.end(), std::complex<Real>(0,0)));

这将myMatrix的所有元素加到std::complex<Real>(0,0)上，得到所需的总和。

但这会加上向量的实部和虚部，这是浪费！如何编写最优化的版本，仅添加此矩阵的实部？

---

更新：

虽然我接受了有效的答案，但我发现它比对该矩阵的实部和虚部求和慢。 对于边长为128的矩阵，速度慢了5％-10％。 这很惊人。 如有其他更快的建议，非常感谢。

如果需要其他信息，请提出要求。

- The Quantum Physicist

---

如果您对此帖子进行了“-1”操作，请解释一下。我非常感激。 - The Quantum Physicist

你是否考虑过基于std::thread的并行或者cuda？此外，对于对称矩阵，你只需要访问n(n-1)/2个元素，而不是n*n个元素来进行累加。 - Feng Wang

@FengWang 实际上我正在使用OpenMP进行并行化，而且我已经访问了n(n-1)/2个元素 :) - The Quantum Physicist

通过对大数组进行简单的聚合，性能瓶颈确实在于内存带宽。一个Real有多少字节？如果Real被定义为4字节浮点数，那么64位机器可能会像仅加载其中的前4个字节一样快地加载8个对齐的字节。您可以通过将myMatrix分成由std::vector<float>表示的实部和虚部，然后仅聚合实部来测试这一点，但是只有在不需要频繁在不同表示之间转换时才有益。 - Maarten Hilferink

Answer 1

Real real_sum = std::accumulate(
    myMatrix.cbegin(), myMatrix.cend(), Real{},
    [](Real const acc, std::complex<Real> const& c) { return acc + std::real(c); }
);

- ildjarn

---

顺便说一下，我刚刚发现这个计算速度较慢 :( - The Quantum Physicist

@TheQuantumPhysicist：尝试使用[](Real& acc, std::complex<Real> const c) -> Real& { return acc += c.real(); } - 它是否有任何不同的表现？ - ildjarn

仍然……这比实部和虚部之和慢5%-10%……我真的很惊讶。 - The Quantum Physicist

可能复数加法版本已经进行了SIMD优化。 - Yakk - Adam Nevraumont

Answer 2

std::accumulate 有两个重载版本，其中一个使用操作符：

template< class InputIt, class T, class BinaryOperation >
T accumulate( InputIt first, InputIt last, T init,
              BinaryOperation op );

请提供您自己的内容，而不是默认使用+：

std::accumulate(myMatrix.begin(), myMatrix.end(), Real{0},
    [](Real const& sum, std::complex<Real> const& next){
        return sum + std::real(next);
    });

---

或者，您可以尝试一些有趣的事情，比如使用boost::transform_iterator：

auto make_real = [](std::complex<Real> const& c) {
    return std::real(c);
};

std::accumulate(
    boost::make_transform_iterator(myMatrix.begin(), make_real),
    boost::make_transform_iterator(myMatrix.end(), make_real),
    Real{0});

或者使用 range-v3:

accumulate(myMatrix,
    Real{0},
    std::plus<>{},
    [](std::complex<Real> const& c) { return c.real(); }
);

如果real没有重载，那么以上两个示例将更易于理解。在Boost示例中，您可以提供 std::real<Real>，在第二个示例中，您可以提供&std::complex<Real>::real。

- Barry

Answer 3

使用std::accumulate是绝对必要的吗？

    Real acc = 0;
    for(auto c : myMatrix)
       acc += real(c);

别误会，我支持在适当的情况下使用标准算法，但是这个循环似乎很难在可读性方面被超越。

与我安装的g++-4.8.4附带的实现相比如此。

  template<typename _InputIterator, typename _Tp>
    inline _Tp
    accumulate(_InputIterator __first, _InputIterator __last, _Tp __init)
    {
      // concept requirements
      __glibcxx_function_requires(_InputIteratorConcept<_InputIterator>)
      __glibcxx_requires_valid_range(__first, __last);

      for (; __first != __last; ++__first)
    __init = __init + *__first;
      return __init;
    }

所以你可以看到它们基本上在做同样的事情。

- SirGuy

---

标准算法是针对处理器进行优化的。因此，如果您想要最佳结果，应尽可能使用标准算法。 - The Quantum Physicist

你是不是想要复制对象？如果不是，建议使用auto& c而不是auto c。 - kfsone

@kfsone 我相信编译器会进行避免复制的优化。 - SirGuy