 }
 
 int run_tests(int argc, char **argv) {
-    test_func(add_row_vec);
-    test_func(mul_row_vec);
-    test_func(div_row_vec);
+    test_func(thrust::plus<float>());
+    test_func(thrust::multiplies<float>());
+    test_func(thrust::divides<float>());
     return 0;
 }
+

    }
}

// Performs the operation matrix[i,:] = matrix[i,:] - alpha*vector for each row 0 <= i < N
// This is in place addition
__global__ void _add_row_vec(float * __restrict__ matrix, const float * __restrict__ vector, const unsigned int N, const unsigned int M, const float alpha) {
    const unsigned int TID = threadIdx.x + blockIdx.x * blockDim.x;
    const unsigned int i = TID % N;
    const unsigned int j = TID / N;

    if (j < M) matrix[j * N + i] = matrix[j * N + i] + alpha*vector[j];
}

void add_row_vec(thrust::device_vector<float> &matrix, thrust::device_vector<float> &vector, const unsigned int N, const unsigned int M, const float alpha) {
    const unsigned int BLOCKSIZE = 32;
    const unsigned int NBLOCKS = iDivUp(N * M, BLOCKSIZE);
    _add_row_vec<<<NBLOCKS,BLOCKSIZE>>>(thrust::raw_pointer_cast(matrix.data()), thrust::raw_pointer_cast(vector.data()), N, M, alpha);
}

// Performs the operation matrix[i,:] = alpha*matrix[i,:]*vector for each row 0 <= i < N
__global__ void _mul_row_vec(float * __restrict__ matrix, const float * __restrict__ vector, const unsigned int N, const unsigned int M, const float alpha) {
    const unsigned int TID = threadIdx.x + blockIdx.x * blockDim.x;
    const unsigned int i = TID % N;
    const unsigned int j = TID / N;

    if (j < M) matrix[j * N + i] = alpha*matrix[j * N + i]*vector[j];
}

void mul_row_vec(thrust::device_vector<float> &matrix, thrust::device_vector<float> &vector, const unsigned int N, const unsigned int M, const float alpha) {
    const unsigned int BLOCKSIZE = 32;
    const unsigned int NBLOCKS = iDivUp(N * M, BLOCKSIZE);
    std::cout << "array should be less than " << BLOCKSIZE * NBLOCKS << std::endl;
    _mul_row_vec<<<NBLOCKS,BLOCKSIZE>>>(thrust::raw_pointer_cast(matrix.data()), thrust::raw_pointer_cast(vector.data()), N, M, alpha);
}

// Performs the operation matrix[i,:] = alpha*matrix[i,:]/vector for each row 0 <= i < N
__global__ void _div_row_vec(float * __restrict__ matrix, const float * __restrict__ vector, const unsigned int N, const unsigned int M, const float alpha) {
    const unsigned int TID = threadIdx.x + blockIdx.x * blockDim.x;
    const unsigned int i = TID % N;
    const unsigned int j = TID / N;
    if (j < M) matrix[j * N + i] = alpha*matrix[j * N + i]/vector[j];
}

void div_row_vec(thrust::device_vector<float> &matrix, thrust::device_vector<float> &vector, const unsigned int N, const unsigned int M, const float alpha) {
    const unsigned int BLOCKSIZE = 32;
    const unsigned int NBLOCKS = iDivUp(N * M, BLOCKSIZE);
    _div_row_vec<<<NBLOCKS,BLOCKSIZE>>>(thrust::raw_pointer_cast(matrix.data()), thrust::raw_pointer_cast(vector.data()), N, M, alpha);
}

// Code from https://github.com/OrangeOwlSolutions/cuBLAS/blob/master/All_pairs_distances.cu
// Expects N x NDIMS matrix in points
void pairwise_dist(cublasHandle_t &handle, thrust::device_vector<float> &distances, const thrust::device_vector<float> &points, const unsigned int N, const unsigned int NDIMS) {

    __host__ __device__ double operator()(const float &x) const { return exp(x); }
};

struct func_exp_no_zero {
    __host__ __device__ double operator()(const float &x) const { return x < -1e-4 ? exp(x) : 0; }
};

struct func_inv {
    __host__ __device__ double operator()(const float &x) const { return pow(x, -1.0); }
};
