OpenMP 在使用 Armadillo 的代码中生成段错误

Posted 2023-02-17

【中文标题】OpenMP 在使用 Armadillo 的代码中生成段错误

【英文标题】：OpenMP generates segfault in a code that uses Armadillo

【发布时间】：2014-11-18 18:13:10

【问题描述】：

首先，我对 OpenMP 不是很熟悉。我想使用 OpenMP 减少我的 C++ 代码的执行时间，这涉及到几百次矩阵对角化迭代。我不试图平行每个对角化（根据 Armadillo 的文档，可以通过强制 Armadillo 使用 OpenBLAS 库来实现）；相反，我想在 8 核机器上的线程之间分配负载。

当我得到“分段错误”时，访问内存似乎有问题。我想知道是我做的不对还是问题出在犰狳创建和操作矩阵的方式上。

这里有一个最小的代码，它抓住了我遇到的问题的本质。这个想法是将 1000 个 100x100 矩阵对角化，并将它们的特征值存储在一个文件中。

#include<iostream>
#include<armadillo>
#include <fstream>
#include<omp.h>         

int main()
   
    std::ofstream  File;
    File.open("./RESULTS.dat");

    arma::mat M;   //THE MATRIX TO BE DIAGONALIZED
    arma::mat Eigenvecs;  //EIGENVECTORS
    arma::vec Eigenval;   //EIGENVALUES

    arma::mat RESULTS; //STORING EIGENVALUES TEMPORARILY 

    //DISTRIBUTING THE ITTERATIONS AMONG CORES USING OpenMP
    #pragma omp parallel shared(RESULTS) private(M,Eigenvecs,Eigenval)
    
    #pragma omp parallel for ordered schedule(guided)
    for( int i = 0 ; i < 1000; i++ )
    
        M = arma::randu<arma::mat>(200,200); //CREATING A RANDOM MATRIX
        M = 0.5*(M + M.t() );  //TO GUARANTEE THAT THE M IS NORMAL

        arma::eig_sym( Eigenval , Eigenvecs , M ); //DIAGONALIZING "M"

        RESULTS = arma::join_vert(RESULTS,Eigenval.t()); //CONCATENATING EIGENVALUES TO THE MATRIX "RESULTS" 
    
    

    File << RESULTS;  //WRITING "RESULTS" TO THE FILE
    File.close();

    return 0;

当我运行此代码时，似乎负载分配正确（我使用 htop 监控机器上的内核），但最后我得到“分段错误”。

【问题讨论】：

如果在调试器下运行会发生什么？段错误发生在哪里？

所以，我很确定你是实际上在使用某种 BLAS+LAPACK 库：without them, Armadillo cannot perform diagonalization。那么问题就变成了：底层的 BLAS+LAPACK 库是线程安全的吗？这在很大程度上取决于您使用的是哪种 BLAS 实现（例如 OpenBLAS？MKL？）以及它的编译方式。

我明白了。感谢您的评论。

【参考方案1】：

我想我明白了，虽然我不知道为什么这能解决问题。这是修改后的代码：

#include<iostream>
#include<ctime>
#include<armadillo>
#include <fstream>
#include<omp.h> 


int main()
   
    std::ofstream  File;
    File.open("./RESULTS.dat");

    arma::mat M;
    arma::mat RESULTS;
    arma::mat Eigenvecs;
    arma::vec Eigenval;

    int i,tid=0;

    //DISTRIBUTING THE ITTERATIONS AMONG CORES USING OpenMP
    #pragma omp parallel for ordered schedule(guided) shared(RESULTS) private(M,Eigenvecs,Eigenval,tid)
    for(i = 0 ; i < 10; i++ )
    
        tid = omp_get_thread_num();

        M = arma::randu<arma::mat>(400,400);
        M = 0.5*(M + M.t() );
        //std::cout << M << std::endl;

        arma::eig_sym( Eigenval , Eigenvecs , M );

        #pragma omp ordered
        RESULTS = arma::join_vert(RESULTS,Eigenval.t());

        std::cout << i << " , " << tid  << " --- " << Eigenval(0) << std::endl;
    

    File << RESULTS;   
    File.close();

    return 0;

【讨论】：

与其在#pragma omp 中将变量标记为私有，不如在循环中显式声明它们可能更安全。这样你就可以完全确保每个线程都有自己的变量（以及内存），而不是依赖于 OpenMP（它可能有一些关于非平凡矩阵对象私有的巴洛克规则）