             }
         }
 
-        kernelUpdateAlphaAndLambda<<<1, 1>>>(d_alpha, d_lambda, d_y, d_g, d_K, C, d_workingset, num_vec, num_vec_aligned);
-        kernelUpdateg<<<dimGrid, dimBlock>>>(d_g, d_lambda, d_y, d_K, d_workingset, num_vec, num_vec_aligned);
+        //kernelUpdateAlphaAndLambda<<<1, 1>>>(d_alpha, d_lambda, d_y, d_g, d_K, C, d_workingset, num_vec, num_vec_aligned);
+        //kernelUpdateg<<<dimGrid, dimBlock>>>(d_g, d_lambda, d_y, d_K, d_workingset, num_vec, num_vec_aligned);
+        kernelUpdateAlphaAndLambdaCached<<<1, 1>>>(d_alpha, d_lambda, d_y, d_g, d_K, C, d_workingset, num_vec, num_vec_aligned, d_KDiag, d_KCacheRemapIdx);
+        kernelUpdategCached<<<dimGrid, dimBlock>>>(d_g, d_lambda, d_y, d_K, d_workingset, num_vec, num_vec_aligned, d_KCacheRemapIdx);
 
         //float lambda;
         //int ws[2];

    }
}

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)

    assert_cuda(cudaFree(d_ones));
}

void computeKDiag(float * d_KDiag, int num_vec)
{
    //K[i,i] is always 1 for RBF kernel, let's just use memset here
    memsetCuda<float>(d_KDiag, 1, num_vec);
}

__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;
