c - 慢 OMP 与串行

Question

我正在尝试优化从 R 调用的 C 子例程，该子例程占用了我试图解决的问题的约 60% 的计算时间。这比纯粹用 R 编码时的 86% 有所下降。我的 C 代码中的绝大多数执行时间发生在嵌套的 for 循环中，因此这似乎是尝试使用 OpenMP 进行并行化的明显候选者。我尝试过这样做，但结果不同 - 最好的情况是所用时间比不使用 OMP 稍差一些，最坏的情况是性能与线程数量成反比。最快版本的代码如下:

#include <R.h>
#include <Rmath.h>
#include <omp.h>
void gradNegLogLik_c(double *param, double *delta, double *X, double *M, int *nBeta, int *nEpsilon, int *nObs, double *gradient){
  // ========================================================================================
  // param:     double[nBeta + nEpsilon]        values of parameters at which to evaluate gradient
  // delta:     double[nObs]                    satellite - buoy differences
  // X:         double[nObs * (nBeta + nEpsilon)]  design matrix for mean components (i.e. beta terms)
  // M:         double[nObs * (nBeta + nEpsilon)]  design matrix for variance components (i.e. epsilon terms)
  // nBeta:     int                             number of mean terms
  // nEpsilon:  int                             number of variance terms
  // nObs:      int                             number of observations
  // gradient:  double[nBeta + nEpsilon]        output array of gradients
  // ========================================================================================

  // ========================================================================================
  // local variables
  size_t i, j, ind;
  size_t nterms = *nBeta + *nEpsilon;  
  size_t nbeta  = *nBeta;
  size_t nepsilon  = *nEpsilon;
  size_t nobs = *nObs; 
  // allocate local memory and set to zero
  double *sigma2        = calloc( nobs     , sizeof(double) );
  double *fittedValues  = calloc( nobs     , sizeof(double) );
  double *residuals     = calloc( nobs     , sizeof(double) );
  double *beta          = calloc( nbeta    , sizeof(double) );
  double *epsilon2      = calloc( nepsilon , sizeof(double) );  
  double *residuals2    = calloc( nobs     , sizeof(double) );
  double gradBeta, gradEpsilon;

  // extract beta and epsilon terms from param
  // =========================================
  for(i = 0 ; i < nbeta ; i++){
    beta[i]    = param[ i ];
    epsilon2[i] = param[ nbeta + i ];
  }

  // Initialise gradient to zero for return value
  // =========================================  
  for( i = 0 ; i < nterms ; i++){
    gradient[i] = 0;
  }

  // calculate sigma, fitted values and residuals
  // ============================================  
  for( i = 0 ; i < nbeta ; i++){
    for( j = 0 ; j < nobs ; j++){
      ind = i * nobs + j;
      sigma2[j] += M[ind] * epsilon2[i];
      fittedValues[j] += X[ind] * beta[i];
    }
  }

  for( j = 0 ; j < nobs ; j++){
    // calculate reciprocal as this is what we actually use and 
    // we only want to do it once.        
    sigma2[j] = 1 / sigma2[j]; 
    residuals[j] = delta[j] - fittedValues[j];    
    residuals2[j] = residuals[j]*residuals[j];
  } 

  // Loop over all observations and calculate value of (negative) derivative  
  // =======================================================================    
#pragma omp parallel for private(i, j, ind, gradBeta, gradEpsilon)\
  shared(gradient, nbeta, nobs, X, M, sigma2, fittedValues, delta, residuals2) \
  default(none) 
  for( i = 0 ; i < nbeta ; i++){
    gradBeta = 0.0;
    gradEpsilon = 0.0;
    for(j = 0 ; j < nobs ; j++){
      ind = i * nobs + j;
      gradBeta    -= -1.0*X[ind] * sigma2[j]*(fittedValues[j] - delta[j]);
      gradEpsilon -= 0.5*M[ind] * sigma2[j]*(residuals2[j] * sigma2[j] - 1);
    }
    gradient[i] = gradBeta;
    gradient[nbeta + i] = gradEpsilon;
  }

  // End of function
  // free local memory
  free(sigma2);
  free(fittedValues);
  free(residuals);
  free(beta);
  free(epsilon2);
  free(residuals2);
}

nObs 是订单 10000。

nBeta 的范围是 20 到几百。

nEpsilon = nBeta，当前未使用。

在搜索了这个网站、一个下午的谷歌搜索并尝试了不同的事情之后，我似乎无法做出任何进一步的改进。我的第一个想法是错误共享 – 我尝试了各种方法，例如展开外循环以一次设置 8 个渐变[] 元素，以创建一个临时填充数组来存储结果。我还尝试了不同的组合共享、私有(private)和第一私有(private)。这些似乎都没有改善事情，而且我的最快执行时间并行比串行稍微差一些。在我花更多时间讨论这个问题之前，这引出了两个问题:

