 }
 
 //second order search
-__global__ void kernelSelectJ(float * valbuf, int * idxbuf, const float * y, const float * g, const float * alpha, float C, int num_vec, int num_vec_aligned, const int * pi, const float * K)
+__global__ void kernelSelectJ(float * valbuf, int * idxbuf, const float * y, const float * g, const float * alpha, float C, int num_vec, int num_vec_aligned, const int * i_ptr, const float * K)
 {
     int k = blockDim.x * blockIdx.x + threadIdx.x;
 
     if (k < num_vec)
     {
-        int i = *pi;
+        int i = *i_ptr;
         float y_ = y[k];
         float a_ = alpha[k];
         float th = y[i] * g[i];
         if (((y_ == 1 && a_ > 0) || (y_ == -1 && a_ < C)) && th > y[k] * g[k])
         {
             float den = K[num_vec_aligned * i + i] + K[num_vec_aligned * k + k] - 2 * K[num_vec_aligned * i + k];
+            float v = th - y[k] * g[k];
+            valbuf[k] = v * v / den;
+        }
+        else
+            valbuf[k] = -FLT_MAX;
+        idxbuf[k] = k;
+    }
+    else
+        valbuf[k] = -FLT_MAX;
+}
+
+//second order search with cached K
+__global__ void kernelSelectJCached(float * valbuf, int * idxbuf, const float * y, const float * g, const float * alpha, float C, int num_vec, int num_vec_aligned, const int * i_ptr, const float * K, const float * KDiag, const int * KCacheRemapIdx)
+{
+    int k = blockDim.x * blockIdx.x + threadIdx.x;
+
+    if (k < num_vec)
+    {
+        int i = *i_ptr;
+        int cache_row = KCacheRemapIdx[i];
+        float y_ = y[k];
+        float a_ = alpha[k];
+        float th = y[i] * g[i];
+        if (((y_ == 1 && a_ > 0) || (y_ == -1 && a_ < C)) && th > y[k] * g[k])
+        {
+            float den = KDiag[i] + KDiag[k] - 2 * K[num_vec_aligned * cache_row + k];
             float v = th - y[k] * g[k];
             valbuf[k] = v * v / den;
         }

    assert_cuda(cudaFree(d_ones));
}

__global__ void kernelSelectI(float * valbuf, int * idxbuf, const float * y, const float * g, const float * alpha, float C, int num_vec)
{
    int k = blockDim.x * blockIdx.x + threadIdx.x;

    }
}

template<typename T>
void memsetCuda(T * d_mem, T v, int n)
{
    dim3 dimBlock(256);
    dim3 dimGrid(std::min(2048, getgriddim<int>(n, dimBlock.x)));
    kernelMemset<T><<<dimGrid, dimBlock>>>(d_mem, v, n);
}

__global__ void kernelComputeK(float * K, const float * x, float gamma, int num_vec, int num_vec_aligned, int dim, int dim_aligned)
