  */
 cusp::ell_matrix<int, double, cusp::device_memory> ell_interpolation( const thrust::device_vector<double>& x, const thrust::device_vector<double>& y, const Grid2d<double>& g  )
 {
+    Timer t;
+    t.tic();
     assert( x.size() == y.size());
     //allocate ell matrix storage
     cusp::ell_matrix<int, double, cusp::device_memory> A( x.size(), g.size(), x.size()*g.n()*g.n(), g.n()*g.n());

    double * Av = thrust::raw_pointer_cast(&A.values(0,0));

    //compute normalized x and y values and cell numbers
    thrust::device_vector<double> xn(x.size()), yn(y.size()), cellnh(x.size()), cellmh(x.size());
    dg::blas1::transform( x, xn, dg::PLUS<double>( -g.x0()));
    dg::blas1::transform( y, yn, dg::PLUS<double>( -g.y0()));
    dg::blas1::scal( xn, 1./g.hx());

    double * Av = thrust::raw_pointer_cast(&A.values(0,0));

    //compute normalized x values and cell numbers
    thrust::device_vector<double> xn(x.size()), cellxh(x.size());
    dg::blas1::transform( x, xn, dg::PLUS<double>( -g.x0()));
    dg::blas1::scal( xn, 1./g.h());
    thrust::device_vector<int> cellx( x.size());
    thrust::transform( xn.begin(), xn.end(), cellx.begin(), dg::FLOOR());
    thrust::transform( cellx.begin(), cellx.end(),cellxh.begin(), dg::PLUS<double>(0.5));
    const int* cellx_ptr = thrust::raw_pointer_cast( &cellx[0]);
    const double* xn_ptr = thrust::raw_pointer_cast( &xn[0]);
    //xn = 2*xn - 2*(cellx+0.5)
    const thrust::device_vector<double> forward(std::vector<double> ( g.dlt().forward()));
    const double* forward_ptr = thrust::raw_pointer_cast( &forward[0]);
    dg::blas1::axpby( 2., xn, -2., cellxh, xn);
    detail::interpolation_kernel1d<BLOCK_SIZE> <<<NUM_BLOCKS, BLOCK_SIZE>>> ( A.num_rows, g.n(), cellx_ptr, xn_ptr, pitch, Aj, Av, forward_ptr);
    return A;

}
