         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 = xi[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);

        assert_cuda(cudaMemcpy(d_workingset + 1, d_reduceidx, sizeof(int), cudaMemcpyDeviceToDevice));

        //check if J is cached
        kernelCheckCache<<<dimGridCache, dimBlockCache, kernelCheckCacheSMSize>>>(d_workingset + 1, d_K, d_KCacheRemapIdx, d_KCacheVecIdx, cache_rows, d_xT, gamma, num_vec, num_vec_aligned, dim, dim_aligned, cacheLastPtrIdx);
        cacheLastPtrIdx = (cacheLastPtrIdx + 1) % 2;
        //workaround if caching J deleted I out of cache
        kernelCheckCache<<<dimGridCache, dimBlockCache, kernelCheckCacheSMSize>>>(d_workingset, d_K, d_KCacheRemapIdx, d_KCacheVecIdx, cache_rows, d_xT, gamma, num_vec, num_vec_aligned, dim, dim_aligned, cacheLastPtrIdx);
        cacheLastPtrIdx = (cacheLastPtrIdx + 1) % 2;

        if (iter > 0 && iter % 1000 == 0)

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

