 // Use grid strided loops, descriped here:
 // https://devblogs.nvidia.com/cuda-pro-tip-write-flexible-kernels-grid-stride-loops/
+// This pattern ensures that all of the loop values are visited once, no matter
+// what grid parameters are used for the function.
 	
 extern "C" __global__
 void seq2col(float* output,



extern "C" __global__
void maxout(float* best, int* which,
        const float* cands, int B, int O, int P)
{
    for (int b = blockIdx.x * blockDim.x + threadIdx.x; b < B; b += blockDim.x * gridDim.x)
    {
        // Go to the regions we're working on
        for (int i=0; i < b; ++i) {
            best += O;
            which += O;
            cands += O * P;
        }

        for (int i=0; i < O; ++i)
        {
            which[i] = 0
            best[i] = cands[0];
            for (int p=1; p < P; ++p)
            {
                if (cands[i+p] > best[0])
                {
                    which[i] = p;
                    best[i] = cands[i+p];
                }
            }
        }
    }
}


extern "C" __global__
void max_pool(float* maxes, int* which,
    const float* X, const int* lengths, int B, int T, int O)
{

    for (int b = blockIdx.x * blockDim.x + threadIdx.x; b < B; b += blockDim.x * gridDim.x)
    {

        // Go to the regions we're working on
        for (int i=0; i < b; ++i) {
            maxes += O;
            which += O;
            X += lengths[i] * O;
        }
 
        // Each invocation of the kernel maxes one batch.
        // Start by assuming maxes are at i=0
        for (int j=0; j < O; ++j) {
            maxes[j] = X[j];
            which[j] = 0;
        }
        X += O;
        
        int length = lengths[b];
        for (int i=1; i < length; ++i) // Iterate over rows
        {
            for (int j=0; j < O; ++j)
            {
                if (X[j] > maxes[j])
                {
                    maxes[j] = X[j];
                    which[j] = i;
                }
            }
            X += O;
        }
    }
}


    //    d_seq[i] += d_cols[i+1, 1]
    //    d_seq[i] += d_cols[i+2, 0]

    int _loop_start = blockIdx.x * blockDim.x + threadIdx.x;
    int _loop_stride = blockDim.x * gridDim.x;
    int nF = nW * 2 + 1;
    int end_d_cols = B * I * nF;
    for (int b = _loop_start; b < B; b += _loop_stride)
    {
        float* d_seqs_b = &d_seqs[b*I];
        int col_feat = nF * I;
        for (int f=-nW; f < (nW+1); ++f)
        {
            int col_row = (b+f) * (I*nF);
            col_feat -= I;
            if ((col_row >= 0) && (col_row < end_d_cols))
            {
                int start = col_row + col_feat;
                if ((start >= 0) && ((start+I) < end_d_cols))
                {
                    for (int i=0; i < I; ++i)
                        d_seqs_b[i] += d_cols[start+i];
                }
            }
        }
