c++ - 将 openmp 与 odeint 和自适应步长一起使用

Question

我正在尝试使用 openmp 来并行化我的代码。当我使用恒定步长时一切正常，但是当我使用自适应步进器运行相同的代码时，我会收到我不理解的错误。

以下是代码的基本部分:

using namespace std;
using namespace boost::numeric::odeint;
const int jmax = 10;
typedef  double value_type;
typedef boost::array<value_type ,2*(jmax+1) > state_type;


//The step function
void rhs( const state_type A , state_type &dAdt , const  value_type t )

{

 value_type RHStemp0;
 value_type RHStemp1;
//We will write the RHS of the equations a  big sum
#pragma omp parallel for schedule(runtime)
for(int j = 0; j < jmax+1 ; j++ ) //Real part
{
    RHStemp0 = value_type(0.0);
    RHStemp1 = value_type(0.0);
    for (int k = 0; k< jmax+1 ;k++)
    {
        for (int l = max(0,j+k-jmax); l < 1 + min(jmax,j+k);l++)
        {
            RHStemp0 = RHStemp0 + S[j*SIZE_S*SIZE_S + k*SIZE_S  + l]*(-A[k+jmax+1]*A[l]*A[j+k-l] + A[k]*A[l+jmax+1]*A[j+k-l]
                                                  + A[k]*A[l]*A[j+k-l+jmax+1] +A[k+jmax+1]*A[l+jmax+1]*A[j+k-l+jmax+1]);
            RHStemp1 = RHStemp1 + S[j*SIZE_S*SIZE_S + k*SIZE_S  + l]*(A[k]*A[l]*A[j+k-l] - A[k]*A[l+jmax+1]*A[j+k-l+jmax+1]
                                                  + A[k+jmax+1]*A[l]*A[j+k-l+jmax+1] +A[k+jmax+1]*A[l+jmax+1]*A[j+k-l]);
        }
    }
    dAdt[j] = (-1/(value_type((2*(2*j+3)))))*RHStemp0;
    dAdt[j+jmax+1] = (1/(value_type((2*(2*j+3)))))*RHStemp1;
}

int main()
{
const state_type initial = loadInitialData();    //Initial condition
omp_set_num_threads(jmax+1);
int chunk_size = jmax/omp_get_max_threads();
omp_set_schedule( omp_sched_dynamic, chunk_size );

//I define my controlled error steppers
typedef runge_kutta_fehlberg78< state_type , value_type ,
                  state_type , value_type,openmp_range_algebra> error_stepper_type;

typedef controlled_runge_kutta< error_stepper_type > controlled_stepper_type;
controlled_stepper_type controlled_stepper;

int steps =  integrate_adaptive( controlled_stepper  ,rhs ,
                   initial, TINITIAL  , TFINAL,INITIAL_STEP ,  push_back_state_and_time( A_vec , times ) );

}

我没有显示所有变量的定义，但我怀疑它们是问题所在，因为如果我只是从 error_stepper_type 的定义中删除 openmp_range_algebra 选项，它就可以正常工作。如果我使用具有恒定步进器大小的 openmp_range_algebra，这也很好用，例如 4 阶的 Runge Kutta。

但是，使用这段代码我得到以下错误:

invalid conversion from 'boost::range_iterator<const boost::array<double, 22ull>, void>::type {aka const double*}' to 'boost::range_iterator<boost::array<double, 22ull>, void>::type {aka double*}' [-fpermissive]|

所以看起来我试图分配一些常量。此错误出现在文件 openmp_range_algebra.hpp 中，代码如下:

template< class S >
static typename norm_result_type< S >::type norm_inf( const S &s )
{
    using std::max;
    using std::abs;
    typedef typename norm_result_type< S >::type result_type;
    result_type init = static_cast< result_type >( 0 );
    const size_t len = boost::size(s);
    typename boost::range_iterator<S>::type beg = boost::begin(s);
    #pragma omp parallel for reduction(max: init) schedule(dynamic)
    for( size_t i = 0 ; i < len ; ++i )
        init = max( init , abs( beg[i] ) );
    return init;
}

我希望我已经足够清楚了，我只是希望能够在我的并行化代码中使用自适应步进器。

非常感谢您的帮助。

Answer 1

这是 odeint 中的一个错误，我已经在 github 上提交了它:https://github.com/headmyshoulder/odeint-v2/issues/166 我会尽快修复它。 感谢发帖。

编辑:已修复，您的程序现在应该可以编译了。另外，我更改了一个示例以使用自适应步进器: https://github.com/headmyshoulder/odeint-v2/blob/master/examples/openmp/lorenz_ensemble_simple.cpp