         //std::cout << "i: " << ws[0] << ", j: " << ws[1] << ", lambda: " << lambda << std::endl;
     }
 
+    assert_cuda(cudaMemcpyFromSymbol(&cacheUpdateCnt, d_cacheUpdateCnt, sizeof(int), 0));
+    std::cout << "Cache row updates: " << cacheUpdateCnt << std::endl;
+
     assert_cuda(cudaMemcpy(alpha, d_alpha, num_vec * sizeof(float), cudaMemcpyDeviceToHost));
 
     assert_cuda(cudaFree(d_K));
     assert_cuda(cudaFree(d_KDiag));
     assert_cuda(cudaFree(d_KCacheRemapIdx));
-    assert_cuda(cudaFree(d_KCacheVecIdx));
+    assert_cuda(cudaFree(d_KCacheRowIdx));
+    assert_cuda(cudaFree(d_KCacheRowPriority));
     assert_cuda(cudaFree(d_alpha));
     assert_cuda(cudaFree(d_x));
     assert_cuda(cudaFree(d_xT));

    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 = 200;
    size_t cache_rows = cache_size_mb * 1024 * 1024 / (num_vec_aligned * sizeof(float));
    cache_rows = std::min(cache_rows, num_vec);


        for (int d = 0; d < dim; d++)
        {
            //float diff = xi[d] - x[dim_aligned * j + d];
            float diff = sx[d] - xT[num_vec_aligned * d + j];
            //float diff = x[dim_aligned * i + d] - x[dim_aligned * j + d];
            //float diff = xT[num_vec_aligned * d + i] - xT[num_vec_aligned * d + j];
            sum += diff * diff;
        }
        K[num_vec_aligned * last + j] = exp(-gamma * sum);
        j += gridDim.x * blockDim.x;
    }
}

__global__ void kernelCheckCache(const int * i_ptr, float * K, int * KCacheRemapIdx, int * KCacheRowIdx, int * KCacheRowPriority, int cache_rows, const float * x, const float * xT, float gamma, int num_vec, int num_vec_aligned, int dim, int dim_aligned, int lastPtrIdx)
{
    extern __shared__ int2 spriority[];
    int i = *i_ptr;
    int j = blockDim.x * blockIdx.x + threadIdx.x;
    if (KCacheRemapIdx[i] >= 0)
    {
        if (j == 0)
            KCacheRowPriority[KCacheRemapIdx[i]] = d_cacheUpdateCnt;  // refresh priority
        return;  //item already in cache
    }
    //int last = (KCacheRowIdx[cache_rows + lastPtrIdx] + 1) % cache_rows;
    int2 minpriority = make_int2(INT_MAX, 0);
    for (int k = 0; k < cache_rows; k += blockDim.x)
    {
        int idx = k + threadIdx.x;
        if (idx < cache_rows)
        {
            int v = KCacheRowPriority[idx];
            if (v < minpriority.x)
                minpriority = make_int2(v, idx);
        }
    }
    spriority[threadIdx.x] = minpriority;
    __syncthreads();
    for (int s = blockDim.x / 2; s > 0; s >>= 1)
    {
        if (threadIdx.x < s)
        {
            if (spriority[threadIdx.x + s].x < spriority[threadIdx.x].x)
                spriority[threadIdx.x] = spriority[threadIdx.x + s];
        }
        __syncthreads();
    }
    int last = spriority[0].y;
    if (j == 0)
    {
        //KCacheRowIdx[cache_rows + (1 - lastPtrIdx)] = last;
        int del_i = KCacheRowIdx[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;
        KCacheRowIdx[last] = i;
        KCacheRowPriority[last] = ++d_cacheUpdateCnt;
    }

    //calculate cache matrix row [last], original index is [i]
    float * sx = (float *)spriority;
    for (int idxshift = 0; idxshift < dim; idxshift += blockDim.x)
    {
        int idx = idxshift + threadIdx.x;
        if (idx < dim)
            //xi[idx] = xT[num_vec_aligned * idx + i];
            sx[idx] = x[dim_aligned * i + idx];
    }
    __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 = sx[d] - xT[num_vec_aligned * d + j];
            //float diff = x[dim_aligned * i + d] - x[dim_aligned * j + d];
            //float diff = xT[num_vec_aligned * d + i] - xT[num_vec_aligned * d + j];
            sum += diff * diff;