• 我的问题(重复约 9000 次同一组计算 20 - 900 次)是否太小，不值得使用 OMP？
• 我是否遗漏或做错了什么？

我怀疑是后者，因为我在使用 C 和 OMP 时相对缺乏经验。任何帮助/想法将不胜感激。

(仅供引用，我在具有 16 核和 192GB 内存的 SLED11 服务器上运行，并使用 GCC 4.7.2 编译我的 C 代码)。其他用户正在使用该服务器，但 OMP 与串行代码的相对性能似乎与其他用户无关。

提前致谢，

戴夫。

编辑:有关信息，我使用的编译命令是

gcc -I/RHOME/R/3.0.1/lib64/R/include -DNDEBUG -I/usr/local/include -fpic \
    -std=c99 -Wall -pedantic –O3 -fopenmp -c src/gradNegLogLik_call.c \
     -o src/gradNegLogLik_call.o

大多数标志是由 R CMD SHLIB 命令设置的 - 我已手动添加了 -O3 -fopenmp。

Answer 1

在回答我为加速代码所做的事情之前，为我上面的问题提供一些背景可能会很有用(尽管这是在不使用 OMP 的情况下实现的)。

我最初编写的 C 函数是为了计算与 R optim() 命令和 L-BFGS-B 方法一起使用的对数似然函数的梯度。对于每次调用 optim，我的对数似然函数和梯度函数都会被调用约 100 次，因为 optim 会找到最佳解决方案。结果，正如 Rprof 所预期和报告的那样，这两个函数占用了我的大部分执行时间，转换为 C 以提高代码效率的两个目标也是如此。

将我的两个函数转换为 C 语言并优化该代码，使我的 optim 调用时间从每次调用 1.88 秒的平均运行时间减少到每次调用 0.25 秒。这将我的处理时间从大约 1 个月减少到几天。影响最大的更改(除了调用 C 之外)是更改嵌套循环的顺序。选择原始顺序是由于 R 存储矩阵的方式，并且选择它是为了避免每次调用 C 函数时都必须转置矩阵。认识到每次调用 optim() 时只需执行一次转置，而不是像我最初编码的那样每次 C 调用都需要执行一次转置，与更改 C 函数中的顺序的影响/好处相比，这是一个很小的开销。

考虑到速度的提高，我们必须证明在这方面花费更多时间是合理的。下面给出了我的梯度函数的最终版本(根据我的原始帖子)。

请注意，虽然我在 R 中从使用 .C 更改为 .Call(因此对函数参数等进行了更改)，但这本身并不能说明速度的提高。

#include <R.h>
#include <Rmath.h>
#include <Rinternals.h>
#include <omp.h>
SEXP gradNegLogLik_call(SEXP param ,SEXP delta, SEXP X, SEXP M, SEXP nBeta, SEXP nEpsilon){

  // local variables
  double *par, *d;
  double *sigma2, *fittedValues, *residuals, *grad, *Xuse, *Muse;
  double val, sig2, gradBeta, gradEpsilon;
  int n, m, ind, nterms, i, j;

  SEXP gradient;

  // get / associate parameters with local pointer
  par = REAL(param);
  Xuse = REAL(X);
  Muse = REAL(M);
  d = REAL(delta);
  n = LENGTH(delta);
  m = INTEGER(nBeta)[0];
  nterms = m + m;


  // allocate memory
  PROTECT( gradient = allocVector(REALSXP, nterms ));
  // set pointer to real portion of gradient
  grad = REAL(gradient);

  // set all gradient terms to zero
  for(i = 0 ; i < nterms ; i++){
      grad[i] = 0.0;
  }


  sigma2 = Calloc(n, double );
  fittedValues = Calloc(n, double );
  residuals = Calloc(n, double );

  // calculate sigma, fitted values and residuals
  for(i = 0 ; i < n ; i++){
    val = 0.0;
    sig2 = 0.0;
    for(j = 0 ; j < m ; j++){
      ind = i*m + j;
      val +=  Xuse[ind]*par[j];
      sig2 += Muse[ind]*par[j+m];
    }
    // calculate reciprocal of sigma as this is what we actually use
    // and we only want to do it once
    sigma2[i] = 1.0 / sig2;
    fittedValues[i] = val;
    residuals[i] = d[i] - val;
  }


  // now loop over each paramter and calculate derivative
  for(i = 0 ; i < n ; i++){
    gradBeta = -1.0*sigma2[i]*(fittedValues[i] - d[i]);
    gradEpsilon = 0.5*sigma2[i]*(residuals[i]*residuals[i]*sigma2[i] - 1);     
    for(j = 0 ; j < m ; j++){
        ind = i*m + j;
        grad[j] -= Xuse[ind]*gradBeta;
        grad[j+m] -= Muse[ind]*gradEpsilon;
    }            
  }

  UNPROTECT(1);
  Free(sigma2);
  Free(residuals);
  Free(fittedValues);
  // return array of gradients  
  return gradient;
}