text: /******************************************************************************
   * $Id: gdalwarpkernel_opencl.cpp 36763 2016-12-09 22:10:55Z rouault $
   *
   * Project:  OpenCL Image Reprojector
   * Purpose:  Implementation of the GDALWarpKernel reprojector in OpenCL.
   * Author:   Seth Price, seth@pricepages.org
   *
   ******************************************************************************
   * Copyright (c) 2010, Seth Price <seth@pricepages.org>
   * Copyright (c) 2010-2012, Even Rouault <even dot rouault at mines-paris dot org>
   *
   * Permission is hereby granted, free of charge, to any person obtaining a
   * copy of this software and associated documentation files (the "Software"),
   * to deal in the Software without restriction, including without limitation
   * the rights to use, copy, modify, merge, publish, distribute, sublicense,
   * and/or sell copies of the Software, and to permit persons to whom the
   * Software is furnished to do so, subject to the following conditions:
   *
   * The above copyright notice and this permission notice shall be included
   * in all copies or substantial portions of the Software.
   *
   * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
   * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
   * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
   * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
   * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
   * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
   * DEALINGS IN THE SOFTWARE.
   ****************************************************************************/
  
  #if defined(HAVE_OPENCL)
  
  /* The following line may be uncommented for increased debugging traces and profiling */
  /* #define DEBUG_OPENCL 1 */
  
  #include <assert.h>
  #include <stdlib.h>
  #include <stdio.h>
  #include <limits.h>
  #include <float.h>
  #include "cpl_string.h"
  #include "gdalwarpkernel_opencl.h"
  
  CPL_CVSID("$Id: gdalwarpkernel_opencl.cpp 36763 2016-12-09 22:10:55Z rouault $");
  
  #define handleErr(err) if((err) != CL_SUCCESS) { \
      CPLError(CE_Failure, CPLE_AppDefined, "Error at file %s line %d: %s", __FILE__, __LINE__, getCLErrorString(err)); \
      return err; \
  }
  
  #define handleErrRetNULL(err) if((err) != CL_SUCCESS) { \
      (*clErr) = err; \
      CPLError(CE_Failure, CPLE_AppDefined, "Error at file %s line %d: %s", __FILE__, __LINE__, getCLErrorString(err)); \
      return NULL; \
  }
  
  #define handleErrGoto(err, goto_label) if((err) != CL_SUCCESS) { \
      (*clErr) = err; \
      CPLError(CE_Failure, CPLE_AppDefined, "Error at file %s line %d: %s", __FILE__, __LINE__, getCLErrorString(err)); \
      goto goto_label; \
  }
  
  #define freeCLMem(clMem, fallBackMem)  do { \
      if ((clMem) != NULL) { \
          handleErr(err = clReleaseMemObject(clMem)); \
          clMem = NULL; \
          fallBackMem = NULL; \
      } else if ((fallBackMem) != NULL) { \
          CPLFree(fallBackMem); \
          fallBackMem = NULL; \
      } \
  } while( false )
  
  static const char* getCLErrorString(cl_int err)
  {
      switch (err)
      {
          case CL_SUCCESS:
              return("CL_SUCCESS");
              break;
          case CL_DEVICE_NOT_FOUND:
              return("CL_DEVICE_NOT_FOUND");
              break;
          case CL_DEVICE_NOT_AVAILABLE:
              return("CL_DEVICE_NOT_AVAILABLE");
              break;
          case CL_COMPILER_NOT_AVAILABLE:
              return("CL_COMPILER_NOT_AVAILABLE");
              break;
          case CL_MEM_OBJECT_ALLOCATION_FAILURE:
              return("CL_MEM_OBJECT_ALLOCATION_FAILURE");
              break;
          case CL_OUT_OF_RESOURCES:
              return("CL_OUT_OF_RESOURCES");
              break;
          case CL_OUT_OF_HOST_MEMORY:
              return("CL_OUT_OF_HOST_MEMORY");
              break;
          case CL_PROFILING_INFO_NOT_AVAILABLE:
              return("CL_PROFILING_INFO_NOT_AVAILABLE");
              break;
          case CL_MEM_COPY_OVERLAP:
              return("CL_MEM_COPY_OVERLAP");
              break;
          case CL_IMAGE_FORMAT_MISMATCH:
              return("CL_IMAGE_FORMAT_MISMATCH");
              break;
          case CL_IMAGE_FORMAT_NOT_SUPPORTED:
              return("CL_IMAGE_FORMAT_NOT_SUPPORTED");
              break;
          case CL_BUILD_PROGRAM_FAILURE:
              return("CL_BUILD_PROGRAM_FAILURE");
              break;
          case CL_MAP_FAILURE:
              return("CL_MAP_FAILURE");
              break;
          case CL_INVALID_VALUE:
              return("CL_INVALID_VALUE");
              break;
          case CL_INVALID_DEVICE_TYPE:
              return("CL_INVALID_DEVICE_TYPE");
              break;
          case CL_INVALID_PLATFORM:
              return("CL_INVALID_PLATFORM");
              break;
          case CL_INVALID_DEVICE:
              return("CL_INVALID_DEVICE");
              break;
          case CL_INVALID_CONTEXT:
              return("CL_INVALID_CONTEXT");
              break;
          case CL_INVALID_QUEUE_PROPERTIES:
              return("CL_INVALID_QUEUE_PROPERTIES");
              break;
          case CL_INVALID_COMMAND_QUEUE:
              return("CL_INVALID_COMMAND_QUEUE");
              break;
          case CL_INVALID_HOST_PTR:
              return("CL_INVALID_HOST_PTR");
              break;
          case CL_INVALID_MEM_OBJECT:
              return("CL_INVALID_MEM_OBJECT");
              break;
          case CL_INVALID_IMAGE_FORMAT_DESCRIPTOR:
              return("CL_INVALID_IMAGE_FORMAT_DESCRIPTOR");
              break;
          case CL_INVALID_IMAGE_SIZE:
              return("CL_INVALID_IMAGE_SIZE");
              break;
          case CL_INVALID_SAMPLER:
              return("CL_INVALID_SAMPLER");
              break;
          case CL_INVALID_BINARY:
              return("CL_INVALID_BINARY");
              break;
          case CL_INVALID_BUILD_OPTIONS:
              return("CL_INVALID_BUILD_OPTIONS");
              break;
          case CL_INVALID_PROGRAM:
              return("CL_INVALID_PROGRAM");
              break;
          case CL_INVALID_PROGRAM_EXECUTABLE:
              return("CL_INVALID_PROGRAM_EXECUTABLE");
              break;
          case CL_INVALID_KERNEL_NAME:
              return("CL_INVALID_KERNEL_NAME");
              break;
          case CL_INVALID_KERNEL_DEFINITION:
              return("CL_INVALID_KERNEL_DEFINITION");
              break;
          case CL_INVALID_KERNEL:
              return("CL_INVALID_KERNEL");
              break;
          case CL_INVALID_ARG_INDEX:
              return("CL_INVALID_ARG_INDEX");
              break;
          case CL_INVALID_ARG_VALUE:
              return("CL_INVALID_ARG_VALUE");
              break;
          case CL_INVALID_ARG_SIZE:
              return("CL_INVALID_ARG_SIZE");
              break;
          case CL_INVALID_KERNEL_ARGS:
              return("CL_INVALID_KERNEL_ARGS");
              break;
          case CL_INVALID_WORK_DIMENSION:
              return("CL_INVALID_WORK_DIMENSION");
              break;
          case CL_INVALID_WORK_GROUP_SIZE:
              return("CL_INVALID_WORK_GROUP_SIZE");
              break;
          case CL_INVALID_WORK_ITEM_SIZE:
              return("CL_INVALID_WORK_ITEM_SIZE");
              break;
          case CL_INVALID_GLOBAL_OFFSET:
              return("CL_INVALID_GLOBAL_OFFSET");
              break;
          case CL_INVALID_EVENT_WAIT_LIST:
              return("CL_INVALID_EVENT_WAIT_LIST");
              break;
          case CL_INVALID_EVENT:
              return("CL_INVALID_EVENT");
              break;
          case CL_INVALID_OPERATION:
              return("CL_INVALID_OPERATION");
              break;
          case CL_INVALID_GL_OBJECT:
              return("CL_INVALID_GL_OBJECT");
              break;
          case CL_INVALID_BUFFER_SIZE:
              return("CL_INVALID_BUFFER_SIZE");
              break;
          case CL_INVALID_MIP_LEVEL:
              return("CL_INVALID_MIP_LEVEL");
              break;
          case CL_INVALID_GLOBAL_WORK_SIZE:
              return("CL_INVALID_GLOBAL_WORK_SIZE");
              break;
      }
  
      return "unknown_error";
  }
  
  static const char* getCLDataTypeString( cl_channel_type dataType )
  {
      switch( dataType )
      {
          case CL_SNORM_INT8: return "CL_SNORM_INT8";
          case CL_SNORM_INT16: return "CL_SNORM_INT16";
          case CL_UNORM_INT8: return "CL_UNORM_INT8";
          case CL_UNORM_INT16: return "CL_UNORM_INT16";
  #if 0
          case CL_UNORM_SHORT_565: return "CL_UNORM_SHORT_565";
          case CL_UNORM_SHORT_555: return "CL_UNORM_SHORT_555";
          case CL_UNORM_INT_101010: return "CL_UNORM_INT_101010";
          case CL_SIGNED_INT8: return "CL_SIGNED_INT8";
          case CL_SIGNED_INT16: return "CL_SIGNED_INT16";
          case CL_SIGNED_INT32: return "CL_SIGNED_INT32";
          case CL_UNSIGNED_INT8: return "CL_UNSIGNED_INT8";
          case CL_UNSIGNED_INT16: return "CL_UNSIGNED_INT16";
          case CL_UNSIGNED_INT32: return "CL_UNSIGNED_INT32";
          case CL_HALF_FLOAT: return "CL_HALF_FLOAT";
  #endif
          case CL_FLOAT: return "CL_FLOAT";
          default: return "unknown";
      }
  }
  
  /*
   Finds an appropriate OpenCL device. For debugging, it's
   always easier to use CL_DEVICE_TYPE_CPU because then */ /*ok*/ /*printf() can be called
   from the kernel. If debugging is on, we can print the name and stats about the
   device we're using.
   */
  static cl_device_id get_device(OCLVendor *peVendor)
  {
      cl_int err = 0;
      size_t returned_size = 0;
      cl_char vendor_name[1024] = {0};
      cl_char device_name[1024] = {0};
  
      cl_platform_id platforms[10];
      cl_uint num_platforms;
      cl_uint i;
      cl_device_id preferred_device_id = NULL;
      int preferred_is_gpu = FALSE;
  
      static bool gbBuggyOpenCL = false;
      if( gbBuggyOpenCL )
          return NULL;
      try
      {
          err = clGetPlatformIDs( 10, platforms, &num_platforms );
          if( err != CL_SUCCESS || num_platforms == 0 )
              return NULL;
      }
      catch( ... )
      {
          gbBuggyOpenCL = true;   
          CPLDebug("OpenCL", "clGetPlatformIDs() threw a C++ exception");
          // This should normally not happen. But that does happen with
          // intel-opencl 0r2.0-54426 when run under xvfb-run 
          return NULL;
      }
      
      bool bUseOpenCLCPU = CPLTestBool( CPLGetConfigOption("OPENCL_USE_CPU", "FALSE") );
  
      // In case we have several implementations, pick up the non Intel one by
      // default, unless the PREFERRED_OPENCL_VENDOR config option is specified.
      for( i=0; i<num_platforms;i++)
      {
          cl_device_id device = NULL;
          const char* pszBlacklistedVendor;
          const char* pszPreferredVendor;
          int is_gpu;
  
          // Find the GPU CL device, this is what we really want
          // If there is no GPU device is CL capable, fall back to CPU
          if( bUseOpenCLCPU )
              err = CL_DEVICE_NOT_FOUND;
          else
              err = clGetDeviceIDs(platforms[i], CL_DEVICE_TYPE_GPU, 1,
                                   &device, NULL);
          is_gpu = (err == CL_SUCCESS);
          if (err != CL_SUCCESS)
          {
              // Find the CPU CL device, as a fallback
              err = clGetDeviceIDs(platforms[i], CL_DEVICE_TYPE_CPU, 1,
                                   &device, NULL);
              if( err != CL_SUCCESS || device == 0 )
                  continue;
          }
  
          // Get some information about the returned device
          err = clGetDeviceInfo(device, CL_DEVICE_VENDOR, sizeof(vendor_name),
                              vendor_name, &returned_size);
          err |= clGetDeviceInfo(device, CL_DEVICE_NAME, sizeof(device_name),
                              device_name, &returned_size);
          assert(err == CL_SUCCESS);
  
          if( num_platforms > 1 )
              CPLDebug( "OpenCL", "Found vendor='%s' / device='%s' (%s implementation).",
                        vendor_name, device_name, (is_gpu) ? "GPU" : "CPU");
          
          pszBlacklistedVendor = CPLGetConfigOption("BLACKLISTED_OPENCL_VENDOR", NULL);
          if( pszBlacklistedVendor &&
              EQUAL( (const char*)vendor_name, pszBlacklistedVendor ) ) 
          {
              CPLDebug("OpenCL", "Blacklisted vendor='%s' / device='%s' implementation skipped",
                       vendor_name, device_name);
              continue;
          }
              
          if( preferred_device_id == NULL || (is_gpu && !preferred_is_gpu) )
          {
              preferred_device_id = device;
              preferred_is_gpu = is_gpu;
          }
              
          pszPreferredVendor = CPLGetConfigOption("PREFERRED_OPENCL_VENDOR", NULL);
          if( pszPreferredVendor )
          {
              if( EQUAL( (const char*)vendor_name, pszPreferredVendor ) )
              {
                  preferred_device_id = device;
                  preferred_is_gpu = is_gpu;
                  break;
              }
          }
          else if( is_gpu && !STARTS_WITH((const char*)vendor_name, "Intel") )
          {
              preferred_device_id = device;
              preferred_is_gpu = is_gpu;
              break;
          }
      }
      if( preferred_device_id == NULL )
      {
          CPLDebug("OpenCL", "No implementation found");
          return NULL;
      }
      
      err = clGetDeviceInfo(preferred_device_id, CL_DEVICE_VENDOR, sizeof(vendor_name),
                              vendor_name, &returned_size);
      err |= clGetDeviceInfo(preferred_device_id, CL_DEVICE_NAME, sizeof(device_name),
                              device_name, &returned_size);
      CPLDebug( "OpenCL", "Connected to vendor='%s' / device='%s' (%s implementation).",
                vendor_name, device_name, (preferred_is_gpu) ? "GPU" : "CPU");
  
      if (STARTS_WITH((const char*)vendor_name, "Advanced Micro Devices"))
          *peVendor = VENDOR_AMD;
      else if (STARTS_WITH((const char*)vendor_name, "Intel"))
          *peVendor = VENDOR_INTEL;
      else
          *peVendor = VENDOR_OTHER;
  
      return preferred_device_id;
  }
  
  /*
   Given that not all OpenCL devices support the same image formats, we need to
   make do with what we have. This leads to wasted space, but as OpenCL matures
   I hope it'll get better.
   */
  static cl_int set_supported_formats(struct oclWarper *warper,
                               cl_channel_order minOrderSize,
                               cl_channel_order *chosenOrder,
                               unsigned int *chosenSize,
                               cl_channel_type dataType )
  {
      cl_image_format *fmtBuf = (cl_image_format *)calloc(256, sizeof(cl_image_format));
      cl_uint numRet;
      cl_uint i;
      cl_uint extraSpace = 9999;
      cl_int err;
      int bFound = FALSE;
  
      //Find what we *can* handle
      handleErr(err = clGetSupportedImageFormats(warper->context,
                                                 CL_MEM_READ_ONLY,
                                                 CL_MEM_OBJECT_IMAGE2D,
                                                 256, fmtBuf, &numRet));
      for (i = 0; i < numRet; ++i) {
          cl_channel_order thisOrderSize = 0;
          switch (fmtBuf[i].image_channel_order)
          {
              //Only support formats which use the channels in order (x,y,z,w)
            case CL_R:
            case CL_INTENSITY:
            case CL_LUMINANCE:
              thisOrderSize = 1;
              break;
            case CL_RG:
              thisOrderSize = 2;
              break;
            case CL_RGB:
              thisOrderSize = 3;
              break;
            case CL_RGBA:
              thisOrderSize = 4;
              break;
          }
  
          //Choose an order with the least wasted space
          if (fmtBuf[i].image_channel_data_type == dataType &&
              minOrderSize <= thisOrderSize &&
              extraSpace > thisOrderSize - minOrderSize ) {
  
              //Set the vector size, order, & remember wasted space
              (*chosenSize) = thisOrderSize;
              (*chosenOrder) = fmtBuf[i].image_channel_order;
              extraSpace = thisOrderSize - minOrderSize;
              bFound = TRUE;
          }
      }
  
      free(fmtBuf);
  
      if( !bFound )
      {
          CPLDebug("OpenCL",
                   "Cannot find supported format for dataType = %s and minOrderSize = %d",
                   getCLDataTypeString(dataType), (int)minOrderSize);
      }
      return (bFound) ? CL_SUCCESS : CL_INVALID_OPERATION;
  }
  
  /*
   Allocate some pinned memory that we can use as an intermediate buffer. We're
   using the pinned memory to assemble the data before transferring it to the
   device. The reason we're using pinned RAM is because the transfer speed from
   host RAM to device RAM is faster than non-pinned. The disadvantage is that
   pinned RAM is a scarce OS resource. I'm making the assumption that the user
   has as much pinned host RAM available as total device RAM because device RAM
   tends to be similarly scarce. However, if the pinned memory fails we fall back
   to using a regular memory allocation.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int alloc_pinned_mem(struct oclWarper *warper, int imgNum, size_t dataSz,
                          void **wrkPtr, cl_mem *wrkCL)
  {
      cl_int err = CL_SUCCESS;
      wrkCL[imgNum] = clCreateBuffer(warper->context,
                                     CL_MEM_READ_ONLY | CL_MEM_ALLOC_HOST_PTR,
                                     dataSz, NULL, &err);
  
      if (err == CL_SUCCESS) {
          wrkPtr[imgNum] = (void *)clEnqueueMapBuffer(warper->queue, wrkCL[imgNum],
                                                      CL_FALSE, CL_MAP_WRITE,
                                                      0, dataSz, 0, NULL, NULL, &err);
          handleErr(err);
      } else {
          wrkCL[imgNum] = NULL;
  #ifdef DEBUG_OPENCL
          CPLDebug("OpenCL", "Using fallback non-pinned memory!");
  #endif
          //Fallback to regular allocation
          wrkPtr[imgNum] = (void *)VSI_MALLOC_VERBOSE(dataSz);
  
          if (wrkPtr[imgNum] == NULL)
              handleErr(err = CL_OUT_OF_HOST_MEMORY);
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Allocates the working host memory for all bands of the image in the warper
   structure. This includes both the source image buffers and the destination
   buffers. This memory is located on the host, so we can assemble the image.
   Reasons for buffering it like this include reading each row from disk and
   de-interleaving bands and parts of bands. Then they can be copied to the device
   as a single operation fit for use as an OpenCL memory object.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int alloc_working_arr(struct oclWarper *warper,
                           size_t ptrSz, size_t dataSz, CPL_UNUSED size_t *fmtSz)
  {
      cl_int err = CL_SUCCESS;
      int i, b;
      size_t srcDataSz1, dstDataSz1, srcDataSz4, dstDataSz4;
      const int numBands = warper->numBands;
  
      //Find the best channel order for this format
      err = set_supported_formats(warper, 1,
                                  &(warper->imgChOrder1), &(warper->imgChSize1),
                                  warper->imageFormat);
      handleErr(err);
      if(warper->useVec) {
          err = set_supported_formats(warper, 4,
                                      &(warper->imgChOrder4), &(warper->imgChSize4),
                                      warper->imageFormat);
          handleErr(err);
      }
  
      //Alloc space for pointers to the main image data
      warper->realWork.v = (void **)VSI_CALLOC_VERBOSE(ptrSz, warper->numImages);
      warper->dstRealWork.v = (void **)VSI_CALLOC_VERBOSE(ptrSz, warper->numImages);
      if (warper->realWork.v == NULL || warper->dstRealWork.v == NULL)
          handleErr(err = CL_OUT_OF_HOST_MEMORY);
  
      if (warper->imagWorkCL != NULL) {
          //Alloc space for pointers to the extra channel, if it exists
          warper->imagWork.v = (void **)VSI_CALLOC_VERBOSE(ptrSz, warper->numImages);
          warper->dstImagWork.v = (void **)VSI_CALLOC_VERBOSE(ptrSz, warper->numImages);
          if (warper->imagWork.v == NULL || warper->dstImagWork.v == NULL)
              handleErr(err = CL_OUT_OF_HOST_MEMORY);
      } else {
          warper->imagWork.v = NULL;
          warper->dstImagWork.v = NULL;
      }
  
      //Calc the sizes we need
      srcDataSz1 = dataSz * warper->srcWidth * warper->srcHeight * warper->imgChSize1;
      dstDataSz1 = dataSz * warper->dstWidth * warper->dstHeight;
      srcDataSz4 = dataSz * warper->srcWidth * warper->srcHeight * warper->imgChSize4;
      dstDataSz4 = dataSz * warper->dstWidth * warper->dstHeight * 4;
  
      //Allocate pinned memory for each band's image
      for (b = 0, i = 0; b < numBands && i < warper->numImages; ++i) {
          if(warper->useVec && b < numBands - numBands % 4) {
              handleErr(err = alloc_pinned_mem(warper, i, srcDataSz4,
                                               warper->realWork.v,
                                               warper->realWorkCL));
  
              handleErr(err = alloc_pinned_mem(warper, i, dstDataSz4,
                                               warper->dstRealWork.v,
                                               warper->dstRealWorkCL));
              b += 4;
          } else {
              handleErr(err = alloc_pinned_mem(warper, i, srcDataSz1,
                                               warper->realWork.v,
                                               warper->realWorkCL));
  
              handleErr(err = alloc_pinned_mem(warper, i, dstDataSz1,
                                               warper->dstRealWork.v,
                                               warper->dstRealWorkCL));
              ++b;
          }
      }
  
      if (warper->imagWorkCL != NULL) {
          //Allocate pinned memory for each band's extra channel, if exists
          for (b = 0, i = 0; b < numBands && i < warper->numImages; ++i) {
              if(warper->useVec && b < numBands - numBands % 4) {
                  handleErr(err = alloc_pinned_mem(warper, i, srcDataSz4,
                                                   warper->imagWork.v,
                                                   warper->imagWorkCL));
  
                  handleErr(err = alloc_pinned_mem(warper, i, dstDataSz4,
                                                   warper->dstImagWork.v,
                                                   warper->dstImagWorkCL));
                  b += 4;
              } else {
                  handleErr(err = alloc_pinned_mem(warper, i, srcDataSz1,
                                                   warper->imagWork.v,
                                                   warper->imagWorkCL));
  
                  handleErr(err = alloc_pinned_mem(warper, i, dstDataSz1,
                                                   warper->dstImagWork.v,
                                                   warper->dstImagWorkCL));
                  ++b;
              }
          }
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Assemble and create the kernel. For optimization, portability, and
   implementation limitation reasons, the program is actually assembled from
   several strings, then compiled with as many invariants as possible defined by
   the preprocessor. There is also quite a bit of error-catching code in here
   because the kernel is where many bugs show up.
  
   Returns CL_SUCCESS on success and other CL_* errors in the error buffer when
   something goes wrong.
   */
  static
  cl_kernel get_kernel(struct oclWarper *warper, char useVec,
                       double dfXScale, double dfYScale, double dfXFilter, double dfYFilter,
                       int nXRadius, int nYRadius, int nFiltInitX, int nFiltInitY,
                       cl_int *clErr )
  {
      cl_program program;
      cl_kernel kernel;
      cl_int err = CL_SUCCESS;
  #define PROGBUF_SIZE 128000
      char *buffer = (char *)CPLCalloc(PROGBUF_SIZE, sizeof(char));
      char *progBuf = (char *)CPLCalloc(PROGBUF_SIZE, sizeof(char));
      float dstMinVal = 0.f, dstMaxVal = 0.0;
  
      const char *outType;
      const char *dUseVec = "";
      const char *dVecf = "float";
      const char *kernGenFuncs =
  // ********************* General Funcs ********************
  "void clampToDst(float fReal,\n"
                  "__global outType *dstPtr,\n"
                  "unsigned int iDstOffset,\n"
                  "__constant float *fDstNoDataReal,\n"
                  "int bandNum);\n"
  "void setPixel(__global outType *dstReal,\n"
                  "__global outType *dstImag,\n"
                  "__global float *dstDensity,\n"
                  "__global int *nDstValid,\n"
                  "__constant float *fDstNoDataReal,\n"
                  "const int bandNum,\n"
                  "vecf fDensity, vecf fReal, vecf fImag);\n"
  "int getPixel(__read_only image2d_t srcReal,\n"
                  "__read_only image2d_t srcImag,\n"
                  "__global float *fUnifiedSrcDensity,\n"
                  "__global int *nUnifiedSrcValid,\n"
                  "__constant char *useBandSrcValid,\n"
                  "__global int *nBandSrcValid,\n"
                  "const int2 iSrc,\n"
                  "int bandNum,\n"
                  "vecf *fDensity, vecf *fReal, vecf *fImag);\n"
  "int isValid(__global float *fUnifiedSrcDensity,\n"
                  "__global int *nUnifiedSrcValid,\n"
                  "float2 fSrcCoords );\n"
  "float2 getSrcCoords(__read_only image2d_t srcCoords);\n"
  
  "#ifdef USE_CLAMP_TO_DST_FLOAT\n"
  "void clampToDst(float fReal,\n"
                  "__global outType *dstPtr,\n"
                  "unsigned int iDstOffset,\n"
                  "__constant float *fDstNoDataReal,\n"
                  "int bandNum)\n"
  "{\n"
      "dstPtr[iDstOffset] = fReal;\n"
  "}\n"
  "#else\n"
  "void clampToDst(float fReal,\n"
                  "__global outType *dstPtr,\n"
                  "unsigned int iDstOffset,\n"
                  "__constant float *fDstNoDataReal,\n"
                  "int bandNum)\n"
  "{\n"
      "fReal *= dstMaxVal;\n"
  
      "if (fReal < dstMinVal)\n"
          "dstPtr[iDstOffset] = (outType)dstMinVal;\n"
      "else if (fReal > dstMaxVal)\n"
          "dstPtr[iDstOffset] = (outType)dstMaxVal;\n"
      "else\n"
          "dstPtr[iDstOffset] = (dstMinVal < 0) ? (outType)floor(fReal + 0.5f) : (outType)(fReal + 0.5f);\n"
  
      "if (useDstNoDataReal && bandNum >= 0 &&\n"
          "fDstNoDataReal[bandNum] == dstPtr[iDstOffset])\n"
      "{\n"
          "if (dstPtr[iDstOffset] == dstMinVal)\n"
              "dstPtr[iDstOffset] = dstMinVal + 1;\n"
          "else\n"
              "dstPtr[iDstOffset] = dstPtr[iDstOffset] - 1;\n"
      "}\n"
  "}\n"
  "#endif\n"
  
  "void setPixel(__global outType *dstReal,\n"
                "__global outType *dstImag,\n"
                "__global float *dstDensity,\n"
                "__global int *nDstValid,\n"
                "__constant float *fDstNoDataReal,\n"
                "const int bandNum,\n"
                "vecf fDensity, vecf fReal, vecf fImag)\n"
  "{\n"
      "unsigned int iDstOffset = get_global_id(1)*iDstWidth + get_global_id(0);\n"
  
  "#ifdef USE_VEC\n"
      "if (fDensity.x < 0.00001f && fDensity.y < 0.00001f &&\n"
          "fDensity.z < 0.00001f && fDensity.w < 0.00001f ) {\n"
  
          "fReal = 0.0f;\n"
          "fImag = 0.0f;\n"
  
      "} else if (fDensity.x < 0.9999f || fDensity.y < 0.9999f ||\n"
                 "fDensity.z < 0.9999f || fDensity.w < 0.9999f ) {\n"
          "vecf fDstReal, fDstImag;\n"
          "float fDstDensity;\n"
  
          "fDstReal.x = dstReal[iDstOffset];\n"
          "fDstReal.y = dstReal[iDstOffset+iDstHeight*iDstWidth];\n"
          "fDstReal.z = dstReal[iDstOffset+iDstHeight*iDstWidth*2];\n"
          "fDstReal.w = dstReal[iDstOffset+iDstHeight*iDstWidth*3];\n"
          "if (useImag) {\n"
              "fDstImag.x = dstImag[iDstOffset];\n"
              "fDstImag.y = dstImag[iDstOffset+iDstHeight*iDstWidth];\n"
              "fDstImag.z = dstImag[iDstOffset+iDstHeight*iDstWidth*2];\n"
              "fDstImag.w = dstImag[iDstOffset+iDstHeight*iDstWidth*3];\n"
          "}\n"
  "#else\n"
      "if (fDensity < 0.00001f) {\n"
  
          "fReal = 0.0f;\n"
          "fImag = 0.0f;\n"
  
      "} else if (fDensity < 0.9999f) {\n"
          "vecf fDstReal, fDstImag;\n"
          "float fDstDensity;\n"
  
          "fDstReal = dstReal[iDstOffset];\n"
          "if (useImag)\n"
              "fDstImag = dstImag[iDstOffset];\n"
  "#endif\n"
  
          "if (useDstDensity)\n"
              "fDstDensity = dstDensity[iDstOffset];\n"
          "else if (useDstValid &&\n"
                   "!((nDstValid[iDstOffset>>5] & (0x01 << (iDstOffset & 0x1f))) ))\n"
              "fDstDensity = 0.0f;\n"
          "else\n"
              "fDstDensity = 1.0f;\n"
  
          "vecf fDstInfluence = (1.0f - fDensity) * fDstDensity;\n"
  
          // Density should be checked for <= 0.0 & handled by the calling function
          "fReal = (fReal * fDensity + fDstReal * fDstInfluence) / (fDensity + fDstInfluence);\n"
          "if (useImag)\n"
              "fImag = (fImag * fDensity + fDstImag * fDstInfluence) / (fDensity + fDstInfluence);\n"
      "}\n"
  
  "#ifdef USE_VEC\n"
      "clampToDst(fReal.x, dstReal, iDstOffset, fDstNoDataReal, bandNum);\n"
      "clampToDst(fReal.y, dstReal, iDstOffset+iDstHeight*iDstWidth, fDstNoDataReal, bandNum);\n"
      "clampToDst(fReal.z, dstReal, iDstOffset+iDstHeight*iDstWidth*2, fDstNoDataReal, bandNum);\n"
      "clampToDst(fReal.w, dstReal, iDstOffset+iDstHeight*iDstWidth*3, fDstNoDataReal, bandNum);\n"
      "if (useImag) {\n"
          "clampToDst(fImag.x, dstImag, iDstOffset, fDstNoDataReal, -1);\n"
          "clampToDst(fImag.y, dstImag, iDstOffset+iDstHeight*iDstWidth, fDstNoDataReal, -1);\n"
          "clampToDst(fImag.z, dstImag, iDstOffset+iDstHeight*iDstWidth*2, fDstNoDataReal, -1);\n"
          "clampToDst(fImag.w, dstImag, iDstOffset+iDstHeight*iDstWidth*3, fDstNoDataReal, -1);\n"
      "}\n"
  "#else\n"
      "clampToDst(fReal, dstReal, iDstOffset, fDstNoDataReal, bandNum);\n"
      "if (useImag)\n"
          "clampToDst(fImag, dstImag, iDstOffset, fDstNoDataReal, -1);\n"
  "#endif\n"
  "}\n"
  
  "int getPixel(__read_only image2d_t srcReal,\n"
               "__read_only image2d_t srcImag,\n"
               "__global float *fUnifiedSrcDensity,\n"
               "__global int *nUnifiedSrcValid,\n"
               "__constant char *useBandSrcValid,\n"
               "__global int *nBandSrcValid,\n"
               "const int2 iSrc,\n"
               "int bandNum,\n"
               "vecf *fDensity, vecf *fReal, vecf *fImag)\n"
  "{\n"
      "int iSrcOffset = 0, iBandValidLen = 0, iSrcOffsetMask = 0;\n"
      "int bHasValid = FALSE;\n"
  
      // Clamp the src offset values if needed
      "if(useUnifiedSrcDensity | useUnifiedSrcValid | useUseBandSrcValid){\n"
          "int iSrcX = iSrc.x;\n"
          "int iSrcY = iSrc.y;\n"
  
          // Needed because the offset isn't clamped in OpenCL hardware
          "if(iSrcX < 0)\n"
              "iSrcX = 0;\n"
          "else if(iSrcX >= iSrcWidth)\n"
              "iSrcX = iSrcWidth - 1;\n"
  
          "if(iSrcY < 0)\n"
              "iSrcY = 0;\n"
          "else if(iSrcY >= iSrcHeight)\n"
              "iSrcY = iSrcHeight - 1;\n"
  
          "iSrcOffset = iSrcY*iSrcWidth + iSrcX;\n"
          "iBandValidLen = 1 + ((iSrcWidth*iSrcHeight)>>5);\n"
          "iSrcOffsetMask = (0x01 << (iSrcOffset & 0x1f));\n"
      "}\n"
  
      "if (useUnifiedSrcValid &&\n"
          "!((nUnifiedSrcValid[iSrcOffset>>5] & iSrcOffsetMask) ) )\n"
          "return FALSE;\n"
  
  "#ifdef USE_VEC\n"
      "if (!useUseBandSrcValid || !useBandSrcValid[bandNum] ||\n"
          "((nBandSrcValid[(iSrcOffset>>5)+iBandValidLen*bandNum    ] & iSrcOffsetMask)) )\n"
          "bHasValid = TRUE;\n"
  
      "if (!useUseBandSrcValid || !useBandSrcValid[bandNum+1] ||\n"
          "((nBandSrcValid[(iSrcOffset>>5)+iBandValidLen*(1+bandNum)] & iSrcOffsetMask)) )\n"
          "bHasValid = TRUE;\n"
  
      "if (!useUseBandSrcValid || !useBandSrcValid[bandNum+2] ||\n"
          "((nBandSrcValid[(iSrcOffset>>5)+iBandValidLen*(2+bandNum)] & iSrcOffsetMask)) )\n"
          "bHasValid = TRUE;\n"
  
      "if (!useUseBandSrcValid || !useBandSrcValid[bandNum+3] ||\n"
          "((nBandSrcValid[(iSrcOffset>>5)+iBandValidLen*(3+bandNum)] & iSrcOffsetMask)) )\n"
          "bHasValid = TRUE;\n"
  "#else\n"
      "if (!useUseBandSrcValid || !useBandSrcValid[bandNum] ||\n"
          "((nBandSrcValid[(iSrcOffset>>5)+iBandValidLen*bandNum    ] & iSrcOffsetMask)) )\n"
          "bHasValid = TRUE;\n"
  "#endif\n"
  
      "if (!bHasValid)\n"
          "return FALSE;\n"
  
      "const sampler_t samp =  CLK_NORMALIZED_COORDS_FALSE |\n"
                              "CLK_ADDRESS_CLAMP_TO_EDGE |\n"
                              "CLK_FILTER_NEAREST;\n"
  
  "#ifdef USE_VEC\n"
      "(*fReal) = read_imagef(srcReal, samp, iSrc);\n"
      "if (useImag)\n"
          "(*fImag) = read_imagef(srcImag, samp, iSrc);\n"
  "#else\n"
      "(*fReal) = read_imagef(srcReal, samp, iSrc).x;\n"
      "if (useImag)\n"
          "(*fImag) = read_imagef(srcImag, samp, iSrc).x;\n"
  "#endif\n"
  
      "if (useUnifiedSrcDensity) {\n"
          "(*fDensity) = fUnifiedSrcDensity[iSrcOffset];\n"
      "} else {\n"
          "(*fDensity) = 1.0f;\n"
          "return TRUE;\n"
      "}\n"
  
  "#ifdef USE_VEC\n"
      "return  (*fDensity).x > 0.0000001f || (*fDensity).y > 0.0000001f ||\n"
              "(*fDensity).z > 0.0000001f || (*fDensity).w > 0.0000001f;\n"
  "#else\n"
      "return (*fDensity) > 0.0000001f;\n"
  "#endif\n"
  "}\n"
  
  "int isValid(__global float *fUnifiedSrcDensity,\n"
              "__global int *nUnifiedSrcValid,\n"
              "float2 fSrcCoords )\n"
  "{\n"
      "if (fSrcCoords.x < 0.0f || fSrcCoords.y < 0.0f)\n"
          "return FALSE;\n"
  
      "int iSrcX = (int) (fSrcCoords.x - 0.5f);\n"
      "int iSrcY = (int) (fSrcCoords.y - 0.5f);\n"
  
      "if( iSrcX < 0 || iSrcX >= iSrcWidth || iSrcY < 0 || iSrcY >= iSrcHeight )\n"
          "return FALSE;\n"
  
      "int iSrcOffset = iSrcX + iSrcY * iSrcWidth;\n"
  
      "if (useUnifiedSrcDensity && fUnifiedSrcDensity[iSrcOffset] < 0.00001f)\n"
          "return FALSE;\n"
  
      "if (useUnifiedSrcValid &&\n"
          "!(nUnifiedSrcValid[iSrcOffset>>5] & (0x01 << (iSrcOffset & 0x1f))) )\n"
          "return FALSE;\n"
  
      "return TRUE;\n"
  "}\n"
  
  "float2 getSrcCoords(__read_only image2d_t srcCoords)\n"
  "{\n"
      // Find an appropriate place to sample the coordinates so we're still
      // accurate after linear interpolation.
      "int nDstX = get_global_id(0);\n"
      "int nDstY = get_global_id(1);\n"
      "float2  fDst = (float2)((0.5f * (float)iCoordMult + nDstX) /\n"
                                  "(float)((ceil((iDstWidth  - 1) / (float)iCoordMult) + 1) * iCoordMult), \n"
                              "(0.5f * (float)iCoordMult + nDstY) /\n"
                                  "(float)((ceil((iDstHeight - 1) / (float)iCoordMult) + 1) * iCoordMult));\n"
  
      // Check & return when the thread group overruns the image size
      "if (nDstX >= iDstWidth || nDstY >= iDstHeight)\n"
          "return (float2)(-99.0f, -99.0f);\n"
  
      "const sampler_t samp =  CLK_NORMALIZED_COORDS_TRUE |\n"
                              "CLK_ADDRESS_CLAMP_TO_EDGE |\n"
                              "CLK_FILTER_LINEAR;\n"
  
      "float4  fSrcCoords = read_imagef(srcCoords,samp,fDst);\n"
  
      "return (float2)(fSrcCoords.x, fSrcCoords.y);\n"
  "}\n";
  
      const char *kernBilinear =
  // ************************ Bilinear ************************
  "__kernel void resamp(__read_only image2d_t srcCoords,\n"
                      "__read_only image2d_t srcReal,\n"
                      "__read_only image2d_t srcImag,\n"
                      "__global float *fUnifiedSrcDensity,\n"
                      "__global int *nUnifiedSrcValid,\n"
                      "__constant char *useBandSrcValid,\n"
                      "__global int *nBandSrcValid,\n"
                      "__global outType *dstReal,\n"
                      "__global outType *dstImag,\n"
                      "__constant float *fDstNoDataReal,\n"
                      "__global float *dstDensity,\n"
                      "__global int *nDstValid,\n"
                      "const int bandNum)\n"
  "{\n"
      "float2  fSrc = getSrcCoords(srcCoords);\n"
      "if (!isValid(fUnifiedSrcDensity, nUnifiedSrcValid, fSrc))\n"
          "return;\n"
  
      "int     iSrcX = (int) floor(fSrc.x - 0.5f);\n"
      "int     iSrcY = (int) floor(fSrc.y - 0.5f);\n"
      "float   fRatioX = 1.5f - (fSrc.x - iSrcX);\n"
      "float   fRatioY = 1.5f - (fSrc.y - iSrcY);\n"
      "vecf    fReal, fImag, fDens;\n"
      "vecf    fAccumulatorReal = 0.0f, fAccumulatorImag = 0.0f;\n"
      "vecf    fAccumulatorDensity = 0.0f;\n"
      "float   fAccumulatorDivisor = 0.0f;\n"
  
      "if ( iSrcY >= 0 && iSrcY < iSrcHeight ) {\n"
          "float fMult1 = fRatioX * fRatioY;\n"
          "float fMult2 = (1.0f-fRatioX) * fRatioY;\n"
  
                  // Upper Left Pixel
                  "if ( iSrcX >= 0 && iSrcX < iSrcWidth\n"
                           "&& getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                                                  "useBandSrcValid, nBandSrcValid, (int2)(iSrcX, iSrcY),\n"
                                                  "bandNum, &fDens, &fReal, &fImag))\n"
                  "{\n"
                          "fAccumulatorDivisor += fMult1;\n"
                          "fAccumulatorReal += fReal * fMult1;\n"
                          "fAccumulatorImag += fImag * fMult1;\n"
                          "fAccumulatorDensity += fDens * fMult1;\n"
                  "}\n"
  
                  // Upper Right Pixel
                  "if ( iSrcX+1 >= 0 && iSrcX+1 < iSrcWidth\n"
                          "&& getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                                                  "useBandSrcValid, nBandSrcValid, (int2)(iSrcX+1, iSrcY),\n"
                                                  "bandNum, &fDens, &fReal, &fImag))\n"
                  "{\n"
                          "fAccumulatorDivisor += fMult2;\n"
                          "fAccumulatorReal += fReal * fMult2;\n"
                          "fAccumulatorImag += fImag * fMult2;\n"
                          "fAccumulatorDensity += fDens * fMult2;\n"
                  "}\n"
      "}\n"
  
      "if ( iSrcY+1 >= 0 && iSrcY+1 < iSrcHeight ) {\n"
          "float fMult1 = fRatioX * (1.0f-fRatioY);\n"
          "float fMult2 = (1.0f-fRatioX) * (1.0f-fRatioY);\n"
  
          // Lower Left Pixel
                  "if ( iSrcX >= 0 && iSrcX < iSrcWidth\n"
                          "&& getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                                                  "useBandSrcValid, nBandSrcValid, (int2)(iSrcX, iSrcY+1),\n"
                                                  "bandNum, &fDens, &fReal, &fImag))\n"
                  "{\n"
                          "fAccumulatorDivisor += fMult1;\n"
                          "fAccumulatorReal += fReal * fMult1;\n"
                          "fAccumulatorImag += fImag * fMult1;\n"
                          "fAccumulatorDensity += fDens * fMult1;\n"
                  "}\n"
  
                  // Lower Right Pixel
                  "if ( iSrcX+1 >= 0 && iSrcX+1 < iSrcWidth\n"
                          "&& getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                                                  "useBandSrcValid, nBandSrcValid, (int2)(iSrcX+1, iSrcY+1),\n"
                                                  "bandNum, &fDens, &fReal, &fImag))\n"
                  "{\n"
                          "fAccumulatorDivisor += fMult2;\n"
                          "fAccumulatorReal += fReal * fMult2;\n"
                          "fAccumulatorImag += fImag * fMult2;\n"
                          "fAccumulatorDensity += fDens * fMult2;\n"
                  "}\n"
      "}\n"
  
      // Compute and save final pixel
      "if ( fAccumulatorDivisor < 0.00001f ) {\n"
          "setPixel(dstReal, dstImag, dstDensity, nDstValid, fDstNoDataReal, bandNum,\n"
                  "0.0f, 0.0f, 0.0f );\n"
      "} else if ( fAccumulatorDivisor < 0.99999f || fAccumulatorDivisor > 1.00001f ) {\n"
          "setPixel(dstReal, dstImag, dstDensity, nDstValid, fDstNoDataReal, bandNum,\n"
                  "fAccumulatorDensity / fAccumulatorDivisor,\n"
                  "fAccumulatorReal / fAccumulatorDivisor,\n"
  "#if useImag != 0\n"
                  "fAccumulatorImag / fAccumulatorDivisor );\n"
  "#else\n"
                  "0.0f );\n"
  "#endif\n"
      "} else {\n"
          "setPixel(dstReal, dstImag, dstDensity, nDstValid, fDstNoDataReal, bandNum,\n"
                  "fAccumulatorDensity, fAccumulatorReal, fAccumulatorImag );\n"
      "}\n"
  "}\n";
  
      const char *kernCubic =
  // ************************ Cubic ************************
  "vecf cubicConvolution(float dist1, float dist2, float dist3,\n"
                          "vecf f0, vecf f1, vecf f2, vecf f3);\n"
  
  "vecf cubicConvolution(float dist1, float dist2, float dist3,\n"
                         "vecf f0, vecf f1, vecf f2, vecf f3)\n"
  "{\n"
     "return (  f1\n"
         "+ 0.5f * (dist1*(f2 - f0)\n"
                 "+ dist2*(2.0f*f0 - 5.0f*f1 + 4.0f*f2 - f3)\n"
                 "+ dist3*(3.0f*(f1 - f2) + f3 - f0)));\n"
  "}\n"
  
  // ************************ Cubic ************************
  "__kernel void resamp(__read_only image2d_t srcCoords,\n"
                       "__read_only image2d_t srcReal,\n"
                       "__read_only image2d_t srcImag,\n"
                       "__global float *fUnifiedSrcDensity,\n"
                       "__global int *nUnifiedSrcValid,\n"
                       "__constant char *useBandSrcValid,\n"
                       "__global int *nBandSrcValid,\n"
                       "__global outType *dstReal,\n"
                       "__global outType *dstImag,\n"
                       "__constant float *fDstNoDataReal,\n"
                       "__global float *dstDensity,\n"
                       "__global int *nDstValid,\n"
                       "const int bandNum)\n"
  "{\n"
      "int i;\n"
      "float2  fSrc = getSrcCoords(srcCoords);\n"
  
      "if (!isValid(fUnifiedSrcDensity, nUnifiedSrcValid, fSrc))\n"
          "return;\n"
  
      "int     iSrcX = (int) floor( fSrc.x - 0.5f );\n"
      "int     iSrcY = (int) floor( fSrc.y - 0.5f );\n"
      "float   fDeltaX = fSrc.x - 0.5f - (float)iSrcX;\n"
      "float   fDeltaY = fSrc.y - 0.5f - (float)iSrcY;\n"
      "float   fDeltaX2 = fDeltaX * fDeltaX;\n"
      "float   fDeltaY2 = fDeltaY * fDeltaY;\n"
      "float   fDeltaX3 = fDeltaX2 * fDeltaX;\n"
      "float   fDeltaY3 = fDeltaY2 * fDeltaY;\n"
      "vecf    afReal[4], afDens[4];\n"
  "#if useImag != 0\n"
      "vecf    afImag[4];\n"
  "#else\n"
      "vecf    fImag = 0.0f;\n"
  "#endif\n"
  
      // Loop over rows
      "for (i = -1; i < 3; ++i)\n"
      "{\n"
          "vecf    fReal1 = 0.0f, fReal2 = 0.0f, fReal3 = 0.0f, fReal4 = 0.0f;\n"
          "vecf    fDens1 = 0.0f, fDens2 = 0.0f, fDens3 = 0.0f, fDens4 = 0.0f;\n"
          "int hasPx;\n"
  "#if useImag != 0\n"
          "vecf    fImag1 = 0.0f, fImag2 = 0.0f, fImag3 = 0.0f, fImag4 = 0.0f;\n"
  
          //Get all the pixels for this row
          "hasPx  = getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                          "useBandSrcValid, nBandSrcValid, (int2)(iSrcX-1, iSrcY+i),\n"
                          "bandNum, &fDens1, &fReal1, &fImag1);\n"
  
          "hasPx |= getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                          "useBandSrcValid, nBandSrcValid, (int2)(iSrcX  , iSrcY+i),\n"
                          "bandNum, &fDens2, &fReal2, &fImag2);\n"
  
          "hasPx |= getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                          "useBandSrcValid, nBandSrcValid, (int2)(iSrcX+1, iSrcY+i),\n"
                          "bandNum, &fDens3, &fReal3, &fImag3);\n"
  
          "hasPx |= getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                          "useBandSrcValid, nBandSrcValid, (int2)(iSrcX+2, iSrcY+i),\n"
                          "bandNum, &fDens4, &fReal4, &fImag4);\n"
  "#else\n"
          //Get all the pixels for this row
          "hasPx  = getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                  "useBandSrcValid, nBandSrcValid, (int2)(iSrcX-1, iSrcY+i),\n"
                  "bandNum, &fDens1, &fReal1, &fImag);\n"
  
          "hasPx |= getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                  "useBandSrcValid, nBandSrcValid, (int2)(iSrcX  , iSrcY+i),\n"
                  "bandNum, &fDens2, &fReal2, &fImag);\n"
  
          "hasPx |= getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                  "useBandSrcValid, nBandSrcValid, (int2)(iSrcX+1, iSrcY+i),\n"
                  "bandNum, &fDens3, &fReal3, &fImag);\n"
  
          "hasPx |= getPixel(srcReal, srcImag, fUnifiedSrcDensity, nUnifiedSrcValid,\n"
                  "useBandSrcValid, nBandSrcValid, (int2)(iSrcX+2, iSrcY+i),\n"
                  "bandNum, &fDens4, &fReal4, &fImag);\n"
  "#endif\n"
  
          // Shortcut if no px
          "if (!hasPx) {\n"
              "afDens[i+1] = 0.0f;\n"
              "afReal[i+1] = 0.0f;\n"
  "#if useImag != 0\n"
              "afImag[i+1] = 0.0f;\n"
  "#endif\n"
              "continue;\n"
          "}\n"
  
          // Process this row
          "afDens[i+1] = cubicConvolution(fDeltaX, fDeltaX2, fDeltaX3, fDens1, fDens2, fDens3, fDens4);\n"
          "afReal[i+1] = cubicConvolution(fDeltaX, fDeltaX2, fDeltaX3, fReal1, fReal2, fReal3, fReal4);\n"
  "#if useImag != 0\n"
          "afImag[i+1] = cubicConvolution(fDeltaX, fDeltaX2, fDeltaX3, fImag1, fImag2, fImag3, fImag4);\n"
  "#endif\n"
      "}\n"
  
      // Compute and save final pixel
      "setPixel(dstReal, dstImag, dstDensity, nDstValid, fDstNoDataReal, bandNum,\n"
               "cubicConvolution(fDeltaY, fDeltaY2, fDeltaY3, afDens[0], afDens[1], afDens[2], afDens[3]),\n"
               "cubicConvolution(fDeltaY, fDeltaY2, fDeltaY3, afReal[0], afReal[1], afReal[2], afReal[3]),\n"
  "#if useImag != 0\n"
               "cubicConvolution(fDeltaY, fDeltaY2, fDeltaY3, afImag[0], afImag[1], afImag[2], afImag[3]) );\n"
  "#else\n"
               "fImag );\n"
  "#endif\n"
  "}\n";
  
      const char *kernResampler =
  // ************************ LanczosSinc ************************
  
  "float lanczosSinc( float fX, float fR );\n"
  "float bSpline( float x );\n"
  
  "float lanczosSinc( float fX, float fR )\n"
  "{\n"
      "if ( fX > fR || fX < -fR)\n"
          "return 0.0f;\n"
      "if ( fX == 0.0f )\n"
          "return 1.0f;\n"
  
      "float fPIX = PI * fX;\n"
      "return ( sin(fPIX) / fPIX ) * ( sin(fPIX / fR) * fR / fPIX );\n"
  "}\n"
  
  // ************************ Bicubic Spline ************************
  
  "float bSpline( float x )\n"
  "{\n"
      "float xp2 = x + 2.0f;\n"
      "float xp1 = x + 1.0f;\n"
      "float xm1 = x - 1.0f;\n"
      "float xp2c = xp2 * xp2 * xp2;\n"
  
      "return (((xp2 > 0.0f)?((xp1 > 0.0f)?((x > 0.0f)?((xm1 > 0.0f)?\n"
                                                       "-4.0f * xm1*xm1*xm1:0.0f) +\n"
                                           "6.0f * x*x*x:0.0f) +\n"
                             "-4.0f * xp1*xp1*xp1:0.0f) +\n"
               "xp2c:0.0f) ) * 0.166666666666666666666f;\n"
  "}\n"
  
  // ************************ General Resampler ************************
  
  "__kernel void resamp(__read_only image2d_t srcCoords,\n"
                       "__read_only image2d_t srcReal,\n"
                       "__read_only image2d_t srcImag,\n"
                       "__global float *fUnifiedSrcDensity,\n"
                       "__global int *nUnifiedSrcValid,\n"
                       "__constant char *useBandSrcValid,\n"
                       "__global int *nBandSrcValid,\n"
                       "__global outType *dstReal,\n"
                       "__global outType *dstImag,\n"
                       "__constant float *fDstNoDataReal,\n"
                       "__global float *dstDensity,\n"
                       "__global int *nDstValid,\n"
                       "const int bandNum)\n"
  "{\n"
      "float2  fSrc = getSrcCoords(srcCoords);\n"
  
      "if (!isValid(fUnifiedSrcDensity, nUnifiedSrcValid, fSrc))\n"
          "return;\n"
  
      "int     iSrcX = floor( fSrc.x - 0.5f );\n"
      "int     iSrcY = floor( fSrc.y - 0.5f );\n"
      "float   fDeltaX = fSrc.x - 0.5f - (float)iSrcX;\n"
      "float   fDeltaY = fSrc.y - 0.5f - (float)iSrcY;\n"
  
      "vecf  fAccumulatorReal = 0.0f, fAccumulatorImag = 0.0f;\n"
      "vecf  fAccumulatorDensity = 0.0f;\n"
      "float fAccumulatorWeight = 0.0f;\n"
      "int   i, j;\n"
  
      // Loop over pixel rows in the kernel
      "for ( j = nFiltInitY; j <= nYRadius; ++j )\n"
      "{\n"
          "float   fWeight1;\n"
          "int2 iSrc = (int2)(0, iSrcY + j);\n"
  
          // Skip sampling over edge of image
          "if ( iSrc.y < 0 || iSrc.y >= iSrcHeight )\n"
              "continue;\n"
  
          // Select the resampling algorithm
          "if ( doCubicSpline )\n"
              // Calculate the Y weight
              "fWeight1 = ( fYScale < 1.0f ) ?\n"
                  "bSpline(((float)j) * fYScale) * fYScale :\n"
                  "bSpline(((float)j) - fDeltaY);\n"
          "else\n"
              "fWeight1 = ( fYScale < 1.0f ) ?\n"
                  "lanczosSinc(j * fYScale, fYFilter) * fYScale :\n"
                  "lanczosSinc(j - fDeltaY, fYFilter);\n"
  
          // Iterate over pixels in row
          "for ( i = nFiltInitX; i <= nXRadius; ++i )\n"
          "{\n"
              "float fWeight2;\n"
              "vecf fDensity = 0.0f, fReal, fImag;\n"
              "iSrc.x = iSrcX + i;\n"
  
              // Skip sampling at edge of image
              // Skip sampling when invalid pixel
              "if ( iSrc.x < 0 || iSrc.x >= iSrcWidth || \n"
                    "!getPixel(srcReal, srcImag, fUnifiedSrcDensity,\n"
                              "nUnifiedSrcValid, useBandSrcValid, nBandSrcValid,\n"
                              "iSrc, bandNum, &fDensity, &fReal, &fImag) )\n"
                  "continue;\n"
  
              // Choose among possible algorithms
              "if ( doCubicSpline )\n"
                  // Calculate & save the X weight
                  "fWeight2 = fWeight1 * ((fXScale < 1.0f ) ?\n"
                      "bSpline((float)i * fXScale) * fXScale :\n"
                      "bSpline(fDeltaX - (float)i));\n"
              "else\n"
                  // Calculate & save the X weight
                  "fWeight2 = fWeight1 * ((fXScale < 1.0f ) ?\n"
                      "lanczosSinc(i * fXScale, fXFilter) * fXScale :\n"
                      "lanczosSinc(i - fDeltaX, fXFilter));\n"
  
              // Accumulate!
              "fAccumulatorReal += fReal * fWeight2;\n"
              "fAccumulatorImag += fImag * fWeight2;\n"
              "fAccumulatorDensity += fDensity * fWeight2;\n"
              "fAccumulatorWeight += fWeight2;\n"
          "}\n"
      "}\n"
  
      "if ( fAccumulatorWeight < 0.000001f ) {\n"
          "setPixel(dstReal, dstImag, dstDensity, nDstValid, fDstNoDataReal, bandNum,\n"
                   "0.0f, 0.0f, 0.0f);\n"
      "} else if ( fAccumulatorWeight < 0.99999f || fAccumulatorWeight > 1.00001f ) {\n"
          // Calculate the output taking into account weighting
          "setPixel(dstReal, dstImag, dstDensity, nDstValid, fDstNoDataReal, bandNum,\n"
                   "fAccumulatorDensity / fAccumulatorWeight,\n"
                   "fAccumulatorReal / fAccumulatorWeight,\n"
  "#if useImag != 0\n"
                   "fAccumulatorImag / fAccumulatorWeight );\n"
  "#else\n"
                   "0.0f );\n"
  "#endif\n"
      "} else {\n"
          "setPixel(dstReal, dstImag, dstDensity, nDstValid, fDstNoDataReal, bandNum,\n"
                   "fAccumulatorDensity, fAccumulatorReal, fAccumulatorImag);\n"
      "}\n"
  "}\n";
  
      //Defines based on image format
      switch (warper->imageFormat) {
          case CL_FLOAT:
              dstMinVal = -FLT_MAX;
              dstMaxVal = FLT_MAX;
              outType = "float";
              break;
          case CL_SNORM_INT8:
              dstMinVal = -128.0;
              dstMaxVal = 127.0;
              outType = "char";
              break;
          case CL_UNORM_INT8:
              dstMinVal = 0.0;
              dstMaxVal = 255.0;
              outType = "uchar";
              break;
          case CL_SNORM_INT16:
              dstMinVal = -32768.0;
              dstMaxVal = 32767.0;
              outType = "short";
              break;
          case CL_UNORM_INT16:
              dstMinVal = 0.0;
              dstMaxVal = 65535.0;
              outType = "ushort";
              break;
      }
  
      //Use vector format?
      if (useVec) {
          dUseVec = "-D USE_VEC";
          dVecf = "float4";
      }
  
      //Assemble the kernel from parts. The compiler is unable to handle multiple
      //kernels in one string with more than a few __constant modifiers each.
      if (warper->resampAlg == OCL_Bilinear)
          snprintf(progBuf, PROGBUF_SIZE, "%s\n%s", kernGenFuncs, kernBilinear);
      else if (warper->resampAlg == OCL_Cubic)
          snprintf(progBuf, PROGBUF_SIZE, "%s\n%s", kernGenFuncs, kernCubic);
      else
          snprintf(progBuf, PROGBUF_SIZE, "%s\n%s", kernGenFuncs, kernResampler);
  
      //Actually make the program from assembled source
      program = clCreateProgramWithSource(warper->context, 1, (const char**)&progBuf,
                                          NULL, &err);
      handleErrGoto(err, error_final);
  
      //Assemble the compiler arg string for speed. All invariants should be defined here.
      snprintf(buffer, PROGBUF_SIZE,
               "-cl-fast-relaxed-math -Werror -D FALSE=0 -D TRUE=1 "
              "%s"
              "-D iSrcWidth=%d -D iSrcHeight=%d -D iDstWidth=%d -D iDstHeight=%d "
              "-D useUnifiedSrcDensity=%d -D useUnifiedSrcValid=%d "
              "-D useDstDensity=%d -D useDstValid=%d -D useImag=%d "
              "-D fXScale=%015.15lff -D fYScale=%015.15lff -D fXFilter=%015.15lff -D fYFilter=%015.15lff "
              "-D nXRadius=%d -D nYRadius=%d -D nFiltInitX=%d -D nFiltInitY=%d "
              "-D PI=%015.15lff -D outType=%s -D dstMinVal=%015.15lff -D dstMaxVal=%015.15lff "
              "-D useDstNoDataReal=%d -D vecf=%s %s -D doCubicSpline=%d "
              "-D useUseBandSrcValid=%d -D iCoordMult=%d ",
              (warper->imageFormat == CL_FLOAT) ? "-D USE_CLAMP_TO_DST_FLOAT=1 " : "",
              warper->srcWidth, warper->srcHeight, warper->dstWidth, warper->dstHeight,
              warper->useUnifiedSrcDensity, warper->useUnifiedSrcValid,
              warper->useDstDensity, warper->useDstValid, warper->imagWorkCL != NULL,
              dfXScale, dfYScale, dfXFilter, dfYFilter,
              nXRadius, nYRadius, nFiltInitX, nFiltInitY,
              M_PI, outType, dstMinVal, dstMaxVal, warper->fDstNoDataRealCL != NULL,
              dVecf, dUseVec, warper->resampAlg == OCL_CubicSpline,
              warper->nBandSrcValidCL != NULL, warper->coordMult);
  
      (*clErr) = err = clBuildProgram(program, 1, &(warper->dev), buffer, NULL, NULL);
  
      //Detailed debugging info
      if (err != CL_SUCCESS)
      {
          const char* pszStatus = "unknown_status";
          err = clGetProgramBuildInfo(program, warper->dev, CL_PROGRAM_BUILD_LOG,
                                      128000*sizeof(char), buffer, NULL);
          handleErrGoto(err, error_free_program);
  
          CPLError(CE_Failure, CPLE_AppDefined, "Error: Failed to build program executable!\nBuild Log:\n%s", buffer);
  
          err = clGetProgramBuildInfo(program, warper->dev, CL_PROGRAM_BUILD_STATUS,
                                      128000*sizeof(char), buffer, NULL);
          handleErrGoto(err, error_free_program);
  
          if(buffer[0] == CL_BUILD_NONE)
              pszStatus = "CL_BUILD_NONE";
          else if(buffer[0] == CL_BUILD_ERROR)
              pszStatus = "CL_BUILD_ERROR";
          else if(buffer[0] == CL_BUILD_SUCCESS)
              pszStatus = "CL_BUILD_SUCCESS";
          else if(buffer[0] == CL_BUILD_IN_PROGRESS)
              pszStatus = "CL_BUILD_IN_PROGRESS";
  
          CPLDebug("OpenCL", "Build Status: %s\nProgram Source:\n%s", pszStatus, progBuf);
          goto error_free_program;
      }
  
      kernel = clCreateKernel(program, "resamp", &err);
      handleErrGoto(err, error_free_program);
  
      err = clReleaseProgram(program);
      handleErrGoto(err, error_final);
  
      CPLFree(buffer);
      CPLFree(progBuf);
      return kernel;
  
  error_free_program:
      err = clReleaseProgram(program);
  
  error_final:
      CPLFree(buffer);
      CPLFree(progBuf);
      return NULL;
  }
  
  /*
   Alloc & copy the coordinate data from host working memory to the device. The
   working memory should be a pinned, linear, array of floats. This allows us to
   allocate and copy all data in one step. The pointer to the device memory is
   saved and set as the appropriate argument number.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int set_coord_data (struct oclWarper *warper, cl_mem *xy)
  {
      cl_int err = CL_SUCCESS;
      cl_image_format imgFmt;
  
      //Copy coord data to the device
      imgFmt.image_channel_order = warper->xyChOrder;
      imgFmt.image_channel_data_type = CL_FLOAT;
  
  #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6) || defined(__clang__)
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
  #endif
      (*xy) = clCreateImage2D(warper->context,
                              CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, &imgFmt,
                              (size_t) warper->xyWidth,
                              (size_t) warper->xyHeight,
                              (size_t) sizeof(float) * warper->xyChSize * warper->xyWidth,
                              warper->xyWork, &err);
  #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6) || defined(__clang__)
  #pragma GCC diagnostic pop
  #endif
      handleErr(err);
  
      //Free the source memory, now that it's copied we don't need it
      freeCLMem(warper->xyWorkCL, warper->xyWork);
  
      //Set up argument
      if (warper->kern1 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern1, 0, sizeof(cl_mem), xy));
      }
      if (warper->kern4 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern4, 0, sizeof(cl_mem), xy));
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Sets the unified density & valid data structures. These are optional structures
   from GDAL, and as such if they are NULL a small placeholder memory segment is
   defined. This is because the spec is unclear on if a NULL value can be passed
   as a kernel argument in place of memory. If it's not NULL, the data is copied
   from the working memory to the device memory. After that, we check if we are
   using the per-band validity mask, and set that as appropriate. At the end, the
   CL mem is passed as the kernel arguments.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int set_unified_data(struct oclWarper *warper,
                          cl_mem *unifiedSrcDensityCL, cl_mem *unifiedSrcValidCL,
                          float *unifiedSrcDensity, unsigned int *unifiedSrcValid,
                          cl_mem *useBandSrcValidCL, cl_mem *nBandSrcValidCL)
  {
      cl_int err = CL_SUCCESS;
      size_t sz = warper->srcWidth * warper->srcHeight;
      int useValid = warper->nBandSrcValidCL != NULL;
      //32 bits in the mask
      int validSz = (int)(sizeof(int) * ((31 + sz) >> 5));
  
      //Copy unifiedSrcDensity if it exists
      if (unifiedSrcDensity == NULL) {
          //Alloc dummy device RAM
          (*unifiedSrcDensityCL) = clCreateBuffer(warper->context, CL_MEM_READ_ONLY, 1, NULL, &err);
          handleErr(err);
      } else {
          //Alloc & copy all density data
          (*unifiedSrcDensityCL) = clCreateBuffer(warper->context,
                                                  CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
                                                  sizeof(float) * sz, unifiedSrcDensity, &err);
          handleErr(err);
      }
  
      //Copy unifiedSrcValid if it exists
      if (unifiedSrcValid == NULL) {
          //Alloc dummy device RAM
          (*unifiedSrcValidCL) = clCreateBuffer(warper->context, CL_MEM_READ_ONLY, 1, NULL, &err);
          handleErr(err);
      } else {
          //Alloc & copy all validity data
          (*unifiedSrcValidCL) = clCreateBuffer(warper->context,
                                                CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
                                                validSz, unifiedSrcValid, &err);
          handleErr(err);
      }
  
      // Set the band validity usage
      if(useValid) {
          (*useBandSrcValidCL) = clCreateBuffer(warper->context,
                                                CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
                                                sizeof(char) * warper->numBands,
                                                warper->useBandSrcValid, &err);
          handleErr(err);
      } else {
          //Make a fake image so we don't have a NULL pointer
          (*useBandSrcValidCL) = clCreateBuffer(warper->context, CL_MEM_READ_ONLY, 1, NULL, &err);
          handleErr(err);
      }
  
      //Do a more thorough check for validity
      if (useValid) {
          int i;
          useValid = FALSE;
          for (i = 0; i < warper->numBands; ++i)
              if (warper->useBandSrcValid[i])
                  useValid = TRUE;
      }
  
      //And the validity mask if needed
      if (useValid) {
          (*nBandSrcValidCL) = clCreateBuffer(warper->context,
                                              CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
                                              warper->numBands * validSz,
                                              warper->nBandSrcValid, &err);
          handleErr(err);
      } else {
          //Make a fake image so we don't have a NULL pointer
          (*nBandSrcValidCL) = clCreateBuffer(warper->context, CL_MEM_READ_ONLY, 1, NULL, &err);
          handleErr(err);
      }
  
      //Set up arguments
      if (warper->kern1 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern1, 3, sizeof(cl_mem), unifiedSrcDensityCL));
          handleErr(err = clSetKernelArg(warper->kern1, 4, sizeof(cl_mem), unifiedSrcValidCL));
          handleErr(err = clSetKernelArg(warper->kern1, 5, sizeof(cl_mem), useBandSrcValidCL));
          handleErr(err = clSetKernelArg(warper->kern1, 6, sizeof(cl_mem), nBandSrcValidCL));
      }
      if (warper->kern4 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern4, 3, sizeof(cl_mem), unifiedSrcDensityCL));
          handleErr(err = clSetKernelArg(warper->kern4, 4, sizeof(cl_mem), unifiedSrcValidCL));
          handleErr(err = clSetKernelArg(warper->kern4, 5, sizeof(cl_mem), useBandSrcValidCL));
          handleErr(err = clSetKernelArg(warper->kern4, 6, sizeof(cl_mem), nBandSrcValidCL));
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Here we set the per-band raster data. First priority is the real raster data,
   of course. Then, if applicable, we set the additional image channel. Once this
   data is copied to the device, it can be freed on the host, so that is done
   here. Finally the appropriate kernel arguments are set.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int set_src_rast_data (struct oclWarper *warper, int iNum, size_t sz,
                            cl_channel_order chOrder, cl_mem *srcReal, cl_mem *srcImag)
  {
      cl_image_format imgFmt;
      cl_int err = CL_SUCCESS;
      int useImagWork = warper->imagWork.v != NULL && warper->imagWork.v[iNum] != NULL;
  
      //Set up image vars
      imgFmt.image_channel_order = chOrder;
      imgFmt.image_channel_data_type = warper->imageFormat;
  
      //Create & copy the source image
  #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6) || defined(__clang__)
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
  #endif
  
      (*srcReal) = clCreateImage2D(warper->context,
                                   CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, &imgFmt,
                                   (size_t) warper->srcWidth, (size_t) warper->srcHeight,
                                   sz * warper->srcWidth, warper->realWork.v[iNum], &err);
      handleErr(err);
  
      //And the source image parts if needed
      if (useImagWork) {
          (*srcImag) = clCreateImage2D(warper->context,
                                       CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, &imgFmt,
                                       (size_t) warper->srcWidth, (size_t) warper->srcHeight,
                                       sz * warper->srcWidth, warper->imagWork.v[iNum], &err);
          handleErr(err);
      } else {
          //Make a fake image so we don't have a NULL pointer
  
          char dummyImageData[16];
          (*srcImag) = clCreateImage2D(warper->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, &imgFmt,
                                      1, 1, sz, dummyImageData, &err);
  
          handleErr(err);
      }
  #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6) || defined(__clang__)
  #pragma GCC diagnostic pop
  #endif
  
      //Free the source memory, now that it's copied we don't need it
      freeCLMem(warper->realWorkCL[iNum], warper->realWork.v[iNum]);
      if (warper->imagWork.v != NULL) {
          freeCLMem(warper->imagWorkCL[iNum], warper->imagWork.v[iNum]);
      }
  
      //Set up per-band arguments
      if (warper->kern1 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern1, 1, sizeof(cl_mem), srcReal));
          handleErr(err = clSetKernelArg(warper->kern1, 2, sizeof(cl_mem), srcImag));
      }
      if (warper->kern4 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern4, 1, sizeof(cl_mem), srcReal));
          handleErr(err = clSetKernelArg(warper->kern4, 2, sizeof(cl_mem), srcImag));
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Set the destination data for the raster. Although it's the output, it still
   is copied to the device because some blending is done there. First the real
   data is allocated and copied, then the imag data is allocated and copied if
   needed. They are then set as the appropriate arguments to the kernel.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int set_dst_rast_data(struct oclWarper *warper, int iImg, size_t sz,
                           cl_mem *dstReal, cl_mem *dstImag)
  {
      cl_int err = CL_SUCCESS;
      sz *= warper->dstWidth * warper->dstHeight;
  
      //Copy the dst real data
      (*dstReal) = clCreateBuffer(warper->context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
                                  sz, warper->dstRealWork.v[iImg], &err);
      handleErr(err);
  
      //Copy the dst imag data if exists
      if (warper->dstImagWork.v != NULL && warper->dstImagWork.v[iImg] != NULL) {
          (*dstImag) = clCreateBuffer(warper->context,
                                      CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
                                      sz, warper->dstImagWork.v[iImg], &err);
          handleErr(err);
      } else {
          (*dstImag) = clCreateBuffer(warper->context, CL_MEM_READ_WRITE, 1, NULL, &err);
          handleErr(err);
      }
  
      //Set up per-band arguments
      if (warper->kern1 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern1, 7, sizeof(cl_mem), dstReal));
          handleErr(err = clSetKernelArg(warper->kern1, 8, sizeof(cl_mem), dstImag));
      }
      if (warper->kern4 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern4, 7, sizeof(cl_mem), dstReal));
          handleErr(err = clSetKernelArg(warper->kern4, 8, sizeof(cl_mem), dstImag));
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Read the final raster data back from the graphics card to working memory. This
   copies both the real memory and the imag memory if appropriate.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int get_dst_rast_data(struct oclWarper *warper, int iImg, size_t wordSz,
                           cl_mem dstReal, cl_mem dstImag)
  {
      cl_int err = CL_SUCCESS;
      size_t sz = warper->dstWidth * warper->dstHeight * wordSz;
  
      //Copy from dev into working memory
      handleErr(err = clEnqueueReadBuffer(warper->queue, dstReal,
                                          CL_FALSE, 0, sz, warper->dstRealWork.v[iImg],
                                          0, NULL, NULL));
  
      //If we are expecting the imag channel, then copy it back also
      if (warper->dstImagWork.v != NULL && warper->dstImagWork.v[iImg] != NULL) {
          handleErr(err = clEnqueueReadBuffer(warper->queue, dstImag,
                                              CL_FALSE, 0, sz, warper->dstImagWork.v[iImg],
                                              0, NULL, NULL));
      }
  
      //The copy requests were non-blocking, so we'll need to make sure they finish.
      handleErr(err = clFinish(warper->queue));
  
      return CL_SUCCESS;
  }
  
  /*
   Set the destination image density & validity mask on the device. This is used
   to blend the final output image with the existing buffer. This handles the
   unified structures that apply to all bands. After the buffers are created and
   copied, they are set as kernel arguments.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int set_dst_data(struct oclWarper *warper,
                      cl_mem *dstDensityCL, cl_mem *dstValidCL, cl_mem *dstNoDataRealCL,
                      float *dstDensity, unsigned int *dstValid, float *dstNoDataReal)
  {
      cl_int err = CL_SUCCESS;
      size_t sz = warper->dstWidth * warper->dstHeight;
  
      //Copy the no-data value(s)
      if (dstNoDataReal == NULL) {
          (*dstNoDataRealCL) = clCreateBuffer(warper->context, CL_MEM_READ_ONLY, 1, NULL, &err);
          handleErr(err);
      } else {
          (*dstNoDataRealCL) = clCreateBuffer(warper->context,
                                           CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
                                           sizeof(float) * warper->numBands, dstNoDataReal, &err);
          handleErr(err);
      }
  
      //Copy unifiedSrcDensity if it exists
      if (dstDensity == NULL) {
          (*dstDensityCL) = clCreateBuffer(warper->context, CL_MEM_READ_ONLY, 1, NULL, &err);
          handleErr(err);
      } else {
          (*dstDensityCL) = clCreateBuffer(warper->context,
                                           CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
                                           sizeof(float) * sz, dstDensity, &err);
          handleErr(err);
      }
  
      //Copy unifiedSrcValid if it exists
      if (dstValid == NULL) {
          (*dstValidCL) = clCreateBuffer(warper->context, CL_MEM_READ_ONLY, 1, NULL, &err);
          handleErr(err);
      } else {
          (*dstValidCL) = clCreateBuffer(warper->context,
                                         CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
                                         sizeof(int) * ((31 + sz) >> 5), dstValid, &err);
          handleErr(err);
      }
  
      //Set up arguments
      if (warper->kern1 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern1,  9, sizeof(cl_mem), dstNoDataRealCL));
          handleErr(err = clSetKernelArg(warper->kern1, 10, sizeof(cl_mem), dstDensityCL));
          handleErr(err = clSetKernelArg(warper->kern1, 11, sizeof(cl_mem), dstValidCL));
      }
      if (warper->kern4 != NULL) {
          handleErr(err = clSetKernelArg(warper->kern4,  9, sizeof(cl_mem), dstNoDataRealCL));
          handleErr(err = clSetKernelArg(warper->kern4, 10, sizeof(cl_mem), dstDensityCL));
          handleErr(err = clSetKernelArg(warper->kern4, 11, sizeof(cl_mem), dstValidCL));
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Go ahead and execute the kernel. This handles some housekeeping stuff like the
   run dimensions. When running in debug mode, it times the kernel call and prints
   the execution time.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int execute_kern(struct oclWarper *warper, cl_kernel kern, size_t loc_size)
  {
      cl_int err = CL_SUCCESS;
      cl_event ev;
      size_t ceil_runs[2];
      size_t group_size[2];
  #ifdef DEBUG_OPENCL
      size_t start_time = 0;
      size_t end_time;
      char *vecTxt = "";
  #endif
  
      // Use a likely X-dimension which is a power of 2
      if (loc_size >= 512)
          group_size[0] = 32;
      else if (loc_size >= 64)
          group_size[0] = 16;
      else if (loc_size > 8)
          group_size[0] = 8;
      else
          group_size[0] = 1;
  
      if (group_size[0] > loc_size)
          group_size[1] = group_size[0]/loc_size;
      else
          group_size[1] = 1;
  
      //Round up num_runs to find the dim of the block of pixels we'll be processing
      if(warper->dstWidth % group_size[0])
          ceil_runs[0] = warper->dstWidth + group_size[0] - warper->dstWidth % group_size[0];
      else
          ceil_runs[0] = warper->dstWidth;
  
      if(warper->dstHeight % group_size[1])
          ceil_runs[1] = warper->dstHeight + group_size[1] - warper->dstHeight % group_size[1];
      else
          ceil_runs[1] = warper->dstHeight;
  
  #ifdef DEBUG_OPENCL
      handleErr(err = clSetCommandQueueProperty(warper->queue, CL_QUEUE_PROFILING_ENABLE, CL_TRUE, NULL));
  #endif
  
      // Run the calculation by enqueuing it and forcing the
      // command queue to complete the task
      handleErr(err = clEnqueueNDRangeKernel(warper->queue, kern, 2, NULL,
                                             ceil_runs, group_size, 0, NULL, &ev));
      handleErr(err = clFinish(warper->queue));
  
  #ifdef DEBUG_OPENCL
      handleErr(err = clGetEventProfilingInfo(ev, CL_PROFILING_COMMAND_START,
                                              sizeof(size_t), &start_time, NULL));
      handleErr(err = clGetEventProfilingInfo(ev, CL_PROFILING_COMMAND_END,
                                              sizeof(size_t), &end_time, NULL));
      assert(end_time != 0);
      assert(start_time != 0);
      handleErr(err = clReleaseEvent(ev));
      if (kern == warper->kern4)
          vecTxt = "(vec)";
  
      CPLDebug("OpenCL", "Kernel Time: %6s %10lu", vecTxt, (long int)((end_time-start_time)/100000));
  #endif
      return CL_SUCCESS;
  }
  
  /*
   Copy data from a raw source to the warper's working memory. If the imag
   channel is expected, then the data will be de-interlaced into component blocks
   of memory.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  static
  cl_int set_img_data(struct oclWarper *warper, void *srcImgData,
                      unsigned int width, unsigned int height, int isSrc,
                      unsigned int bandNum, void **dstRealImgs, void **dstImagImgs)
  {
      unsigned int imgChSize = warper->imgChSize1;
      unsigned int iSrcY, i;
      unsigned int vecOff = 0;
      unsigned int imgNum = bandNum;
      void *dstReal = NULL;
      void *dstImag = NULL;
  
      // Handle vector if needed
      if (warper->useVec && (int)bandNum < warper->numBands - warper->numBands % 4) {
          imgChSize = warper->imgChSize4;
          vecOff = bandNum % 4;
          imgNum = bandNum / 4;
      } else if(warper->useVec) {
          imgNum = bandNum / 4 + bandNum % 4;
      }
  
      // Set the images as needed
      dstReal = dstRealImgs[imgNum];
      if(dstImagImgs == NULL)
          dstImag = NULL;
      else
          dstImag = dstImagImgs[imgNum];
  
      // Set stuff for dst imgs
      if (!isSrc) {
          vecOff *= height * width;
          imgChSize = 1;
      }
  
      // Copy values as needed
      if (warper->imagWorkCL == NULL && !(warper->useVec && isSrc)) {
          //Set memory size & location depending on the data type
          //This is the ideal code path for speed
          switch (warper->imageFormat) {
              case CL_UNORM_INT8:
              {
                  unsigned char *realDst = &(((unsigned char *)dstReal)[vecOff]);
                  memcpy(realDst, srcImgData, width*height*sizeof(unsigned char));
                  break;
              }
              case CL_SNORM_INT8:
              {
                  char *realDst = &(((char *)dstReal)[vecOff]);
                  memcpy(realDst, srcImgData, width*height*sizeof(char));
                  break;
              }
              case CL_UNORM_INT16:
              {
                  unsigned short *realDst = &(((unsigned short *)dstReal)[vecOff]);
                  memcpy(realDst, srcImgData, width*height*sizeof(unsigned short));
                  break;
              }
              case CL_SNORM_INT16:
              {
                  short *realDst = &(((short *)dstReal)[vecOff]);
                  memcpy(realDst, srcImgData, width*height*sizeof(short));
                  break;
              }
              case CL_FLOAT:
              {
                  float *realDst = &(((float *)dstReal)[vecOff]);
                  memcpy(realDst, srcImgData, width*height*sizeof(float));
                  break;
              }
          }
      } else if (warper->imagWorkCL == NULL) {
          //We need to space the values due to OpenCL implementation reasons
          for( iSrcY = 0; iSrcY < height; iSrcY++ )
          {
              int pxOff = width*iSrcY;
              int imgOff = imgChSize*pxOff + vecOff;
              //Copy & deinterleave interleaved data
              switch (warper->imageFormat) {
                  case CL_UNORM_INT8:
                  {
                      unsigned char *realDst = &(((unsigned char *)dstReal)[imgOff]);
                      unsigned char *dataSrc = &(((unsigned char *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i)
                          realDst[imgChSize*i] = dataSrc[i];
                  }
                      break;
                  case CL_SNORM_INT8:
                  {
                      char *realDst = &(((char *)dstReal)[imgOff]);
                      char *dataSrc = &(((char *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i)
                          realDst[imgChSize*i] = dataSrc[i];
                  }
                      break;
                  case CL_UNORM_INT16:
                  {
                      unsigned short *realDst = &(((unsigned short *)dstReal)[imgOff]);
                      unsigned short *dataSrc = &(((unsigned short *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i)
                          realDst[imgChSize*i] = dataSrc[i];
                  }
                      break;
                  case CL_SNORM_INT16:
                  {
                      short *realDst = &(((short *)dstReal)[imgOff]);
                      short *dataSrc = &(((short *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i)
                          realDst[imgChSize*i] = dataSrc[i];
                  }
                      break;
                  case CL_FLOAT:
                  {
                      float *realDst = &(((float *)dstReal)[imgOff]);
                      float *dataSrc = &(((float *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i)
                          realDst[imgChSize*i] = dataSrc[i];
                  }
                      break;
              }
          }
      } else {
          //Copy, deinterleave, & space interleaved data
          for( iSrcY = 0; iSrcY < height; iSrcY++ )
          {
              int pxOff = width*iSrcY;
              int imgOff = imgChSize*pxOff + vecOff;
              //Copy & deinterleave interleaved data
              switch (warper->imageFormat) {
                  case CL_FLOAT:
                  {
                      float *realDst = &(((float *)dstReal)[imgOff]);
                      float *imagDst = &(((float *)dstImag)[imgOff]);
                      float *dataSrc = &(((float *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i) {
                          realDst[imgChSize*i] = dataSrc[i*2  ];
                          imagDst[imgChSize*i] = dataSrc[i*2+1];
                      }
                  }
                      break;
                  case CL_SNORM_INT8:
                  {
                      char *realDst = &(((char *)dstReal)[imgOff]);
                      char *imagDst = &(((char *)dstImag)[imgOff]);
                      char *dataSrc = &(((char *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i) {
                          realDst[imgChSize*i] = dataSrc[i*2  ];
                          imagDst[imgChSize*i] = dataSrc[i*2+1];
                      }
                  }
                      break;
                  case CL_UNORM_INT8:
                  {
                      unsigned char *realDst = &(((unsigned char *)dstReal)[imgOff]);
                      unsigned char *imagDst = &(((unsigned char *)dstImag)[imgOff]);
                      unsigned char *dataSrc = &(((unsigned char *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i) {
                          realDst[imgChSize*i] = dataSrc[i*2  ];
                          imagDst[imgChSize*i] = dataSrc[i*2+1];
                      }
                  }
                      break;
                  case CL_SNORM_INT16:
                  {
                      short *realDst = &(((short *)dstReal)[imgOff]);
                      short *imagDst = &(((short *)dstImag)[imgOff]);
                      short *dataSrc = &(((short *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i) {
                          realDst[imgChSize*i] = dataSrc[i*2  ];
                          imagDst[imgChSize*i] = dataSrc[i*2+1];
                      }
                  }
                      break;
                  case CL_UNORM_INT16:
                  {
                      unsigned short *realDst = &(((unsigned short *)dstReal)[imgOff]);
                      unsigned short *imagDst = &(((unsigned short *)dstImag)[imgOff]);
                      unsigned short *dataSrc = &(((unsigned short *)srcImgData)[pxOff]);
                      for (i = 0; i < width; ++i) {
                          realDst[imgChSize*i] = dataSrc[i*2  ];
                          imagDst[imgChSize*i] = dataSrc[i*2+1];
                      }
                  }
                      break;
              }
          }
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Creates the struct which inits & contains the OpenCL context & environment.
   Inits wired(?) space to buffer the image in host RAM. Chooses the OpenCL
   device, perhaps the user can choose it later? This would also choose the
   appropriate OpenCL image format (R, RG, RGBA, or multiples thereof). Space
   for metadata can be allocated as required, though.
  
   Supported image formats are:
   CL_FLOAT, CL_SNORM_INT8, CL_UNORM_INT8, CL_SNORM_INT16, CL_UNORM_INT16
   32-bit int formats won't keep precision when converted to floats internally
   and doubles are generally not supported on the GPU image formats.
   */
  struct oclWarper* GDALWarpKernelOpenCL_createEnv(int srcWidth, int srcHeight,
                                                   int dstWidth, int dstHeight,
                                                   cl_channel_type imageFormat,
                                                   int numBands, int coordMult,
                                                   int useImag, int useBandSrcValid,
                                                   CPL_UNUSED float *fDstDensity,
                                                   double *dfDstNoDataReal,
                                                   OCLResampAlg resampAlg, cl_int *clErr)
  {
      struct oclWarper *warper;
      int i;
      size_t maxWidth = 0, maxHeight = 0;
      cl_int err = CL_SUCCESS;
      size_t fmtSize, sz;
      cl_device_id device;
      cl_bool bool_flag;
      OCLVendor eCLVendor = VENDOR_OTHER;
  
      // Do we have a suitable OpenCL device?
      device = get_device(&eCLVendor);
      if( device == NULL )
          return NULL;
  
      err = clGetDeviceInfo(device, CL_DEVICE_IMAGE_SUPPORT,
                            sizeof(cl_bool), &bool_flag, &sz);
      if( err != CL_SUCCESS || !bool_flag )
      {
          CPLDebug( "OpenCL", "No image support on selected device." );
          return NULL;
      }
  
      // Set up warper environment.
      warper = (struct oclWarper *)CPLCalloc(1, sizeof(struct oclWarper));
  
      warper->eCLVendor = eCLVendor;
  
      //Init passed vars
      warper->srcWidth = srcWidth;
      warper->srcHeight = srcHeight;
      warper->dstWidth = dstWidth;
      warper->dstHeight = dstHeight;
  
      warper->coordMult = coordMult;
      warper->numBands = numBands;
      warper->imageFormat = imageFormat;
      warper->resampAlg = resampAlg;
  
      warper->useUnifiedSrcDensity = FALSE;
      warper->useUnifiedSrcValid = FALSE;
      warper->useDstDensity = FALSE;
      warper->useDstValid = FALSE;
  
      warper->imagWorkCL = NULL;
      warper->dstImagWorkCL = NULL;
      warper->useBandSrcValidCL = NULL;
      warper->useBandSrcValid = NULL;
      warper->nBandSrcValidCL = NULL;
      warper->nBandSrcValid = NULL;
      warper->fDstNoDataRealCL = NULL;
      warper->fDstNoDataReal = NULL;
      warper->kern1 = NULL;
      warper->kern4 = NULL;
  
      warper->dev = device;
  
      warper->context = clCreateContext(0, 1, &(warper->dev), NULL, NULL, &err);
      handleErrGoto(err, error_label);
  #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6) || defined(__clang__)
  #pragma GCC diagnostic push
  #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
  #endif
      warper->queue = clCreateCommandQueue(warper->context, warper->dev, 0, &err);
  #if __GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 6) || defined(__clang__)
  #pragma GCC diagnostic pop
  #endif
      handleErrGoto(err, error_label);
  
      //Ensure that we hand handle imagery of these dimensions
      err = clGetDeviceInfo(warper->dev, CL_DEVICE_IMAGE2D_MAX_WIDTH, sizeof(size_t), &maxWidth, &sz);
      handleErrGoto(err, error_label);
      err = clGetDeviceInfo(warper->dev, CL_DEVICE_IMAGE2D_MAX_HEIGHT, sizeof(size_t), &maxHeight, &sz);
      handleErrGoto(err, error_label);
      if (maxWidth < (size_t)srcWidth || maxHeight < (size_t)srcHeight) {
          err = CL_INVALID_IMAGE_SIZE;
          handleErrGoto(err, error_label);
      }
  
      // Split bands into sets of four when possible
      // Cubic runs slower as vector, so don't use it (probably register pressure)
      // Feel free to do more testing and come up with more precise case statements
      if(numBands < 4 || resampAlg == OCL_Cubic) {
          warper->numImages = numBands;
          warper->useVec = FALSE;
      } else {
          warper->numImages = numBands/4 + numBands % 4;
          warper->useVec = TRUE;
      }
  
      //Make the pointer space for the real images
      warper->realWorkCL = (cl_mem *)CPLCalloc(sizeof(cl_mem), warper->numImages);
      warper->dstRealWorkCL = (cl_mem *)CPLCalloc(sizeof(cl_mem), warper->numImages);
  
      //Make space for the per-channel Imag data (if exists)
      if (useImag) {
          warper->imagWorkCL = (cl_mem *)CPLCalloc(sizeof(cl_mem), warper->numImages);
          warper->dstImagWorkCL = (cl_mem *)CPLCalloc(sizeof(cl_mem), warper->numImages);
      }
  
      //Make space for the per-band BandSrcValid data (if exists)
      if (useBandSrcValid) {
          //32 bits in the mask
          sz = warper->numBands * ((31 + warper->srcWidth * warper->srcHeight) >> 5);
  
          //Allocate some space for the validity of the validity mask
          err = alloc_pinned_mem(warper, 0, warper->numBands*sizeof(char),
                                 (void **)&(warper->useBandSrcValid),
                                 &(warper->useBandSrcValidCL));
          handleErrGoto(err, error_label);
  
          for (i = 0; i < warper->numBands; ++i)
              warper->useBandSrcValid[i] = FALSE;
  
          // Allocate one array for all the band validity masks.
          // Remember that the masks don't use much memory (they're bitwise).
          err = alloc_pinned_mem(warper, 0, sz * sizeof(int),
                                 (void **)&(warper->nBandSrcValid),
                                 &(warper->nBandSrcValidCL));
          handleErrGoto(err, error_label);
      }
  
      //Make space for the per-band
      if (dfDstNoDataReal != NULL) {
          alloc_pinned_mem(warper, 0, warper->numBands,
                           (void **)&(warper->fDstNoDataReal), &(warper->fDstNoDataRealCL));
  
          //Copy over values
          for (i = 0; i < warper->numBands; ++i)
              warper->fDstNoDataReal[i] = (float) dfDstNoDataReal[i];
      }
  
      //Alloc working host image memory
      //We'll be copying into these buffers soon
      switch (imageFormat) {
        case CL_FLOAT:
          err = alloc_working_arr(warper, sizeof(float *), sizeof(float), &fmtSize);
          break;
        case CL_SNORM_INT8:
          err = alloc_working_arr(warper, sizeof(char *), sizeof(char), &fmtSize);
          break;
        case CL_UNORM_INT8:
          err = alloc_working_arr(warper, sizeof(unsigned char *), sizeof(unsigned char), &fmtSize);
          break;
        case CL_SNORM_INT16:
          err = alloc_working_arr(warper, sizeof(short *), sizeof(short), &fmtSize);
          break;
        case CL_UNORM_INT16:
          err = alloc_working_arr(warper, sizeof(unsigned short *), sizeof(unsigned short), &fmtSize);
          break;
      }
      handleErrGoto(err, error_label);
  
      // Find a good & compatible image channel order for the Lat/Long array.
      err = set_supported_formats(warper, 2,
                                  &(warper->xyChOrder), &(warper->xyChSize),
                                  CL_FLOAT);
      handleErrGoto(err, error_label);
  
      //Set coordinate image dimensions
      warper->xyWidth  = (int)ceil(((float)warper->dstWidth  + (float)warper->coordMult-1)/(float)warper->coordMult);
      warper->xyHeight = (int)ceil(((float)warper->dstHeight + (float)warper->coordMult-1)/(float)warper->coordMult);
  
      //Alloc coord memory
      sz = sizeof(float) * warper->xyChSize * warper->xyWidth * warper->xyHeight;
      err = alloc_pinned_mem(warper, 0, sz, (void **)&(warper->xyWork),
                             &(warper->xyWorkCL));
      handleErrGoto(err, error_label);
  
      //Ensure everything is finished allocating, copying, & mapping
      err = clFinish(warper->queue);
      handleErrGoto(err, error_label);
  
      (*clErr) = CL_SUCCESS;
      return warper;
  
  error_label:
      GDALWarpKernelOpenCL_deleteEnv(warper);
      return NULL;
  }
  
  /*
   Copy the validity mask for an image band to the warper.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  cl_int GDALWarpKernelOpenCL_setSrcValid(struct oclWarper *warper,
                                          int *bandSrcValid, int bandNum)
  {
      //32 bits in the mask
      int stride = (31 + warper->srcWidth * warper->srcHeight) >> 5;
  
      //Copy bandSrcValid
      assert(warper->nBandSrcValid != NULL);
      memcpy(&(warper->nBandSrcValid[bandNum*stride]), bandSrcValid, sizeof(int) * stride);
      warper->useBandSrcValid[bandNum] = TRUE;
  
      return CL_SUCCESS;
  }
  
  /*
   Sets the source image real & imag into the host memory so that it is
   permuted (ex. RGBA) for better graphics card access.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  cl_int GDALWarpKernelOpenCL_setSrcImg(struct oclWarper *warper, void *imgData,
                                        int bandNum)
  {
      void **imagWorkPtr = NULL;
  
      if (warper->imagWorkCL != NULL)
          imagWorkPtr = warper->imagWork.v;
  
      return set_img_data(warper, imgData, warper->srcWidth, warper->srcHeight,
                          TRUE, bandNum, warper->realWork.v, imagWorkPtr);
  }
  
  /*
   Sets the destination image real & imag into the host memory so that it is
   permuted (ex. RGBA) for better graphics card access.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  cl_int GDALWarpKernelOpenCL_setDstImg(struct oclWarper *warper, void *imgData,
                                        int bandNum)
  {
      void **dstImagWorkPtr = NULL;
  
      if (warper->dstImagWorkCL != NULL)
          dstImagWorkPtr = warper->dstImagWork.v;
  
      return set_img_data(warper, imgData, warper->dstWidth, warper->dstHeight,
                          FALSE, bandNum, warper->dstRealWork.v, dstImagWorkPtr);
  }
  
  /*
   Inputs the source coordinates for a row of the destination pixels. Invalid
   coordinates are set as -99.0, which should be out of the image bounds. Sets
   the coordinates as ready to be used in OpenCL image memory: interleaved and
   minus the offset. By using image memory, we can use a smaller texture for
   coordinates and use OpenCL's built-in interpolation to save memory.
  
   What it does: generates a smaller matrix of X/Y coordinate transformation
   values from an original matrix. When bilinearly sampled in the GPU hardware,
   the generated values are as close as possible to the original matrix.
  
   Complication: matrices have arbitrary dimensions and the sub-sampling factor
   is an arbitrary integer greater than zero. Getting the edge cases right is
   difficult.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  cl_int GDALWarpKernelOpenCL_setCoordRow(struct oclWarper *warper,
                                          double *rowSrcX, double *rowSrcY,
                                          double srcXOff, double srcYOff,
                                          int *success, int rowNum)
  {
      int coordMult = warper->coordMult;
      int width = warper->dstWidth;
      int height = warper->dstHeight;
      int xyWidth = warper->xyWidth;
      int i;
      int xyChSize = warper->xyChSize;
      float *xyPtr, *xyPrevPtr = NULL;
      int lastRow = rowNum == height - 1;
      double dstHeightMod = 1.0, dstWidthMod = 1.0;
  
      //Return if we're at an off row
      if(!lastRow && rowNum % coordMult != 0)
          return CL_SUCCESS;
  
      //Standard row, adjusted for the skipped rows
      xyPtr = &(warper->xyWork[xyWidth * xyChSize * rowNum / coordMult]);
  
      //Find our row
      if(lastRow){
          //Setup for the final row
          xyPtr     = &(warper->xyWork[xyWidth * xyChSize * (warper->xyHeight - 1)]);
          xyPrevPtr = &(warper->xyWork[xyWidth * xyChSize * (warper->xyHeight - 2)]);
  
          if((height-1) % coordMult)
              dstHeightMod = (double)coordMult / (double)((height-1) % coordMult);
      }
  
      //Copy selected coordinates
      for (i = 0; i < width; i += coordMult) {
          if (success[i]) {
              xyPtr[0] = (float) (rowSrcX[i] - srcXOff);
              xyPtr[1] = (float) (rowSrcY[i] - srcYOff);
  
              if(lastRow) {
                  //Adjust bottom row so interpolator returns correct value
                  xyPtr[0] = (float) (dstHeightMod * (xyPtr[0] - xyPrevPtr[0]) + xyPrevPtr[0]);
                  xyPtr[1] = (float) (dstHeightMod * (xyPtr[1] - xyPrevPtr[1]) + xyPrevPtr[1]);
              }
          } else {
              xyPtr[0] = -99.0f;
              xyPtr[1] = -99.0f;
          }
  
          xyPtr += xyChSize;
          xyPrevPtr += xyChSize;
      }
  
      //Copy remaining coordinate
      if((width-1) % coordMult){
          dstWidthMod = (double)coordMult / (double)((width-1) % coordMult);
          xyPtr -= xyChSize;
          xyPrevPtr -= xyChSize;
      } else {
          xyPtr -= xyChSize*2;
          xyPrevPtr -= xyChSize*2;
      }
  
      if(lastRow) {
          double origX = rowSrcX[width-1] - srcXOff;
          double origY = rowSrcY[width-1] - srcYOff;
          double a = 1.0, b = 1.0;
  
          // Calculate the needed x/y values using an equation from the OpenCL Spec
          // section 8.2, solving for Ti1j1
          if((width -1) % coordMult)
              a = ((width -1) % coordMult)/(double)coordMult;
  
          if((height-1) % coordMult)
              b = ((height-1) % coordMult)/(double)coordMult;
  
          xyPtr[xyChSize  ] = (float) ((((1.0 - a) * (1.0 - b) * xyPrevPtr[0]
                                + a * (1.0 - b) * xyPrevPtr[xyChSize]
                                + (1.0 - a) * b * xyPtr[0]) - origX)/(-a * b));
  
          xyPtr[xyChSize+1] = (float) ((((1.0 - a) * (1.0 - b) * xyPrevPtr[1]
                                + a * (1.0 - b) * xyPrevPtr[xyChSize+1]
                                + (1.0 - a) * b * xyPtr[1]) - origY)/(-a * b));
      } else {
          //Adjust last coordinate so interpolator returns correct value
          xyPtr[xyChSize  ] = (float) (dstWidthMod * (rowSrcX[width-1] - srcXOff - xyPtr[0]) + xyPtr[0]);
          xyPtr[xyChSize+1] = (float) (dstWidthMod * (rowSrcY[width-1] - srcYOff - xyPtr[1]) + xyPtr[1]);
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Copies all data to the device RAM, frees the host RAM, runs the
   appropriate resampling kernel, mallocs output space, & copies the data
   back from the device RAM for each band. Also check to make sure that
   setRow*() was called the appropriate number of times to init all image
   data.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  cl_int GDALWarpKernelOpenCL_runResamp(struct oclWarper *warper,
                                        float *unifiedSrcDensity,
                                        unsigned int *unifiedSrcValid,
                                        float *dstDensity,
                                        unsigned int *dstValid,
                                        double dfXScale, double dfYScale,
                                        double dfXFilter, double dfYFilter,
                                        int nXRadius, int nYRadius,
                                        int nFiltInitX, int nFiltInitY)
  {
      int i, nextBandNum = 0, chSize = 1;
      cl_int err = CL_SUCCESS;
      cl_mem xy, unifiedSrcDensityCL, unifiedSrcValidCL;
      cl_mem dstDensityCL, dstValidCL, dstNoDataRealCL;
      cl_mem useBandSrcValidCL, nBandSrcValidCL;
      size_t groupSize, wordSize = 0;
      cl_kernel kern = NULL;
      cl_channel_order chOrder;
  
      warper->useUnifiedSrcDensity = unifiedSrcDensity != NULL;
      warper->useUnifiedSrcValid = unifiedSrcValid != NULL;
  
      //Check the word size
      switch (warper->imageFormat) {
          case CL_FLOAT:
              wordSize = sizeof(float);
              break;
          case CL_SNORM_INT8:
              wordSize = sizeof(char);
              break;
          case CL_UNORM_INT8:
              wordSize = sizeof(unsigned char);
              break;
          case CL_SNORM_INT16:
              wordSize = sizeof(short);
              break;
          case CL_UNORM_INT16:
              wordSize = sizeof(unsigned short);
              break;
      }
  
      //Compile the kernel; the invariants are being compiled into the code
      if (!warper->useVec || warper->numBands % 4) {
          warper->kern1 = get_kernel(warper, FALSE,
                                     dfXScale, dfYScale, dfXFilter, dfYFilter,
                                     nXRadius, nYRadius, nFiltInitX, nFiltInitY, &err);
          handleErr(err);
      }
      if (warper->useVec){
          warper->kern4 = get_kernel(warper, TRUE,
                                     dfXScale, dfYScale, dfXFilter, dfYFilter,
                                     nXRadius, nYRadius, nFiltInitX, nFiltInitY, &err);
          handleErr(err);
      }
  
      //Copy coord data to the device
      handleErr(err = set_coord_data(warper, &xy));
  
      //Copy unified density & valid data
      handleErr(err = set_unified_data(warper, &unifiedSrcDensityCL, &unifiedSrcValidCL,
                                       unifiedSrcDensity, unifiedSrcValid,
                                       &useBandSrcValidCL, &nBandSrcValidCL));
  
      //Copy output density & valid data
      handleErr(set_dst_data(warper, &dstDensityCL, &dstValidCL, &dstNoDataRealCL,
                             dstDensity, dstValid, warper->fDstNoDataReal));
  
      //What's the recommended group size?
      if (warper->useVec) {
          // Start with the vector kernel
          handleErr(clGetKernelWorkGroupInfo(warper->kern4, warper->dev,
                                             CL_KERNEL_WORK_GROUP_SIZE,
                                             sizeof(size_t), &groupSize, NULL));
          kern = warper->kern4;
          chSize = warper->imgChSize4;
          chOrder = warper->imgChOrder4;
      } else {
          // We're only using the float kernel
          handleErr(clGetKernelWorkGroupInfo(warper->kern1, warper->dev,
                                             CL_KERNEL_WORK_GROUP_SIZE,
                                             sizeof(size_t), &groupSize, NULL));
          kern = warper->kern1;
          chSize = warper->imgChSize1;
          chOrder = warper->imgChOrder1;
      }
  
      //Loop over each image
      for (i = 0; i < warper->numImages; ++i)
      {
          cl_mem srcImag, srcReal;
          cl_mem dstReal, dstImag;
          int bandNum = nextBandNum;
  
          //Switch kernels if needed
          if (warper->useVec && nextBandNum < warper->numBands - warper->numBands % 4) {
              nextBandNum += 4;
          } else {
              if (kern == warper->kern4) {
                  handleErr(clGetKernelWorkGroupInfo(warper->kern1, warper->dev,
                                                     CL_KERNEL_WORK_GROUP_SIZE,
                                                     sizeof(size_t), &groupSize, NULL));
                  kern = warper->kern1;
                  chSize = warper->imgChSize1;
                  chOrder = warper->imgChOrder1;
              }
              ++nextBandNum;
          }
  
          //Create & copy the source image
          handleErr(err = set_src_rast_data(warper, i, chSize*wordSize, chOrder,
                                            &srcReal, &srcImag));
  
          //Create & copy the output image
          if (kern == warper->kern1) {
              handleErr(err = set_dst_rast_data(warper, i, wordSize, &dstReal, &dstImag));
          } else {
              handleErr(err = set_dst_rast_data(warper, i, wordSize*4, &dstReal, &dstImag));
          }
  
          //Set the bandNum
          handleErr(err = clSetKernelArg(kern, 12, sizeof(int), &bandNum));
  
          //Run the kernel
          handleErr(err = execute_kern(warper, kern, groupSize));
  
          //Free loop CL mem
          handleErr(err = clReleaseMemObject(srcReal));
          handleErr(err = clReleaseMemObject(srcImag));
  
          //Copy the back output results
          if (kern == warper->kern1) {
              handleErr(err = get_dst_rast_data(warper, i, wordSize, dstReal, dstImag));
          } else {
              handleErr(err = get_dst_rast_data(warper, i, wordSize*4, dstReal, dstImag));
          }
  
          //Free remaining CL mem
          handleErr(err = clReleaseMemObject(dstReal));
          handleErr(err = clReleaseMemObject(dstImag));
      }
  
      //Free remaining CL mem
      handleErr(err = clReleaseMemObject(xy));
      handleErr(err = clReleaseMemObject(unifiedSrcDensityCL));
      handleErr(err = clReleaseMemObject(unifiedSrcValidCL));
      handleErr(err = clReleaseMemObject(useBandSrcValidCL));
      handleErr(err = clReleaseMemObject(nBandSrcValidCL));
      handleErr(err = clReleaseMemObject(dstDensityCL));
      handleErr(err = clReleaseMemObject(dstValidCL));
      handleErr(err = clReleaseMemObject(dstNoDataRealCL));
  
      return CL_SUCCESS;
  }
  
  /*
   Sets pointers to the floating point data in the warper. The pointers
   are internal to the warper structure, so don't free() them. If the imag
   channel is in use, it will receive a pointer. Otherwise it'll be set to NULL.
   These are pointers to floating point data, so the caller will need to
   manipulate the output as appropriate before saving the data.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  cl_int GDALWarpKernelOpenCL_getRow(struct oclWarper *warper,
                                     void **rowReal, void **rowImag,
                                     int rowNum, int bandNum)
  {
      int memOff = rowNum * warper->dstWidth;
      int imgNum = bandNum;
  
      if (warper->useVec && bandNum < warper->numBands - warper->numBands % 4) {
          memOff += warper->dstWidth * warper->dstHeight * (bandNum % 4);
          imgNum = bandNum / 4;
      } else if(warper->useVec) {
          imgNum = bandNum / 4 + bandNum % 4;
      }
  
      //Return pointers into the warper's data
      switch (warper->imageFormat) {
          case CL_FLOAT:
              (*rowReal) = &(warper->dstRealWork.f[imgNum][memOff]);
              break;
          case CL_SNORM_INT8:
              (*rowReal) = &(warper->dstRealWork.c[imgNum][memOff]);
              break;
          case CL_UNORM_INT8:
              (*rowReal) = &(warper->dstRealWork.uc[imgNum][memOff]);
              break;
          case CL_SNORM_INT16:
              (*rowReal) = &(warper->dstRealWork.s[imgNum][memOff]);
              break;
          case CL_UNORM_INT16:
              (*rowReal) = &(warper->dstRealWork.us[imgNum][memOff]);
              break;
      }
  
      if (warper->dstImagWorkCL == NULL) {
          (*rowImag) = NULL;
      } else {
          switch (warper->imageFormat) {
              case CL_FLOAT:
                  (*rowImag) = &(warper->dstImagWork.f[imgNum][memOff]);
                  break;
              case CL_SNORM_INT8:
                  (*rowImag) = &(warper->dstImagWork.c[imgNum][memOff]);
                  break;
              case CL_UNORM_INT8:
                  (*rowImag) = &(warper->dstImagWork.uc[imgNum][memOff]);
                  break;
              case CL_SNORM_INT16:
                  (*rowImag) = &(warper->dstImagWork.s[imgNum][memOff]);
                  break;
              case CL_UNORM_INT16:
                  (*rowImag) = &(warper->dstImagWork.us[imgNum][memOff]);
                  break;
          }
      }
  
      return CL_SUCCESS;
  }
  
  /*
   Free the OpenCL warper environment. It should check everything for NULL, so
   be sure to mark free()ed pointers as NULL or it'll be double free()ed.
  
   Returns CL_SUCCESS on success and other CL_* errors when something goes wrong.
   */
  cl_int GDALWarpKernelOpenCL_deleteEnv(struct oclWarper *warper)
  {
      int i;
      cl_int err = CL_SUCCESS;
  
      for (i = 0; i < warper->numImages; ++i) {
          // Run free!!
          void* dummy = NULL;
          if( warper->realWork.v )
              freeCLMem(warper->realWorkCL[i], warper->realWork.v[i]);
          else
              freeCLMem(warper->realWorkCL[i], dummy);
          if( warper->realWork.v )
              freeCLMem(warper->dstRealWorkCL[i], warper->dstRealWork.v[i]);
          else
              freeCLMem(warper->dstRealWorkCL[i], dummy);
  
          //(As applicable)
          if(warper->imagWorkCL != NULL && warper->imagWork.v != NULL && warper->imagWork.v[i] != NULL) {
              freeCLMem(warper->imagWorkCL[i], warper->imagWork.v[i]);
          }
          if(warper->dstImagWorkCL != NULL && warper->dstImagWork.v != NULL && warper->dstImagWork.v[i] != NULL) {
              freeCLMem(warper->dstImagWorkCL[i], warper->dstImagWork.v[i]);
          }
      }
  
      //Free cl_mem
      freeCLMem(warper->useBandSrcValidCL, warper->useBandSrcValid);
      freeCLMem(warper->nBandSrcValidCL, warper->nBandSrcValid);
      freeCLMem(warper->xyWorkCL, warper->xyWork);
      freeCLMem(warper->fDstNoDataRealCL, warper->fDstNoDataReal);
  
      //Free pointers to cl_mem*
      if (warper->realWorkCL != NULL)
          CPLFree(warper->realWorkCL);
      if (warper->dstRealWorkCL != NULL)
          CPLFree(warper->dstRealWorkCL);
  
      if (warper->imagWorkCL != NULL)
          CPLFree(warper->imagWorkCL);
      if (warper->dstImagWorkCL != NULL)
          CPLFree(warper->dstImagWorkCL);
  
      if (warper->realWork.v != NULL)
          CPLFree(warper->realWork.v);
      if (warper->dstRealWork.v != NULL)
          CPLFree(warper->dstRealWork.v);
  
      if (warper->imagWork.v != NULL)
          CPLFree(warper->imagWork.v);
      if (warper->dstImagWork.v != NULL)
          CPLFree(warper->dstImagWork.v);
  
      //Free OpenCL structures
      if (warper->kern1 != NULL)
          clReleaseKernel(warper->kern1);
      if (warper->kern4 != NULL)
          clReleaseKernel(warper->kern4);
      if (warper->queue != NULL)
          clReleaseCommandQueue(warper->queue);
      if (warper->context != NULL)
          clReleaseContext(warper->context);
  
      CPLFree(warper);
  
      return CL_SUCCESS;
  }
  
  #endif /* defined(HAVE_OPENCL) */
  