 	cublasSafeCall(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_T, N, N, NDIMS, &alpha,
 		                       thrust::raw_pointer_cast(points.data()), N, thrust::raw_pointer_cast(points.data()), N, &beta,
 							   thrust::raw_pointer_cast(distances.data()), N));
- 
 	dim3 dimBlock(BLOCKSIZE, BLOCKSIZE);
 	dim3 dimGrid(iDivUp(N, BLOCKSIZE), iDivUp(N, BLOCKSIZE));
 	assemble_final_result<<<dimGrid, dimBlock>>>(thrust::raw_pointer_cast(squared_norms.data()), 

    auto means = reduce_mean(handle, points, N, NDIMS, 0);

    // zero center
    add_row_vec(points, means, N, NDIMS, -1.f);
    
    // compute standard deviation
    auto squared_vals = square(points, N * NDIMS);

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