     //export_cuda_buffer(d_idx, 1, orig_len, sizeof(int), "reduceidx.dat");
 }
 
-__global__ void kernelComputeLambda(float * lambda, const float * y, const float * g, const float * K, const float * alpha, float C, const int * ws, int num_vec_aligned)
+__global__ void kernelUpdateg(float * g, const float * lambda, const float * y, const float * K, const int * ws, int num_vec, int num_vec_aligned)
 {
     int i = ws[0];
     int j = ws[1];
-    float l1 = y[i] > 0 ? C - alpha[i] : alpha[i];
-    float l2 = y[j] > 0 ? alpha[j] : C - alpha[j];
-    float l3 = (y[i] * g[i] - y[j] * g[j]) / (K[num_vec_aligned * i + i] + K[num_vec_aligned * j + j] - 2 * K[num_vec_aligned * i + j]);
-    *lambda = min(l1, min(l2, l3));
-}
-
-__global__ void kernelUpdateg(float * g, const float * lambda, const float * y, const float * K, const int * ws, int num_vec, int num_vec_aligned)
+    int k = blockDim.x * blockIdx.x + threadIdx.x;
+    if (k < num_vec)
+    {
+        g[k] += *lambda * y[k] * (K[num_vec_aligned * j + k] - K[num_vec_aligned * i + k]);
+    }
+}
+
+__global__ void kernelUpdategCached(float * g, const float * lambda, const float * y, const float * K, const int * ws, int num_vec, int num_vec_aligned, const int * KCacheRemapIdx)
 {
     int i = ws[0];
     int j = ws[1];
+    int i_cache_row = KCacheRemapIdx[i];
+    int j_cache_row = KCacheRemapIdx[j];
     int k = blockDim.x * blockIdx.x + threadIdx.x;
     if (k < num_vec)
     {
-        g[k] += *lambda * y[k] * (K[num_vec_aligned * j + k] - K[num_vec_aligned * i + k]);
-    }
-}
-
-__global__ void kernelUpdateAlpha(float * alpha, const float * lambda, const float * y, const int * ws, int num_vec)
-{
-    int i = ws[0];
-    int j = ws[1];
-    float l = *lambda;
-    alpha[i] += l * y[i];
-    alpha[j] -= l * y[j];
+        g[k] += *lambda * y[k] * (K[num_vec_aligned * j_cache_row + k] - K[num_vec_aligned * i_cache_row + k]);
+    }
 }
 
 __global__ void kernelUpdateAlphaAndLambda(float * alpha, float * lambda, const float * y, const float * g, const float * K, float C, const int * ws, int num_vec, int num_vec_aligned)

    }
}

void cudaMemsetFloat(float * d_mem, float v, int n)
{
    dim3 dimBlock(256);
    dim3 dimGrid(std::min(2048, getgriddim<int>(n, dimBlock.x)));
    kernelMemsetFloat<<<dimGrid, dimBlock>>>(d_mem, v, n);
}

__global__ void kernelInitg(float * g, int len)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    if (i < len)
        g[i] = 1;
}

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

}

//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 * i_ptr, const float * K)
{
    int k = blockDim.x * blockIdx.x + threadIdx.x;

    if (k < num_vec)
    {
        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;
        }
