Использование 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)); 

} 

Я не показывать определение всех переменных, но я сомневаюсь, что они проблема, так как это отлично работает, если я просто удалю опцию openmp_range_algebra из определения error_stepper_ty физическое воспитание Это работает также хорошо, если я использую openmp_range_algebra с постоянным размером шагового, как Рунге Кутта порядка 4.

Однако, с этим кодом, я получаю следующие ошибки:

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