Различные результаты при использовании более одного потока - PullRequest
Новая
       14

Различные результаты при использовании более одного потока

0 голосов
/ 21 ноября 2018

Следующий код выполняет транспонирование умножения матрицы на вектор, где матрица является разреженной и сохраняется в формате CSR.В зависимости от количества потоков результат будет разным.Я думаю, что причина заключается в одновременном доступе к памяти и дополнении.Есть ли способ использовать многопоточность, но оставить результат таким же, как и для однопоточности? 

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

Ответы [ 2 ]

0 голосов
/ 30 ноября 2018

Без воспроизводимого примера, я должен был предположить в контексте этого кода: 

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

С некоторыми подходящими переменными мы можем упростить код, чтобы мы видели, что происходит: 

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

Теперь мы можем видеть в теле цикла, который мы присваиваем записи результата, в которую другие потоки также могут записывать.Нам нужно сериализовать доступ к result[cols[j]], чтобы только один поток одновременно выполнял +=.Мы можем сделать это, пометив эту операцию как неделимую, используя ключевое слово OMP atomic: 

#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;
        }
    }
0 голосов
/ 22 ноября 2018

Операция приращения действительно должна быть атомарной, чтобы убедиться, что два потока не обновляют значение одновременно, что будет условием гонки .Это сделает код правильным, но поскольку атомарный инкремент медленнее, чем стандартный инкремент, он может сделать параллельный код медленнее, чем последовательный код в целом.
На производительность также может влиять false-share : еслиразмер вектора не намного больше числа потоков, часто бывает так, что два потока пытаются увеличить значения, принадлежащие одной и той же строке кэша, и много времени уходит на синхронизацию кэша между процессорами. 

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

При этом для данного элемента в векторе результатов порядок добавления различных продуктов для создания этого элемента не будет одинаковым в последовательном и параллельном случаях.Ошибки округления будут складываться по-разному, и не следует ожидать, что последовательные и параллельные результаты будут одинаковыми, а разница между этими результатами будет небольшой по сравнению с точностью float или double формат.

...