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

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

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

    
    return loss;
}


    auto stddev = sqrt(norm_sum_of_squares, N * NDIMS);

    // normalize
    broadcast_matrix_vector(points, stddev, N, NDIMS, thrust::divides<float>(), 1, 1.f);
}

thrust::device_vector<float> square(const thrust::device_vector<float> &vec, const unsigned int N) {
