结合推力和助推器的颂歌;不正确的指针传递被认为会导致段错误

Question

我试图修改lorenz_parameters。来自boost::odeint的一个例子，用于解决许多初始条件的初值问题。程序段错误，我认为这是由于我错误地将指针传递给integrate_adaptive()的调用引起的，因为当我试图打印出用于存储导数的向量的大小时，它是0。

作为参考，这里是示例的文档，以及相同示例的代码(由于某种原因，原始文档中指向代码的链接不起作用)。

我在下面附上了程序的完整源代码。有一条标记为UNEXPECTED OUTPUT!的线，我试图在这里打印导数向量的大小，结果显示为0，它应该与状态向量(512)的大小相同。还有一条线标记LINE I DO NOT UNDERSTAND，这是我认为我可能做错了什么。这一行是对integrate_adaptive()调用的一部分，我认为这个参数是(最终/通过odeint包的方式)传递给parallel_initial_condition_problem中的operator()函数的。但我试着做同样的例子，但它不起作用。任何帮助都非常感谢!

#include <iostream>
#include <vector>
#include <thrust/device_vector.h>
#include <boost/numeric/odeint.hpp>
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <typeinfo>
typedef double value_type;
//change this to device_vector< ... > of you want to run on GPU
typedef thrust::host_vector< value_type > state_type;
using namespace std;
using namespace boost::numeric::odeint;
const size_t D = 16; // how many values of each initial condition to try
const size_t N_ICS = D * D; // number of initial conditions
const size_t N_VARS = 2; // number of variables in the system
// time constants
const value_type t_0 = 0.0;
const value_type t_final = 0.1;    
const value_type dt = 0.01;
struct parallel_initial_condition_problem {
  struct functor
  {
    template< class T >
    __host__ __device__
    void operator()( T t ) const
    {
      const value_type k_b = 0.45/6.0;
      const value_type k_f = 1.0;
      const value_type k_d = 1.0;
      // unpack the variables (and the parameters we want to vary, if any)
      value_type f = thrust::get< 0 >( t );
      value_type a = thrust::get< 1 >( t );
      // value_type df = thrust::get< 2 >( t ); // this also causes a segfault
      // value_type da = thrust::get< 3 >( t );
      // calculate the derivative
      thrust::get< 2 >( t ) = -0.01; // minimal value for debugging
      thrust::get< 3 >( t ) = -0.01; 
    }
  };
  parallel_initial_condition_problem( size_t N )
  : m_N( N ) { }
  template< class State , class Deriv >
  void operator()(const State &x , Deriv &dxdt , value_type t ) const
  {
    //// Printing out debugging info (here is the error..but why?)
    cout << "INSIDE operator()" << endl;
    cout << "t: " << t << endl;    
    cout << "size of x: t" << boost::end(x) - boost::begin(x) << endl;
    cout << "size of dxdt: t" << boost::end(dxdt) - boost::begin(dxdt) << endl; // UNEXPECTED OUTPUT!!
    cout << endl;
    thrust::for_each(
      thrust::make_zip_iterator( 
        thrust::make_tuple(
               boost::begin( x ) + 0*m_N,
               boost::begin( x ) + 1*m_N ,
               boost::begin( dxdt ) + 0*m_N,
               boost::begin( dxdt ) + 1*m_N
               )
                 ) ,
      thrust::make_zip_iterator( 
        thrust::make_tuple(
               boost::begin( x ) + 1*m_N,
               boost::begin( x ) + 2*m_N ,
               boost::begin( dxdt ) + 3*m_N,
               boost::begin( dxdt ) + 4*m_N
               )
                 ) ,
      functor() );
  }
  size_t m_N;
};
int main(int /* argc */ , char** /* argv */ ) {
  // x stores both (initial) state (x) [0:N_VARS*N_ICS] and dxdt [N_VARS*N_ICS:END]
  state_type x( 2 * N_VARS * N_ICS ); 
  // DUMMY INITIAL CONDITIONS FOR DEBUGGING PURPOSES
  // fill each value with its index
  for(int i=0;i<x.size();i++) {
    x[i] = 0.01+i;
  }
  // system function
  parallel_initial_condition_problem init_con_solver( N_ICS );
  //////////////////////////////// Debugging info: (all seems correct here)
  cout << "POINTERS"  << endl;
  cout << "size of x: " << boost::end(x) - boost::begin(x) << endl;
  cout << "x    :" << &(*x.begin()) << endl;
  cout << "dxdt :" << &(*x.begin()) + N_VARS*N_ICS << endl;
  cout << "end  :" << &(*x.end()) << endl;
  cout << endl;
  /////////////////////////// STATE_T      VALUE_T      DERIV_T      TIME_T
  typedef runge_kutta_dopri5< state_type , value_type , state_type , value_type > stepper_type;
  const value_type abs_err = 1.0e-6;
  const value_type rel_err = 1.0e-6;
  integrate_adaptive( make_controlled( abs_err , rel_err , stepper_type() ),
              init_con_solver ,
              std::make_pair( x.begin() , x.begin() + N_VARS*N_ICS  ), /// LINE I DO NOT UNDERSTAND
              t_0, t_final, dt );
}

运行这个命令的输出是:

$ g++ main.c && ./a.out
POINTERS
size of x: 1024
x    :0x9a2010
dxdt :0x9a3010
end  :0x9a4010
INSIDE operator()
t: 0
size of x:  512
size of dxdt:   0
Segmentation fault (core dumped)

Answer 1

我认为有一些包括缺失。试着与

#include <boost/numeric/odeint/external/thrust/thrust.hpp>

实际上我想知道，为什么你的代码甚至可以编译。如果您使用boost(1.55)的当前odeint版本，您还需要指定代数和操作类型:

typedef runge_kutta_dopri5<
    state_type ,
    value_type ,
    state_type ,
    value_type ,
    thrust_algebra ,
    thrust_operations > stepper_type;

这应该在下一个boost版本中改变。对于gihub版本的odeint，您的代码应该没问题。