     auto means = reduce_mean(handle, points, N, NDIMS, 0);
 
     // zero center
-    add_row_vec(points, means, N, NDIMS, -1.f);
+    broadcast_matrix_vector(points, means, N, NDIMS, thrust::minus<float>(), 1, 1.f);
     
     // compute standard deviation
     auto squared_vals = square(points, N * NDIMS);

{
    pairwise_dist(handle, dist, ys, N, PROJDIM);
    // dist = (1 + ||y_i - y_j||^2)^-1
    thrust::transform(dist.begin(), dist.end(), dist.begin(), func_inc_inv());

    thrust::device_vector<float> ones(N, 1.f);
    thrust::device_vector<float> sums(N, -1.f);
    float alpha = 1.f;
    float beta = 1.f;
    cublasSafeCall(cublasSgemv(handle, CUBLAS_OP_N, N, N, &alpha, thrust::raw_pointer_cast(dist.data()), N,
                                thrust::raw_pointer_cast(ones.data()), 1, &beta, thrust::raw_pointer_cast(sums.data()), 1));
    thrust::copy(dist.begin(), dist.end(), qij.begin());
    // qij = (1 + ||y_i - y_j||^2)^-1 / \Sum_{k != i} (1 + ||y_i - y_k||^2)^-1
    div_row_vec(qij, sums, N, N, 1.0f);

    thrust::device_vector<float> loss_(N * N);
    thrust::transform(pij.begin(), pij.end(), qij.begin(), loss_.begin(), func_kl());

    float loss = thrust::reduce(loss_.begin(), loss_.end(), 5.0f, thrust::minimum<float>());

    // printf("dist Min: %0.5f \n", thrust::reduce(dist.begin(), dist.end(), 5000.0f, thrust::minimum<float>()));
    // printf("dist Max: %0.5f \n", thrust::reduce(dist.begin(), dist.end(), 0.0f, thrust::maximum<float>()));

    // printf("Qij Min: %0.5f \n", thrust::reduce(qij.begin(), qij.end(), 5000.0f, thrust::minimum<float>()));
    // printf("Qij Max: %0.5f \n", thrust::reduce(qij.begin(), qij.end(), 0.0f, thrust::maximum<float>()));

    // printf("pij Min: %0.5f \n", thrust::reduce(pij.begin(), pij.end(), 5000.0f, thrust::minimum<float>()));
    // printf("pij Max: %0.5f \n", thrust::reduce(pij.begin(), pij.end(), 0.0f, thrust::maximum<float>()));

    // qij = pij - qij
    thrust::transform(pij.begin(), pij.end(), qij.begin(), qij.begin(), thrust::minus<float>());

    // printf("Qij Min: %0.5f \n", thrust::reduce(qij.begin(), qij.end(), 5000.0f, thrust::minimum<float>()));
    // printf("Qij Max: %0.5f \n", thrust::reduce(qij.begin(), qij.end(), 0.0f, thrust::maximum<float>()));
    
    // qij = (pij - qij) .* (1 + ||y_i - y_j||^2)^-1
    thrust::transform(qij.begin(), qij.end(), dist.begin(), qij.begin(), thrust::multiplies<float>());

    // printf("Qij Min: %0.5f \n", thrust::reduce(qij.begin(), qij.end(), 5000.0f, thrust::minimum<float>()));
    // printf("Qij Max: %0.5f \n", thrust::reduce(qij.begin(), qij.end(), 0.0f, thrust::maximum<float>()));

    cublasSafeCall(cublasSgemv(handle, CUBLAS_OP_N, N, N, &alpha, thrust::raw_pointer_cast(qij.data()), N,
                                thrust::raw_pointer_cast(ones.data()), 1, &beta, thrust::raw_pointer_cast(forces.data()), 1));

    // printf("Qij Min: %0.5f \n", thrust::reduce(qij.begin(), qij.end(), 5000.0f, thrust::minimum<float>()));
    // printf("Qij Max: %0.5f \n", thrust::reduce(qij.begin(), qij.end(), 0.0f, thrust::maximum<float>()));

    // TODO: needs to change for 3 dimensions
    thrust::copy(forces.begin(), forces.begin() + N, forces.begin() + N);

    // forces = A * ones(N, 1) .* ys
    thrust::transform(forces.begin(), forces.end(), ys.begin(), forces.begin(), thrust::multiplies<float>());

    // printf("forces Min: %0.5f \n", thrust::reduce(forces.begin(), forces.end(), 5000.0f, thrust::minimum<float>()));
    // printf("forces Max: %0.5f \n", thrust::reduce(forces.begin(), forces.end(), 0.0f, thrust::maximum<float>()));

    alpha = -4.0f * eta;
    beta = 4.0f * eta;

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    printf("Starting TSNE calculation with %u points.\n", N);
    cudaEventRecord(start);
    do_tsne(handle, d_X, N, NDIMS);
    cudaEventRecord(stop);
