犰狳的 eig_sym 中的线程安全问题

Question

我对通过 eig_sym 在犰狳中的特征分解有疑问。当我尝试并行计算多组特征值和特征向量时，特征向量不时是

• 不正交
• 未标准化
• 甚至没有相关矩阵的特征向量。

如果每次只运行一个计算，这个问题就会消失（所以这似乎是一些线程安全问题）。一旦两个计算并行运行，问题就会再次出现。奇怪的是，特征值似乎在每种情况下都是正确的。

//compile with: g++ -std=c++11 -pthread -larmadillo -o evecs armadillo_evecs.cpp
#include <iostream>
#include <armadillo>
#include <assert.h>
#include <future>
#include <vector>
using namespace std;

void testcalc() {
    // set up random symmetric matrix
    arma::mat r = arma::randu<arma::mat> (100, 100);
    r  = r.t() * r;

    arma::vec eval;
    arma::mat evec;
    // calculate eigenvalues and -vectors
    assert(arma::eig_sym(eval, evec, r));
    arma::mat test = evec.t() * evec;

    // Check whether eigenvectors are orthogonal, (i. e. matrix 'test' is diagonal)
    assert(arma::norm(test - arma::diagmat(test)) < 1.0e-10);
}


int main() {
    // start 100 eigenvalue (+vector) calculations
    vector<future<void>> fus;
    for (size_t i = 0; i < 100; i++) {
        // try parallel evaluation ... fails sometimes 
        fus.push_back(async(launch::async, &testcalc));

        // try sequential evaluation ... works fine
        // future<void> f = async(launch::async, &testcalc);
        // f.get(); // Wait until calculation has finished, before starting new one
    }

    // wait until calculations have finished
    for(auto it = fus.begin(); it != fus.end(); it++) {
        it->get();
    }

    return 0;
}

所以在断言上面的代码中

assert(arma::norm(test - arma::diagmat(test)) < 1.0e-10);

有时会失败。这可能是底层库的问题（我读到 lapack 有一些线程安全问题）？我真的不知道，从哪里开始寻找。

Answer 1

与“滚动您自己的并行化”不同，使用底层库已经提供的并行化更容易、更安全。

因此，不要使用参考 BLAS 和 LAPACK，而是使用多线程版本，例如 OpenBLAS 或 Intel MKL。有关详细信息，请参阅犰狳FAQ。