此OpenMP代码崩溃linux

Question

我正在处理一些代码以减少使用openmp的矩阵。我有两个版本，都让我的Ubuntu和Fedora安装崩溃。通过硬我的意思是我的鼠标和键盘没有响应，即使我点击PC塔上的重置按钮，它也不会重新启动。我必须按住电源按钮。奇怪的是，代码在运行几分钟后崩溃了。它并没有消耗大量的内存（我认为750 MB是很小的，因为我有16gb的内存）。

#include <iostream> 
#include <cstddef> 
#include <cstring> 
#include <iomanip> 
#include <cstdlib> 
#include <ctime> 
#include <cmath> 

using namespace std; 

class Matrix 
{ 
    public: 

    Matrix(size_t rows, size_t cols): 
     data(0), w(rows), h(cols) 
    { 
     data = new double[w * h]; 
     memset(data, 0, sizeof(double) * w * h); 
    } 

    ~Matrix() 
    { 
     if(data) 
     { 
      delete[] data; 
      w = h = 0; 
      data = 0; 
     } 
    } 

    double* operator[](size_t row) 
    { 
     return data + row * w; 
    } 

    const double* operator[](size_t row) const 
    { 
     return data + row * w; 
    } 

    size_t width() const 
    { 
     return w; 
    } 

    size_t height() const 
    { 
     return h; 
    } 

    void scale_row(size_t row, double x) 
    { 
     double* prow = (*this)[row]; 

     for(size_t i = 0; i < w; i++) 
      prow[i] *= x; 
    } 

    void add_row(size_t dest_row, size_t source_row, double scaling = 1.0) 
    { 
     if(dest_row == source_row) 
     { 
      scale_row(dest_row, 1.0 + scaling); 
      return; 
     } 

     double* __restrict__ drow = (*this)[dest_row]; 
     double* __restrict__ srow = (*this)[source_row]; 

     for(size_t i = 0; i < w; i++) 
      drow[i] += srow[i] * scaling; 
    } 

    void swap_rows(size_t r1, size_t r2) 
    { 
     if(r1 == r2) 
      return; 

     double* __restrict__ a = (*this)[r1]; 
     double* __restrict__ b = (*this)[r2]; 

     #pragma omp parallel for simd 
     for(size_t i = 0; i < w; i++) 
     { 
      double tmp = a[i]; 
      a[i] = b[i]; 
      b[i] = tmp; 
     } 
    } 

    double* find_leading(size_t row) 
    { 
     double* ptr = (*this)[row]; 
     for(size_t i = 0; i < w; i++) 
      if(ptr[i]) 
       return ptr + i; 
     return 0; 
    } 

    void clamp_zeros(double threshold = 1e-12) 
    { 
     #pragma omp parallel for simd 
     for(size_t i = 0; i < w * h; i++) 
     { 
      if(fabs(data[i]) < threshold) 
       data[i] = 0; 
     } 
    } 

