     fflush(stdout);
 }
 
+
+
+/**
+ * Run tests
+ */
+template <
+    typename ValueT,
+    typename OffsetT>
+void RunTest(
+    bool                        rcm_relabel,
+    ValueT                      alpha,
+    ValueT                      beta,
+    CooMatrix<ValueT, OffsetT>& coo_matrix,
+    int                         timing_iterations,
+    CommandLineArgs&            args)
+{
+    // Adaptive timing iterations: run 16 billion nonzeros through
+    if (timing_iterations == -1)
+    {
+        timing_iterations = std::min(50000ull, std::max(100ull, ((16ull << 30) / coo_matrix.num_nonzeros)));
+        if (!g_quiet)
+            printf("\t%d timing iterations\n", timing_iterations);
+    }
+
+    CsrMatrix<ValueT, OffsetT> csr_matrix;
+    csr_matrix.FromCoo(coo_matrix);
+    coo_matrix.Clear();
+
+    // Relabel
+    if (rcm_relabel)
+    {
+        if (!g_quiet)
+        {
+            csr_matrix.Stats().Display();
+            printf("\n");
+            csr_matrix.DisplayHistogram();
+            printf("\n");
+            if (g_verbose2)
+                csr_matrix.Display();
+            printf("\n");
+        }
+
+        RcmRelabel(csr_matrix, !g_quiet);
+
+        if (!g_quiet)
+        printf("\n");
+    }
+
+    // Display matrix info
+    csr_matrix.Stats().Display(!g_quiet);
+    if (!g_quiet)
+    {
+        printf("\n");
+        csr_matrix.DisplayHistogram();
+        printf("\n");
+        if (g_verbose2)
+            csr_matrix.Display();
+        printf("\n");
+    }
+    fflush(stdout);
+
+    // Allocate input and output vectors
+    ValueT* vector_x        = new ValueT[csr_matrix.num_cols];
+    ValueT* vector_y_in     = new ValueT[csr_matrix.num_rows];
+    ValueT* vector_y_out    = new ValueT[csr_matrix.num_rows];
+
+    for (int col = 0; col < csr_matrix.num_cols; ++col)
+        vector_x[col] = 1.0;
+
+    for (int row = 0; row < csr_matrix.num_rows; ++row)
+        vector_y_in[row] = 1.0;
+
+    // Compute reference answer
+    SpmvGold(csr_matrix, vector_x, vector_y_in, vector_y_out, alpha, beta);
+
+    float avg_millis;
+
+    {
+        if (g_quiet) {
+            printf("%s, %s, ", args.deviceProp.name, (sizeof(ValueT) > 4) ? "fp64" : "fp32"); fflush(stdout);
+        }
+
+        // Get GPU device bandwidth (GB/s)
+        float device_giga_bandwidth = args.device_giga_bandwidth;
+
+        if (((size_t(csr_matrix.num_rows) * (sizeof(OffsetT) + (sizeof(ValueT) * 2))) +
+            (size_t(csr_matrix.num_nonzeros) * (sizeof(OffsetT) + sizeof(ValueT)))) > args.device_free_physmem)
+        {
+            // Won't fit
+            printf("Too big\n"); fflush(stdout);
+        }
+        else
+        {
+            // Will fit
+
+            // Allocate and initialize GPU problem
+            SpmvParams<ValueT, OffsetT> params;
+
+            g_allocator.DeviceAllocate((void **) &params.d_values,          sizeof(ValueT) * csr_matrix.num_nonzeros);
+            g_allocator.DeviceAllocate((void **) &params.d_row_end_offsets, sizeof(OffsetT) * (csr_matrix.num_rows + 1));
+            g_allocator.DeviceAllocate((void **) &params.d_column_indices,  sizeof(OffsetT) * csr_matrix.num_nonzeros);
+            g_allocator.DeviceAllocate((void **) &params.d_vector_x,        sizeof(ValueT) * csr_matrix.num_cols);
+            g_allocator.DeviceAllocate((void **) &params.d_vector_y,        sizeof(ValueT) * csr_matrix.num_rows);
+            params.num_rows = csr_matrix.num_rows;
+            params.num_cols = csr_matrix.num_cols;
+            params.num_nonzeros = csr_matrix.num_nonzeros;
+            params.alpha = alpha;
+            params.beta = beta;
+
+            CubDebugExit(cudaMemcpy(params.d_values,            csr_matrix.values,          sizeof(ValueT) * csr_matrix.num_nonzeros, cudaMemcpyHostToDevice));
+            CubDebugExit(cudaMemcpy(params.d_row_end_offsets,   csr_matrix.row_offsets,     sizeof(OffsetT) * (csr_matrix.num_rows + 1), cudaMemcpyHostToDevice));
+            CubDebugExit(cudaMemcpy(params.d_column_indices,    csr_matrix.column_indices,  sizeof(OffsetT) * csr_matrix.num_nonzeros, cudaMemcpyHostToDevice));
+            CubDebugExit(cudaMemcpy(params.d_vector_x,          vector_x,                   sizeof(ValueT) * csr_matrix.num_cols, cudaMemcpyHostToDevice));
+
+            if (!g_quiet) {
+                printf("\n\nGPU CSR I/O Prox: "); fflush(stdout);
+            }
+            avg_millis = TestGpuCsrIoProxy(params, timing_iterations);
+            DisplayPerf(device_giga_bandwidth, avg_millis, csr_matrix);
+
+            if (!g_quiet) {
+                printf("\n\nCUB: "); fflush(stdout);
+            }
+            avg_millis = TestGpuMergeCsrmv(vector_y_in, vector_y_out, params, timing_iterations);
+            DisplayPerf(device_giga_bandwidth, avg_millis, csr_matrix);
+
+            // Initalize cuSparse
+            cusparseHandle_t cusparse;
+            AssertEquals(CUSPARSE_STATUS_SUCCESS, cusparseCreate(&cusparse));
+
+            if (!g_quiet) {
+                printf("\n\nCusparse: "); fflush(stdout);
+            }
+            avg_millis = TestCusparseCsrmv(vector_y_in, vector_y_out, params, timing_iterations, cusparse);
+            DisplayPerf(device_giga_bandwidth, avg_millis, csr_matrix);
+
+            // Cleanup
+            if (params.d_values)            g_allocator.DeviceFree(params.d_values);
+            if (params.d_row_end_offsets)   g_allocator.DeviceFree(params.d_row_end_offsets);
+            if (params.d_column_indices)    g_allocator.DeviceFree(params.d_column_indices);
+            if (params.d_vector_x)          g_allocator.DeviceFree(params.d_vector_x);
+            if (params.d_vector_y)          g_allocator.DeviceFree(params.d_vector_y);
+        }
+    }
+
+    if (vector_x)                   delete[] vector_x;
+    if (vector_y_in)                delete[] vector_y_in;
+    if (vector_y_out)               delete[] vector_y_out;
+}
 
 /**
  * Run tests

bool                    g_verbose   = false;        // Whether to display output to console
bool                    g_verbose2  = false;        // Whether to display input to console
CachingDeviceAllocator  g_allocator(true);          // Caching allocator for device memory
int                     g_omp_threads = -1;         // Number of openMP threads


//---------------------------------------------------------------------

        exit(1);
    }

    RunTest(
        rcm_relabel,
        alpha,
        beta,
        coo_matrix,
        timing_iterations,
        args);
}


