 
 }
 
-cusp::ell_matrix<double, int, cusp::device_memory> interpolation( thrust::device_vector<double>& x, const Grid1d<double>& g  )
+} //namespace detail 
+
+cusp::ell_matrix<double, int, cusp::device_memory> interpolation( const thrust::device_vector<double>& x, const Grid1d<double>& g  )
 {
-    cusp::ell_matrix<int, double, cusp::device_memory> A( x.size(), x.size(), g.size(), g.n());
+    //allocate ell matrix storage
+    cusp::ell_matrix<int, double, cusp::device_memory> A( x.size(), g.size(), x.size()*g.n(), g.n());
 
+    //set up kernel parameters
+    const size_t BLOCK_SIZE = 256;
+    const size_t MAX_BLOCKS = cusp::system::cuda::detail::max_active_blocks( interpolation_kernal<BLOCK_SIZE>, BLOCK_SIZE, (size_t) 0  );
+    const size_t NUM_BLOCKS = std::min<size_t>( 
+            MAX_BLOCKS, 
+            cusp::system::cuda::DIVIDE_INTO( A.num_rows, BLOCK_SIZE));
+    const int pitch = A.column_indices.pitch;
+    const int* Aj = thrust::raw_pointer_cast(&A.column_indices(0,0));
+    const double * Av = thrust::raw_pointer_cast(&A.values(0,0));
+
+    //compute normalized x values and cell numbers
     thrust::device_vector<double> xn(x.size()), cellnh(x.size());
     dg::blas1::transform( x, xn, dg::PLUS<double>( -g.x0()));
     dg::blas1::scal( xn, 1./g.h());
     thrust::device_vector<int> celln( x.size());
     thrust::transform( xn.begin(), xn.end(), celln.begin(), dg::FLOOR());
     thrust::transform( celln.begin(), celln.end(),cellnh.begin(), dg::PLUS<double>(0.5));
+    const int* celln_ptr thrust::raw_pointer_cast( &celln[0]);
+    const double* xn_ptr = thrust::raw_pointer_cast( &xn[0]);
+    //xn = 2*xn - 2*(celln+0.5)
+    thrust::device_vector<double> forward(std::vector<double> ( g.dlt().forward()));
+    const double forward_ptr = thrust::raw_pointer_cast( forward);
+    dg::blas1::axpby( 2., xn, -2., cellnh, xn);
+    interpolation_kernel<BLOCK_SIZE> <<<NUM_BLOCKS, BLOCK_SIZE, 0, 0>>> ( A.num_rows, g.n(), celln_ptr, xn_ptr, pitch, Aj, Av, forward_ptr);
+    return A;
+
+}
+
+cusp::ell_matrix<double, int, cusp::device_memory> interpolation( const thrust::device_vector<double>& x, const thrust::device_vector<double>& y, const Grid2d<double>& g  )
+{
+    //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());
+
+    //set up kernel parameters
+    const size_t BLOCK_SIZE = 256;
+    const size_t MAX_BLOCKS = cusp::system::cuda::detail::max_active_blocks( interpolation_kernal<BLOCK_SIZE>, BLOCK_SIZE, (size_t) 0  );
+    const size_t NUM_BLOCKS = std::min<size_t>( 
+            MAX_BLOCKS, 
+            cusp::system::cuda::DIVIDE_INTO( A.num_rows, BLOCK_SIZE));
+    const int pitch = A.column_indices.pitch;
+    const int* Aj = thrust::raw_pointer_cast(&A.column_indices(0,0));
+    const 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());
+    dg::blas1::scal( yn, 1./g.hy());
+    thrust::device_vector<int> celln( x.size());
+    thrust::device_vector<int> cellm( x.size());
+    thrust::transform( xn.begin(), xn.end(), celln.begin(), dg::FLOOR());
+    thrust::transform( yn.begin(), yn.end(), cellm.begin(), dg::FLOOR());
+    thrust::transform( celln.begin(), celln.end(),cellnh.begin(), dg::PLUS<double>(0.5));
+    thrust::transform( cellm.begin(), cellm.end(),cellmh.begin(), dg::PLUS<double>(0.5));
+    const int* celln_ptr thrust::raw_pointer_cast( &celln[0]);
+    const int* cellm_ptr thrust::raw_pointer_cast( &cellm[0]);
+    const double* xn_ptr = thrust::raw_pointer_cast( &xn[0]);
+    const double* xm_ptr = thrust::raw_pointer_cast( &xm[0]);
     //xn = 2*xn - 2*(celln+0.5)
     dg::blas1::axpby( 2., xn, -2., cellnh, xn);
-
-    
-
+    dg::blas1::axpby( 2., yn, -2., cellmh, yn);
+    thrust::device_vector<double> forward(std::vector<double> ( g.dlt().forward()));
+    const double forward_ptr = thrust::raw_pointer_cast( forward);
+    interpolation_kernel<BLOCK_SIZE> <<<NUM_BLOCKS, BLOCK_SIZE, 0, 0>>> ( A.num_rows, g.n(), g.Nx(), celln_ptr, cellm_ptr, xn_ptr, yn_ptr, pitch, Aj, Av, forward);
+    return A;
 
 }
-} //namespace detail 
 } //namespace create
 } //namespace dg

#pragma once




    {
        //evaluate 2d Legendre polynomials at (xn, yn)...
        double px[n], py[n], pxy[n*n];
        detail::legendre( px, &xn[row], n, forward), 
        detail::legendre( py, &yn[row], n, forward), 
        
        int offset = row;
        int col_begin = cellm[row]*Nx*n*n + celln*n;
