 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);
    }

    CsrMatrix<ValueT, OffsetT> csr_matrix;
    csr_matrix.FromCoo(coo_matrix);

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

    // Adaptive timing iterations: run two billion nonzeros through
    if (timing_iterations == -1)
    {
        timing_iterations = std::max(5, std::min(OffsetT(1e4), (OffsetT(2e9) / csr_matrix.num_nonzeros)));
        if (!g_quiet)
            printf("\t%d timing iterations\n", timing_iterations);
    }

    // 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;
}



    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
