     alpha[j] -= l * y[j];
 }
 
+__global__ void kernelUpdateAlphaAndLambdaCached(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, const float * KDiag, const int * KCacheRemapIdx)
+{
+    int i = ws[0];
+    int j = ws[1];
+    int cache_row = KCacheRemapIdx[i];
+    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]) / (KDiag[i] + KDiag[j] - 2 * K[num_vec_aligned * cache_row + j]);
+    float l = min(l1, min(l2, l3));
+
+    *lambda = l;
+    alpha[i] += l * y[i];
+    alpha[j] -= l * y[j];
+}
+
+__global__ void kernelCheckCache(const int * i_ptr, float * K, int * KCacheRemapIdx, int * KCacheVecIdx, int cache_rows, const float * x, float gamma, int num_vec, int num_vec_aligned, int dim, int dim_aligned, int lastPtrIdx)
+{
+    int i = *i_ptr;
+    if (KCacheRemapIdx[i] >= 0)
+        return;  //item already in cache
+    int j = blockDim.x * blockIdx.x + threadIdx.x;
+    int last = (KCacheVecIdx[cache_rows + lastPtrIdx] + 1) % cache_rows;
+    if (j == 0)
+    {
+        KCacheVecIdx[cache_rows + (lastPtrIdx + 1) % 2] = last;
+        int del_i = KCacheVecIdx[last];
+        if (del_i >= 0)
+            KCacheRemapIdx[del_i] = -1;  //cache row for vector [del_i] will be overwritten, remove it from RemapIdx array
+        //set correct indices
+        KCacheRemapIdx[i] = last;
+        KCacheVecIdx[last] = i;
+    }
+    //calculate cache matrix row [last], original index is [i]
+    extern __shared__ float xi[];
+    int idxshift = 0;
+    while (idxshift < dim)
+    {
+        int idx = idxshift + threadIdx.x;
+        if (idx < dim)
+            xi[idx] = x[num_vec_aligned * idx + i];
+            //xi[idx] = x[dim_aligned * i + idx];
+        idxshift += blockDim.x;
+    }
+    __syncthreads();
+    while (j < num_vec)
+    {
+        float sum = 0;
+        for (int d = 0; d < dim; d++)
+        {
+            //float diff = xi[d] - x[dim_aligned * j + d];
+            float diff = xi[d] - x[num_vec_aligned * d + j];
+            //float diff = x[dim_aligned * i + d] - x[dim_aligned * j + d];
+            sum += diff * diff;
+        }
+        K[num_vec_aligned * last + j] = exp(-gamma * sum);
+        j += gridDim.x * blockDim.x;
+    }
+}
+
 void OrcusSvmTrain(float * alpha, float * rho, const float * x, const float * y, size_t num_vec, size_t num_vec_aligned, size_t dim, size_t dim_aligned, float C, float gamma, float eps)
 {
     float *d_alpha = nullptr,
         *d_x = nullptr,
+        *d_xT = nullptr,
         *d_y = nullptr,
         *d_g = nullptr,
         *d_K = nullptr,
+        *d_KDiag = nullptr,
         *d_reduceval = nullptr,
         *d_reduceval2 = nullptr;
     int *d_reduceidx = nullptr,
         *d_reduceidx2 = nullptr;
-    float *d_lambda = nullptr,
-        *d_ones = nullptr;
-    int *d_workingset = nullptr;
+    float *d_lambda = nullptr;
+    int *d_workingset = nullptr,
+        *d_KCacheRemapIdx = nullptr,
+        *d_KCacheVecIdx = nullptr;  // item at index [cache_rows] is index of last inserted item
 
     size_t reduce_block_size = 256;
     size_t reduce_buff_size = rounduptomult(num_vec, reduce_block_size);
     size_t ones_size = std::max(num_vec_aligned, dim_aligned);
+    size_t cache_size_mb = 50;
+    size_t cache_rows = cache_size_mb * 1024 * 1024 / (num_vec_aligned * sizeof(float));
+    cache_rows = std::min(cache_rows, num_vec);
+
+    std::cout << "Cache size: " << cache_rows << " rows (" << (100.f * cache_rows / (float)num_vec) << " % of data set)" << std::endl;
 
     cublasHandle_t cublas;
     assert_cublas(cublasCreate(&cublas));
 
     assert_cuda(cudaMalloc(&d_alpha, num_vec_aligned * sizeof(float)));
     assert_cuda(cudaMalloc(&d_x, num_vec_aligned * dim_aligned * sizeof(float)));
+    assert_cuda(cudaMalloc(&d_xT, num_vec_aligned * dim_aligned * sizeof(float)));
     assert_cuda(cudaMalloc(&d_y, num_vec_aligned * sizeof(float)));
     assert_cuda(cudaMalloc(&d_g, num_vec_aligned * sizeof(float)));
     assert_cuda(cudaMalloc(&d_reduceval, reduce_buff_size * sizeof(float)));

    assert_cuda(cudaMalloc(&d_reduceidx2, reduce_buff_size / reduce_block_size * sizeof(int)));
    assert_cuda(cudaMalloc(&d_lambda, sizeof(float)));
    assert_cuda(cudaMalloc(&d_workingset, 2 * sizeof(int)));
    assert_cuda(cudaMalloc(&d_K, num_vec_aligned * num_vec_aligned * sizeof(float)));

    assert_cuda(cudaMemset(d_alpha, 0, num_vec_aligned * sizeof(float)));
    assert_cuda(cudaMemcpy(d_x, x, num_vec_aligned * dim_aligned * sizeof(float), cudaMemcpyHostToDevice));
    assert_cuda(cudaMemcpy(d_y, y, num_vec_aligned * sizeof(float), cudaMemcpyHostToDevice));

    //export_cuda_buffer(d_x, dim_aligned, num_vec_aligned, sizeof(float), "x.dat");

    dim3 dimBlock(256);
    dim3 dimGrid(getgriddim(num_vec_aligned, (size_t)dimBlock.x));
    kernelInitg<<<dimGrid, dimBlock>>>(d_g, num_vec_aligned);
    assert_cuda(cudaDeviceSynchronize());

    //export_cuda_buffer(d_g, num_vec_aligned, 1, sizeof(float), "g.dat");

    cudaEvent_t evstart, evend;
    cudaEventCreate(&evstart);
    cudaEventCreate(&evend);
    cudaEventRecord(evstart);

    computeK(d_K, d_x, gamma, num_vec, num_vec_aligned, dim, dim_aligned, cublas);
    //dimBlock = dim3(16, 16);
    //dimGrid = dim3(getgriddim(num_vec_aligned, (size_t)dimBlock.x), getgriddim(num_vec_aligned, (size_t)dimBlock.y));
    //kernelComputeK<<<dimGrid, dimBlock>>>(d_K, d_x, gamma, num_vec, num_vec_aligned, dim, dim_aligned);

    cudaEventRecord(evend);
    assert_cuda(cudaDeviceSynchronize());
    float t;
    cudaEventElapsedTime(&t, evstart, evend);
    std::cout << "Kernel time: " << t << std::endl;
    cudaEventDestroy(evstart);
    cudaEventDestroy(evend);

    //export_cuda_buffer(d_K, num_vec_aligned, num_vec_aligned, sizeof(float), "K.dat");

    dimBlock = dim3(reduce_block_size);
    dimGrid = dim3(getgriddim(num_vec_aligned, (size_t)dimBlock.x));
    for (int iter = 0;; iter++)
    {
        kernelSelectI<<<dimGrid, dimBlock>>>(d_reduceval, d_reduceidx, d_y, d_g, d_alpha, C, num_vec);

    }
}

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)
