performance/graph/symmetric_rcm.cu

#include <cusp/csr_matrix.h>
#include <cusp/print.h>

#include <cusp/gallery/poisson.h>
#include <cusp/graph/symmetric_rcm.h>
#include <cusp/io/matrix_market.h>

#include <thrust/functional.h>
#include "../timer.h"

template<typename MatrixType>
size_t bandwidth(const MatrixType& G)
{
    typedef typename MatrixType::index_type IndexType;
    typedef typename MatrixType::value_type ValueType;
    typedef typename MatrixType::memory_space MemorySpace;

    cusp::coo_matrix<IndexType,ValueType,MemorySpace> G_coo(G);

    cusp::array1d<IndexType, MemorySpace> min_column(G.num_rows, 0);
    cusp::array1d<IndexType, MemorySpace> max_column(G.num_rows, 0);
    cusp::array1d<IndexType, MemorySpace> rowwise_bandwidth(G.num_rows, 0);

    thrust::reduce_by_key(G_coo.row_indices.begin(),
                          G_coo.row_indices.end(),
                          G_coo.column_indices.begin(),
                          thrust::make_discard_iterator(),
                          min_column.begin(),
                          thrust::equal_to<IndexType>(),
                          thrust::minimum<IndexType>());
    thrust::reduce_by_key(G_coo.row_indices.begin(),
                          G_coo.row_indices.end(),
                          G_coo.column_indices.begin(),
                          thrust::make_discard_iterator(),
                          max_column.begin(),
                          thrust::equal_to<IndexType>(),
                          thrust::maximum<IndexType>());

    thrust::transform(max_column.begin(), max_column.end(), min_column.begin(), rowwise_bandwidth.begin(), thrust::minus<IndexType>());
    return *thrust::max_element(rowwise_bandwidth.begin(), rowwise_bandwidth.end()) + 1;
}

template<typename MemorySpace, typename MatrixType>
void RCM(const MatrixType& G)
{
    typedef typename MatrixType::index_type IndexType;
    typedef cusp::csr_matrix<IndexType,IndexType,MemorySpace> GraphType;
    typedef cusp::array1d<IndexType,MemorySpace> Array;

    GraphType G_rcm(G);
    Array permutation(G.num_rows);

    timer t;
    cusp::graph::symmetric_rcm(G_rcm, permutation);
    std::cout << " RCM time : " << t.milliseconds_elapsed() << " (ms)." << std::endl;
    std::cout << " Bandwidth after RCM : " << bandwidth(G_rcm) << std::endl;
}

int main(int argc, char*argv[])
{
    srand(time(NULL));

    typedef int   IndexType;
    typedef float ValueType;
    typedef cusp::device_memory MemorySpace;

    cusp::coo_matrix<IndexType, ValueType, MemorySpace> A;
    size_t size = 1024;

    if (argc == 1)
    {
        // no input file was specified, generate an example
        std::cout << "Generated matrix (poisson5pt) ";
        cusp::gallery::poisson5pt(A, size, size);
    }
    else if (argc == 2)
    {
        // an input file was specified, read it from disk
        cusp::io::read_matrix_market_file(A, argv[1]);
        std::cout << "Read matrix (" << argv[1] << ") ";
    }

    std::cout << "with shape ("  << A.num_rows << "," << A.num_cols << ") and "
              << A.num_entries << " entries" << "\n\n";

    std::cout << "Bandwidth before RCM : " << bandwidth(A) << std::endl;

    std::cout << " Device ";
    RCM<cusp::device_memory>(A);

    std::cout << " Host ";
    RCM<cusp::host_memory>(A);

    return EXIT_SUCCESS;
}

1	#include <cusp/csr_matrix.h>
2	#include <cusp/print.h>
3
4	#include <cusp/gallery/poisson.h>
5	#include <cusp/graph/symmetric_rcm.h>
6	#include <cusp/io/matrix_market.h>
7
8	#include <thrust/functional.h>
9	#include "../timer.h"
10
11	template<typename MatrixType>
12	size_t bandwidth(const MatrixType& G)
13	{
14	typedef typename MatrixType::index_type IndexType;
15	typedef typename MatrixType::value_type ValueType;
16	typedef typename MatrixType::memory_space MemorySpace;
17
18	cusp::coo_matrix<IndexType,ValueType,MemorySpace> G_coo(G);
19
20	cusp::array1d<IndexType, MemorySpace> min_column(G.num_rows, 0);
21	cusp::array1d<IndexType, MemorySpace> max_column(G.num_rows, 0);
22	cusp::array1d<IndexType, MemorySpace> rowwise_bandwidth(G.num_rows, 0);
23
24	thrust::reduce_by_key(G_coo.row_indices.begin(),
25	G_coo.row_indices.end(),
26	G_coo.column_indices.begin(),
27	thrust::make_discard_iterator(),
28	min_column.begin(),
29	thrust::equal_to<IndexType>(),
30	thrust::minimum<IndexType>());
31	thrust::reduce_by_key(G_coo.row_indices.begin(),
32	G_coo.row_indices.end(),
33	G_coo.column_indices.begin(),
34	thrust::make_discard_iterator(),
35	max_column.begin(),
36	thrust::equal_to<IndexType>(),
37	thrust::maximum<IndexType>());
38
39	thrust::transform(max_column.begin(), max_column.end(), min_column.begin(), rowwise_bandwidth.begin(), thrust::minus<IndexType>());
40	return *thrust::max_element(rowwise_bandwidth.begin(), rowwise_bandwidth.end()) + 1;
41	}
42
43	template<typename MemorySpace, typename MatrixType>
44	void RCM(const MatrixType& G)
45	{
46	typedef typename MatrixType::index_type IndexType;
47	typedef cusp::csr_matrix<IndexType,IndexType,MemorySpace> GraphType;
48	typedef cusp::array1d<IndexType,MemorySpace> Array;
49
50	GraphType G_rcm(G);
51	Array permutation(G.num_rows);
52
53	timer t;
54	cusp::graph::symmetric_rcm(G_rcm, permutation);
55	std::cout << " RCM time : " << t.milliseconds_elapsed() << " (ms)." << std::endl;
56	std::cout << " Bandwidth after RCM : " << bandwidth(G_rcm) << std::endl;
57	}
58
59	int main(int argc, char*argv[])
60	{
61	srand(time(NULL));
62
63	typedef int IndexType;
64	typedef float ValueType;
65	typedef cusp::device_memory MemorySpace;
66
67	cusp::coo_matrix<IndexType, ValueType, MemorySpace> A;
68	size_t size = 1024;
69
70	if (argc == 1)
71	{
72	// no input file was specified, generate an example
73	std::cout << "Generated matrix (poisson5pt) ";
74	cusp::gallery::poisson5pt(A, size, size);
75	}
76	else if (argc == 2)
77	{
78	// an input file was specified, read it from disk
79	cusp::io::read_matrix_market_file(A, argv[1]);
80	std::cout << "Read matrix (" << argv[1] << ") ";
81	}
82
83	std::cout << "with shape (" << A.num_rows << "," << A.num_cols << ") and "
84	<< A.num_entries << " entries" << "\n\n";
85
86	std::cout << "Bandwidth before RCM : " << bandwidth(A) << std::endl;
87
88	std::cout << " Device ";
89	RCM<cusp::device_memory>(A);
90
91	std::cout << " Host ";
92	RCM<cusp::host_memory>(A);
93
94	return EXIT_SUCCESS;
95	}
96