    void row_reduce(Matrix* mirror = 0) 
    { 
     for(size_t r1 = 0; r1 < h; r1++) 
     { 
      double* lead = find_leading(r1); 
      if(!lead) 
       continue; 

      size_t rank = lead - (*this)[r1]; 
      if(mirror) 
       mirror->scale_row(r1, 1.0/*lead); 
      scale_row(r1, 1.0/*lead); 

      #pragma omp parallel for 
      for(size_t r2 = 0; r2 < h; r2++) 
      { 
       if(r2 == r1 || (*this)[r2][rank] == 0) 
        continue; 
       if(mirror) 
        mirror->add_row(r2, r1, -(*this)[r2][rank]); 
       add_row(r2, r1, -(*this)[r2][rank]); 
      } 
      clamp_zeros(); 
     } 

     size_t zero_count = 0; 
     for(size_t r = 0; r < h; r++) 
     { 
      double* lead = find_leading(r); 
      if(lead) 
      { 
       size_t rank = lead - (*this)[r]; 
       swap_rows(rank, r); 
       if(mirror) 
        mirror->swap_rows(rank, r); 
      } 
      else 
      { 
       size_t with = h - ++zero_count; 
       swap_rows(r, with); 
       if(mirror) 
        mirror->swap_rows(r, with); 
      } 
     } 
    } 

    private: 

    double* data; 
    size_t w, h; 
}; 

ostream& operator<<(ostream& o, const Matrix& m) 
{ 
    o << setprecision(2); 
    for(size_t j = 0; j < m.width(); j++) 
    { 
     o << "----------"; 
    } 
    o << "--\n"; 
    for(size_t i = 0; i < m.height(); i++) 
    { 
     o << "|"; 
     for(size_t j = 0; j < m.width(); j++) 
     { 
      o << setw(10) << m[i][j]; 
     } 
     o << "|\n"; 
    } 
    for(size_t j = 0; j < m.width(); j++) 
    { 
     o << "----------"; 
    } 
    o << "--"; 
    return o; 
} 

int main() 
{ 
    srand(time(0)); 
    Matrix m (10000, 10000); 

    for(int i = 0; i < m.height(); i++) 
    { 
     for(int j = 0; j < m.width(); j++) 
     { 
      m[i][j] = rand() % 100; 
     } 
    } 

    time_t start = time(0); 
    m.row_reduce(); 
    time_t end = time(0); 
    cout << m[0][2] << endl; 
    cout << "dt = " << (end - start) << endl; 
    return 0; 
} 

我也尝试了另一种愚蠢的简单OMP程序，看它是否会崩溃我的系统，这一次没有。

double sum = 0.0; 

double start = omp_get_wtime(); 
#pragma omp parallel for reduction(+:sum) 
for(long long i = 1; i < 100000000000000LL; i++) 
{ 
    sum += 1.0/((double)i * i); 
} 
printf("%lf %lf\n", omp_get_wtime() - start, sum); 

我尝试的第一个跑进了同样的问题，当我跑了它在Ubuntu 15.04使用gcc 4.9和Fedora 22用gcc编译5.1编译。

当我运行它没有openmp它工作正常。另外，如果我尝试使用像2000x2000矩阵那样的较小数据，它可以正常工作（当我尝试使用10,000x10,000矩阵时发生崩溃）。

似乎工作正常我的笔记本电脑，这也是运行的Ubuntu 15.04。

Answer 1

为了支持与OpenMP 2.0的兼容性，我做了一些代码修改，我可以告诉你，你的代码运行良好（Windows 7，Visual Studio 2008）。内存消耗约800MB。

输出：

DT = 2881

这是你改变代码。

//////////////////////////////////////////////////////////////// 
// OpenMP test function 
#include <iostream> 
#include <cstddef> 
#include <cstring> 
#include <iomanip> 
#include <cstdlib> 
#include <ctime> 
#include <cmath> 
#include <omp.h> 

using namespace std; 

class Matrix 
{ 
    public: 

    Matrix(size_t rows, size_t cols): 
     data(0), w(rows), h(cols) 
    { 
     data = new double[w * h]; 
     memset(data, 0, sizeof(double) * w * h); 
    } 

    ~Matrix() 
    { 
     if(data) 
     { 
      delete[] data; 
      w = h = 0; 
      data = 0; 
     } 
    } 

    double* operator[](size_t row) 
    { 
     return data + row * w; 
    } 

    const double* operator[](size_t row) const 
    { 
     return data + row * w; 
    } 

    size_t width() const 
    { 
     return w; 
    } 

    size_t height() const 
    { 
     return h; 
    } 

    void scale_row(size_t row, double x) 
    { 
     double* prow = (*this)[row]; 

     for(size_t i = 0; i < w; i++) 
      prow[i] *= x; 
    } 

    void add_row(size_t dest_row, size_t source_row, double scaling = 1.0) 
    { 
     if(dest_row == source_row) 
     { 
      scale_row(dest_row, 1.0 + scaling); 
      return; 
     } 

     double* drow = (*this)[dest_row]; 
     double* srow = (*this)[source_row]; 

     for(size_t i = 0; i < w; i++) 
      drow[i] += srow[i] * scaling; 
    } 

    void swap_rows(size_t r1, size_t r2) 
    { 
     if(r1 == r2) 
      return; 

     double* a = (*this)[r1]; 
     double* b = (*this)[r2]; 

     #pragma omp parallel for schedule(dynamic) 
     for(int i = 0; i < w; i++) 
     { 
      double tmp = a[i]; 
      a[i] = b[i]; 
      b[i] = tmp; 
     } 
    } 

    double* find_leading(size_t row) 
    { 
     double* ptr = (*this)[row]; 
     for(int i = 0; i < w; i++) 
      if(ptr[i]) 
       return ptr + i; 
     return 0; 
    } 

    void clamp_zeros(double threshold = 1e-12) 
    { 
     #pragma omp parallel for schedule(dynamic) 
     for(int i = 0; i < w * h; i++) 
     { 
      if(fabs(data[i]) < threshold) 
       data[i] = 0; 
     } 
    } 

    void row_reduce(Matrix* mirror = 0) 
    { 
     for(size_t r1 = 0; r1 < h; r1++) 
     { 
      double* lead = find_leading(r1); 
      if(!lead) 
       continue; 

      size_t rank = lead - (*this)[r1]; 
      if(mirror) 
       mirror->scale_row(r1, 1.0/*lead); 
      scale_row(r1, 1.0/*lead); 

      #pragma omp parallel for schedule(dynamic) 
      for(int r2 = 0; r2 < h; r2++) 
      { 
       if(r2 == r1 || (*this)[r2][rank] == 0) 
        continue; 
       if(mirror) 
        mirror->add_row(r2, r1, -(*this)[r2][rank]); 
       add_row(r2, r1, -(*this)[r2][rank]); 
      } 
      clamp_zeros(); 
     } 

     size_t zero_count = 0; 
     for(size_t r = 0; r < h; r++) 
     { 
      double* lead = find_leading(r); 
      if(lead) 
      { 
       size_t rank = lead - (*this)[r]; 
       swap_rows(rank, r); 
       if(mirror) 
        mirror->swap_rows(rank, r); 
      } 
      else 
      { 
       size_t with = h - ++zero_count; 
       swap_rows(r, with); 
       if(mirror) 
        mirror->swap_rows(r, with); 
      } 
     } 
    } 

    private: 

    double* data; 
    size_t w, h; 
}; 

ostream& operator<<(ostream& o, const Matrix& m) 
{ 
    o << setprecision(2); 
    for(size_t j = 0; j < m.width(); j++) 
    { 
     o << "----------"; 
    } 
    o << "--\n"; 
    for(size_t i = 0; i < m.height(); i++) 
    { 
     o << "|"; 
     for(size_t j = 0; j < m.width(); j++) 
     { 
      o << setw(10) << m[i][j]; 
     } 
     o << "|\n"; 
    } 
    for(size_t j = 0; j < m.width(); j++) 
    { 
     o << "----------"; 
    } 
    o << "--"; 
    return o; 
} 

int main() 
{ 
    int iMaxThreads = omp_get_max_threads(); 
    omp_set_num_threads(iMaxThreads); 

    omp_set_dynamic(false); 
    omp_set_nested(true); 

    srand(time(0)); 
    Matrix m (10000, 10000); 

    for(int i = 0; i < m.height(); i++) 
    { 
     for(int j = 0; j < m.width(); j++) 
     { 
      m[i][j] = rand() % 100; 
     } 
    } 

    time_t start = time(0); 
    m.row_reduce(); 
    time_t end = time(0); 
    cout << m[0][2] << endl; 
    cout << "dt = " << (end - start) << endl; 
    return 0; 
} 

Answer 2

我使用GCC 4.9.2在Linux 3.19.0-26-generic＃28-Ubuntu 64位上测试了你的代码。

如下图所示，它使用了一些RAM，但是我还没有崩溃，CPU时间为11分钟，RAM保持稳定。

Answer 3

事实证明，这是不是一个真正的编程问题都没有。我的程序显然在我的笔记本电脑和其他人的系统上运行良好。我刚刚运行了25位pi基准测试数据。这会导致我的电脑以完全相同的方式崩溃。在此之后，我尝试了y cruncher的windows版本。它在大约30秒后引起蓝屏。我想这是一个硬件问题，发生在内存或cpu被推到一段时间后真的很难。现在我有我的借口升级到skylake cpu。

更新： 我设法解决它。前一阵子，我在华硕主板上的EZ XMP开关上轻轻一按。这是为了自动超过时钟记忆。之前我曾尝试过在主板上设置时钟的CPU，并且它们总是使我的系统不稳定。记忆似乎工作，所以我离开它，忘了它。我想这不是这种情况，它导致我的崩溃。现在我把它关掉了，我的两个cruncher和openmp都可以运行完成。