/*
 Copyright 2011-2012 Karsten Ahnert
 Copyright 2013 Mario Mulansky

 Distributed under the Boost Software License, Version 1.0.
 (See accompanying file LICENSE_1_0.txt or
 copy at http://www.boost.org/LICENSE_1_0.txt)
 */


/*
 * Solves many relaxation equations dxdt = - a * x in parallel and for different values of a.
 * The relaxation equations are completely uncoupled.
 */

#include <thrust/device_vector.h>

#include <boost/ref.hpp>

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


using namespace std;
using namespace boost::numeric::odeint;

// change to float if your GPU does not support doubles
typedef double value_type;
typedef thrust::device_vector< value_type > state_type;
typedef runge_kutta4< state_type , value_type , state_type , value_type > stepper_type;

struct relaxation
{
    struct relaxation_functor
    {
        template< class T >
        __host__ __device__
        void operator()( T t ) const
        {
            // unpack the parameter we want to vary and the Lorenz variables
            value_type a = thrust::get< 1 >( t );
            value_type x = thrust::get< 0 >( t );
            thrust::get< 2 >( t ) = -a * x;
        }
    };

    relaxation( size_t N , const state_type &a )
    : m_N( N ) , m_a( a ) { }

    void operator()(  const state_type &x , state_type &dxdt , value_type t ) const
    {
        thrust::for_each(
            thrust::make_zip_iterator( thrust::make_tuple( x.begin() , m_a.begin() , dxdt.begin() ) ) ,
            thrust::make_zip_iterator( thrust::make_tuple( x.end() , m_a.end() , dxdt.end() ) ) ,
            relaxation_functor() );
    }

    size_t m_N;
    const state_type &m_a;
};

const size_t N = 1024 * 1024;
const value_type dt = 0.01;

int main( int arc , char* argv[] )
{
    // initialize the relaxation constants a
    vector< value_type > a_host( N );
    for( size_t i=0 ; i<N ; ++i ) a_host[i] = drand48();
    state_type a = a_host;

    // initialize the intial state x
    state_type x( N );
    thrust::fill( x.begin() , x.end() , 1.0 );

    // integrate
    relaxation relax( N , a );
    integrate_const( stepper_type() , boost::ref( relax ) , x , 0.0 , 10.0 , dt );

    return 0;
}