 #pragma once
+#include <cusp/system/device/cuda/arch.h>
+#include <cusp/system/device/cuda/utils.h>
+#include "../blas.h"
+#include "../functors.h"
 
 
 

template<size_t BLOCK_SIZE>
__launch_bounds__(BLOCK_SIZE, 1) //cuda performance hint macro, (max_threads_per_block, minBlocksPerMultiprocessor)
 __global__ void 
interpolation_kernel( const int num_rows, const int n, const int Nx, int* celln, int* cellm, const double* xn, const double * yn, const int pitch, int* Aj, double Av*)
{
    const int thread_id = blockDim.x * blockIdx.x + threadIdx.x;
    const int grid_size = gridDim.x*blockDim.x

    {
        //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;
