Ytukyg

The following code performs a transpose matrix-vector multiplication where the matrix is sparse and stored in CSR format. Depending on the number of threads the result is different. I think the reason is the concurrent memory access and addition.
Is there a way to use multithreading but keep the result the same as for single threading?

#pragma omp parallel for num_threads(m_numthreads)

    for (int i = 0; i < matrix.get_rowptr()->size() - 1; ++i)

    {

        for (int j = matrix.get_rowptr()->operator(i); j < matrix.get_rowptr()->operator(i + 1); ++j)

        {

            result[matrix.get_columnindex()->operator(j)] += matrix.get_value()->operator(j) * vector1[i];

        }

    }

The increment operation should indeed be made atomic to make sure two threads are not updating the value at the same time, which would be a race condition. This will make the code correct, but since atomic increment is slower than standard increment, it may make the parallel code slower than the serial one overall.

Performance may also be affected by false-sharing: if the vector size is not much greater than the number of threads, it will often happen that two threads attempt to increment values belonging to the same cache line and much time is spent synchronizing the cache between CPUs.

#pragma omp parallel for num_threads(m_numthreads)

    for (int i = 0; i < matrix.get_rowptr()->size() - 1; ++i)

        for (int j = matrix.get_rowptr()->operator(i); j < matrix.get_rowptr()->operator(i + 1); ++j)

        {

            #pragma omp atomic

            result[matrix.get_columnindex()->operator(j)] += matrix.get_value()->operator(j) * vector1[i];

        }

This being said, for a given element in the result vector, the order in which the various products are added to make that element will not be the same in the serial and parallel cases. Rounding errors will add differently and one should not expect the serial and parallel results to be the same, but rather the difference between these results to be small relative to the precision of the float or double format.

Without a reproducible example, I've had to guess at the context of this code:

#include <vector>



class SparseMatrix

{

    std::vector<double> values = { 5, 8, 3, 6, };

    std::vector<std::size_t> rows = { 0, 0, 2, 3, 4, };

    std::vector<std::size_t> cols = { 1, 2, 1, 1, };



public:

    const std::vector<std::size_t> *get_rowptr() const { return &rows; };

    const std::vector<std::size_t> *get_columnindex() const { return &cols; };

    const std::vector<double> *get_value() const { return &values; };

};



#include <array>

#include <iostream>

int main()

{

    SparseMatrix matrix;



    std::array<double, 4> result{};

    std::array<double, 4> vector1{ 1, 2, 3, 4 };



#pragma omp parallel for

    for (int i = 0; i < matrix.get_rowptr()->size() - 1; ++i)

    {

        for (int j = matrix.get_rowptr()->operator(i); j < matrix.get_rowptr()->operator(i + 1); ++j)

        {

            result[matrix.get_columnindex()->operator(j)] += matrix.get_value()->operator(j) * vector1[i];

        }

    }



    for (auto const& i: result)

        std::cout << i << " ";

    std::cout << 'n';

}

With some suitable variables we can simplify the code so we see what's going on:

    auto const& rows = *matrix.get_rowptr();

    auto const& cols = *matrix.get_columnindex();

    auto const& values = *matrix.get_value();



 #pragma omp parallel for

    for (std::size_t i = 0;  i < rows.size() - 1;  ++i)

    {

        for (std::size_t j = rows[i];  j < rows[i+1];  ++j)

        {

            result[cols[j]] += values[j] * vector1[i];

        }

    }

Now we can see in the body of the loop that we're assigning to a result entry that other threads may also be writing to. We need to serialize the access to result[cols[j]] so that only one thread at a time is executing the +=. We can do that by marking that operation indivisible, using the OMP atomic keyword:

#pragma omp parallel for

    for (std::size_t i = 0;  i < rows.size() - 1;  ++i)

    {

        for (std::size_t j = rows[i];  j < rows[i+1];  ++j)

        {

            auto& dest = result[cols[j]];

            auto const term = values[j] * vector1[i];

#pragma omp atomic

            dest += term;

        }

    }