CODEWITHSUNDEEP
c++c++11matrixvectorcomplex-numbers
The Quantum Physicist
The Quantum Physicist

Reputation: 26256

std::accumulate() only the real part of a complex std::vector

I used to take the sum of a symmetric (Hermitian) matrix (the matrix is in an std::vector), which is a huge waste because the imaginary part will always add to zero (I'm talking about huge matrices with side-length of n>1000, so I think it makes a difference.

So now I only add the upper-triangular part of the matrix. But I want to optimize this further and avoid adding the complex part because I don't need it.

What I currently use is:

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

This adds all the element of myMatrix to std::complex<Real>(0,0), resulting in the sum I need.

But this will add real and imaginary components of my vector, which is a waste! How can I write the most optimized version of this that adds only the real part of this matrix?

UPDATE:

While I accepted the answer that works, I discovered that it's slower than summing the real and imaginary parts of my matrix. It's slower by 5%-10% for a matrix with side-length of 128. This is surprising. Any other faster suggestions are highly appreciated.

Please ask if you require additional information.

Upvotes: 4

Views: 1559

Answers (3)

ildjarn
ildjarn

Reputation: 62975

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

Upvotes: 8

SirGuy
SirGuy

Reputation: 10770

Is using std::accumulate absolutely necessary?

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

Don't get me wrong, I'm in favour of using standard algorithms when they suit, but it seems hard to beat this loop in readability.

This is as compared to the implementation that came with my g++-4.8.4 install:

  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;
    }

So you can see they're pretty much doing the same thing.

Upvotes: 1

Barry
Barry

Reputation: 302852

std::accumulate has two overloads, one of which takes an operator:

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

Just provide your own instead of defaulting to +:

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

Alternatively, you could do something fun like use 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});

Or with range-v3:

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

Both of the above would be much nicer if real didn't have overloads and you could just provide std::real<Real> in the boost example and &std::complex<Real>::real in the second.

Upvotes: 4

Related Questions