         for (int c = 0; c < colorsPerThread; c++) {
             #pragma unroll
             for (int i = 0; i < imgsPerThread; i++) {
-                targets[c * B_Y * imgSizeY * imgSizeX * numImages + i * B_X] = scaleOutputs * prod[c][i];
+                if (!checkCaseBounds || myCaseIdx + i * B_X < numImages) {
+                    targets[c * B_Y * imgPixels * numImages + i * B_X] = scaleOutputs * prod[c][i];
+                }
             }
         }
     }

//    const float* fLoad = conv ? &filters[pxIdxInFilter * numFilters + 0]
//                              : &filters[moduleIdx * numFilterColors * filterPixels * numFilters + pxIdxInFilter * numFilters + 0];
    int filtersLoadOffset = filtersOffset + (conv ? pxIdxInFilter * numFilters + 0
                                                  : moduleIdx * numFilterColors * filterPxlAll + pxIdxInFilter * numFilters);
    // preload the filter's pixel
    #pragma unroll
    for (int i = 0; i < 4; i++) {
        // discarding redundant conditions  
        wPreload[i] = tex1Dfetch<float>(filters, filtersLoadOffset + i * 16 * filterPxlAll);
    }

//    const float* hLoad = &hidActs[(moduleIdx + 0 * numModules) * numImages];
    int hidActsLoadOffset = hidActsOffset + moduleIdx * numImages;
    
    // preload the hidAct's pixel
    #pragma unroll
    for (int x = 0; x < 4; x++) {
        hPreload[x]   = tex1Dfetch<float>(hidActs, hidActsLoadOffset + x * 32); 
        hPreload[x+4] = tex1Dfetch<float>(hidActs, hidActsLoadOffset + 8 * numImgMods + x * 32); 
    }

    for (int my = startY; my < endY; my++) {

        for (int mx = startX; mx < endX; mx++) {
            conv_img_acts_manycolor_preload_ty_8_tx_32_c_8_ff_32_fh_16_setCoords(my+(mx+1 == endX), (mx+1 == endX)?startX:(mx+1), numModulesX, paddingStart, moduleStride, blockPixelIdxY, blockPixelIdxX, filterSize, moduleIdx, pxIdxInFilter);

            for (int f = 0; f < numFiltersPerGroup; f += filterCacheF) { // multiply with filterCacheF filters at a time
                // As 'filterCacheF' is two times of 'filterCacheH', so divide the procedure into two steps, which loads less data during each step
                // with this division, program would not stall at memory access, and get higher gpu occupancy
                for (int k = 0; k < filterCacheF / filterCacheH; k++) {
                    // store the preloaded filters's pixel into shared memory with 64 bit width
                    shFilters[(ty * B_X + tx) % 16][(ty * B_X + tx) / 16].x = wPreload[0];
                    shFilters[(ty * B_X + tx) % 16][(ty * B_X + tx) / 16].y = wPreload[2];
                    shFilters[(ty * B_X + tx) % 16][(ty * B_X + tx) / 16 + 16].x = wPreload[1];
                    shFilters[(ty * B_X + tx) % 16][(ty * B_X + tx) / 16 + 16].y = wPreload[3];

                    // store the preloaded hidacts' pixel into shared memory with 64 bit width
                    shHidActs[ty][tx].x = hPreload[0];
                    shHidActs[ty][tx].y = hPreload[1];
                    shHidActs[ty][tx + B_X].x = hPreload[2];
                    shHidActs[ty][tx + B_X].y = hPreload[3];
                    shHidActs[ty + filterCacheH/2][tx].x = hPreload[4];
                    shHidActs[ty + filterCacheH/2][tx].y = hPreload[5];
                    shHidActs[ty + filterCacheH/2][tx + B_X].x = hPreload[6];
                    shHidActs[ty + filterCacheH/2][tx + B_X].y = hPreload[7];

                    // update to the 2nd 'filterCacheH' 
                    filtersLoadOffset += filterCacheH;
                    hidActsLoadOffset += filterCacheH * numImgMods;

                    // move to the next pixel
                    if ((k == 1) && (f == numFiltersPerGroup - filterCacheF))  {
                        filtersLoadOffset = filtersOffset + (conv ? pxIdxInFilter * numFilters
                                            : moduleIdx * numFilterColors * filterPxlAll + pxIdxInFilter * numFilters);
                        hidActsLoadOffset = hidActsOffset + moduleIdx * numImages;
                    }

                    // preloade filter's pixel 
                    #pragma unroll 4
                    for (int x = 0; x < 4; x++) {
                        wPreload[x] = tex1Dfetch<float>(filters, filtersLoadOffset + x * 16 * filterPxlAll);
                    }

                    __syncthreads();


                    // put together the instructions with same type to improve instruction-level parallelism 
                    // calculate the derivative of the hidAct with respect to image
                    #pragma unroll 
                    for (int w = 0; w < filterCacheH; w++) {
                        #pragma unroll
                        for (int c = 0; c < colorsPerThread / 2; c++) { 
                            prod[c][0] += shFilters[w][c * B_Y + ty].x * shHidActs[w][tx].x; 
                            prod[c][1] += shFilters[w][c * B_Y + ty].x * shHidActs[w][tx].y; 
                            prod[c][2] += shFilters[w][c * B_Y + ty].x * shHidActs[w][tx+B_X].x; 
                            prod[c][3] += shFilters[w][c * B_Y + ty].x * shHidActs[w][tx+B_X].y; 
                            prod[c + colorsPerThread / 2][0] += shFilters[w][c * B_Y + ty].y * shHidActs[w][tx].x;
                            prod[c + colorsPerThread / 2][1] += shFilters[w][c * B_Y + ty].y * shHidActs[w][tx].y; 
                            prod[c + colorsPerThread / 2][2] += shFilters[w][c * B_Y + ty].y * shHidActs[w][tx+B_X].x; 
                            prod[c + colorsPerThread / 2][3] += shFilters[w][c * B_Y + ty].y * shHidActs[w][tx+B_X].y; 
                        } 
                    }

                    // preloade the hidAct's pixel
                    #pragma unroll
                    for (int x = 0; x < 4; x++) {
                        hPreload[x]   = tex1Dfetch<float>(hidActs, hidActsLoadOffset + x * B_X); 
                        hPreload[x+4] = tex1Dfetch<float>(hidActs, hidActsLoadOffset + 8 * numImgMods + x * B_X); 
                    }

                    __syncthreads();
                }
            }
        }
    }
    
    // remove redundant conditions for less registers
    if (scale) {
        #pragma unroll
        for (int c = 0; c < colorsPerThread; c++) {
            #pragma unroll
            for (int i = 0; i < imgsPerThread; i++) {
                targets[c * B_Y * imgSizeY * imgSizeX * numImages + i * B_X] = scaleTargets * targets[c * B_Y * imgSizeY * imgSizeX * numImages + i * B_X] + scaleOutputs * prod[c][i];
            }
        }
    } else {

                                          const int filterSize, const int imgSizeY, const int imgSizeX, const int paddingStart, const int moduleStride,
                                          const int numImgColors, const int numGroups,
                                          const float scaleTargets, const float scaleOutputs) {
    __shared__ float shFilters[colorsPerThread*B_Y][filterCacheF];
    __shared__ float shHidActs[filterCacheH][B_X*imgsPerThread];

    const int numImgBlocks = DIVUP(numImages,B_X*imgsPerThread);
    const int blockCaseIdx = (blockIdx.x % numImgBlocks) * B_X*imgsPerThread;
    const int myCaseIdx = blockCaseIdx + threadIdx.x;

    const int imgColorIdx = (blockIdx.x / numImgBlocks) * B_Y*colorsPerThread; // color idx globally
    const int numFilterColors = numImgColors / numGroups;
    const int blockGroupIdx = imgColorIdx / numFilterColors;
    const int filterColorIdx = imgColorIdx % numFilterColors; // color idx within group
    const int numFiltersPerGroup = numFilters / numGroups;
    const int blockFilterIdx = blockGroupIdx * numFiltersPerGroup;

    const int blockPixelIdx = blockIdx.y;
    const int blockPixelIdxX = blockPixelIdx % imgSizeX;
    const int blockPixelIdxY = blockPixelIdx / imgSizeX;

    const int filterPixels = filterSize * filterSize;
    const int imgPixels = imgSizeY * imgSizeX;
    const int tidx = threadIdx.y * B_X + threadIdx.x;
//    const int hidActLoadY = threadIdx.y % B_Y, hidActLoadX = threadIdx.x % B_X;
    //const int hidActLoadY = tidx / (B_X*imgsPerThread), hidActLoadX = tidx % (B_X*imgsPerThread);
    const int filtersLoadY = tidx / filterCacheF, filtersLoadX = tidx % filterCacheF;
    // nvcc is behaving idiotically again, these useless declarations save registers
    //const int outputY = threadIdx.y, outputX = threadIdx.x;
    //const int ty = threadIdx.y, tx = threadIdx.x;
    const int numModules = numModulesY * numModulesX;
    const int hidActsOffset = (blockFilterIdx + threadIdx.y) * numImages * numModules + myCaseIdx;
    const int filtersOffset = blockFilterIdx + (filterColorIdx + filtersLoadY) * filterPixels * numFilters + filtersLoadX;
//    hidActs += (blockFilterIdx + threadIdx.y) * numImages * numModules + myCaseIdx;
//    filters += blockFilterIdx + (filterColorIdx + filtersLoadY) * filterPixels * numFilters + filtersLoadX;
    targets += (imgColorIdx + threadIdx.y) * imgPixels * numImages + blockPixelIdx * numImages + myCaseIdx;

    float prod[colorsPerThread][imgsPerThread];
    #pragma unroll
