     cudaMalloc((void**) &d_fe
              , sizeof(float) * n_rows);
     for (unsigned int i=0; i < N_STREAMS_NH; ++i) {
-      cudaMalloc((void**) &d_nh[i];
+      cudaMalloc((void**) &d_nh[i]
                , sizeof(float) * n_rows * 2);
-      cudaMalloc((void**) &d_nhhd[i];
+      cudaMalloc((void**) &d_nhhd[i]
                , sizeof(float) * n_rows * 2);
-      cudaMemset(d_nh
+      cudaMemset(&d_nh[i]
                , 0
                , sizeof(float) * n_rows * 2);
-      cudaMemset(d_nhhd
+      cudaMemset(&d_nhhd[i]
                , 0
                , sizeof(float) * n_rows * 2);
     }

                  , const std::size_t n_rows
                  , const std::size_t n_cols
                  , const std::vector<float>& free_energy) {


    //TODO separate work for multiple GPUs



  }

}}} // end Clustering::Density::CUDA

                        , float* nhhd_dist_ndx
                        , unsigned int i_from
                        , unsigned int i_to) {
    extern __shared__ float smem[];
    unsigned int bid = blockIdx.x;
    unsigned int tid = threadIdx.x;
    unsigned int bsize = blockDim.x;
    unsigned int gid = bid * bsize + tid;
    if (gid + offset < n_rows) {
      // load frames for comparison into shared memory
      for (unsigned int j=0; j < n_cols; ++j) {
        smem[tid*n_cols+j] = coords[(gid+offset)*n_cols+j];
      }
