Как эффективно и точно распараллелить Eigen Pardiso Solve - PullRequest
Новая
       50

Как эффективно и точно распараллелить Eigen Pardiso Solve

0 голосов
/ 27 марта 2020

Чтобы ускорить код физического моделирования, я использую Pardiso с Eigen. По сравнению с альтернативными решателями, такими как Eigen SparseLU, Pardiso дает мне значительное увеличение шага факторизации, которое увеличивается с увеличением количества потоков.

В моем коде моделирования, однако, RHS этого решения не один вектор, но матрица. Я бы предположил, что это будет хорошо с параллельной точки зрения, потому что каждое обратное решение является независимым, но я не замечаю никакого ускорения на этапе решения при увеличении числа потоков.

Из-за этого я попытался использовать OpenMP для параллельного решения. Хотя, по-видимому, быстрее, он продолжает возвращать неверный результат. Я настроил тест примерно так ... 

SparseMatrix<double> A = ...;
SparseMatrix<double> B = ...;
mkl_set_num_threads(t);

// Factorize
pardiso.factorize(A);

// Solve1
SparseMatrix<double> X1 = pardiso.solve(B);

// Solve2
SparseMatrix<double> X2(A.rows(), B.cols());
#pragma omp parallel num_threads(t)
{
#pragma omp for
    {
         for (int n = 0; n < B.cols(); ++n)
        {
             X2.col(n) = pardiso.solve(B.col(n));
        }
    }
}

X2.isApprox(X1); // I know X1 is giving me the correct matrix, but this return false

Мне любопытно, если кто-нибудь знает, почему Solve1 имеет одинаковую скорость независимо от num_threads, или если кто-нибудь знает, почему Solve2 возвращает неправильную матрицу .

Я на Windows использую VS2019, компилирую с помощью компилятора Intel C ++, установил Use MKL на параллельный режим и включил OpenMP. Я проверял, что OpenMP работает и работает, Solve2 работает правильно, когда omp for использует один поток, а Solve2 работает правильно, используя OpenMP и Eigen SparseLU. Я также просмотрел все эти предыдущие вопросы, чтобы они не преобладали.

Eigen :: PardisoLU (MKL) в параллельном сценарии

Вызов многопоточного MKL из параллельная область openmp

Количество потоков функций Intel MKL внутри параллельных областей OMP

Я все еще получаю неправильный вывод. Установка mkl_set_dynamic(1) и omp_set_nested(1) до того, как все ничего не изменит.

...