1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
/*
* Copyright 2008-2009 NVIDIA Corporation
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#pragma once
#include <cusp/detail/device/arch.h>
#include <cusp/detail/device/common.h>
#include <cusp/detail/device/utils.h>
#include <cusp/detail/device/texture.h>
#include <thrust/device_ptr.h>
namespace cusp
{
namespace detail
{
namespace device
{
//////////////////////////////////////////////////////////////////////////////
// CSR SpMV kernels based on a vector model (one warp per row)
//////////////////////////////////////////////////////////////////////////////
//
// spmv_csr_vector_device
// Each row of the CSR matrix is assigned to a warp. The warp computes
// y[i] = A[i,:] * x, i.e. the dot product of the i-th row of A with
// the x vector, in parallel. This division of work implies that
// the CSR index and data arrays (Aj and Ax) are accessed in a contiguous
// manner (but generally not aligned). On GT200 these accesses are
// coalesced, unlike kernels based on the one-row-per-thread division of
// work. Since an entire 32-thread warp is assigned to each row, many
// threads will remain idle when their row contains a small number
// of elements. This code relies on implicit synchronization among
// threads in a warp.
//
// spmv_csr_vector_tex_device
// Same as spmv_csr_vector_tex_device, except that the texture cache is
// used for accessing the x vector.
//
// Note: THREADS_PER_VECTOR must be one of [2,4,8,16,32]
template <typename IndexType, typename ValueType, unsigned int VECTORS_PER_BLOCK, unsigned int THREADS_PER_VECTOR, bool UseCache>
__launch_bounds__(VECTORS_PER_BLOCK * THREADS_PER_VECTOR,1)
__global__ void
spmv_csr_vector_kernel(const IndexType num_rows,
const IndexType * Ap,
const IndexType * Aj,
const ValueType * Ax,
const ValueType * x,
ValueType * y)
{
__shared__ volatile ValueType sdata[VECTORS_PER_BLOCK * THREADS_PER_VECTOR + THREADS_PER_VECTOR / 2]; // padded to avoid reduction conditionals
__shared__ volatile IndexType ptrs[VECTORS_PER_BLOCK][2];
const IndexType THREADS_PER_BLOCK = VECTORS_PER_BLOCK * THREADS_PER_VECTOR;
const IndexType thread_id = THREADS_PER_BLOCK * blockIdx.x + threadIdx.x; // global thread index
const IndexType thread_lane = threadIdx.x & (THREADS_PER_VECTOR - 1); // thread index within the vector
const IndexType vector_id = thread_id / THREADS_PER_VECTOR; // global vector index
const IndexType vector_lane = threadIdx.x / THREADS_PER_VECTOR; // vector index within the block
const IndexType num_vectors = VECTORS_PER_BLOCK * gridDim.x; // total number of active vectors
for(IndexType row = vector_id; row < num_rows; row += num_vectors)
{
// use two threads to fetch Ap[row] and Ap[row+1]
// this is considerably faster than the straightforward version
if(thread_lane < 2)
ptrs[vector_lane][thread_lane] = Ap[row + thread_lane];
const IndexType row_start = ptrs[vector_lane][0]; //same as: row_start = Ap[row];
const IndexType row_end = ptrs[vector_lane][1]; //same as: row_end = Ap[row+1];
// initialize local sum
ValueType sum = 0;
if (THREADS_PER_VECTOR == 32 && row_end - row_start > 32)
{
// ensure aligned memory access to Aj and Ax
IndexType jj = row_start - (row_start & (THREADS_PER_VECTOR - 1)) + thread_lane;
// accumulate local sums
if(jj >= row_start && jj < row_end)
sum += Ax[jj] * fetch_x<UseCache>(Aj[jj], x);
// accumulate local sums
for(jj += THREADS_PER_VECTOR; jj < row_end; jj += THREADS_PER_VECTOR)
sum += Ax[jj] * fetch_x<UseCache>(Aj[jj], x);
}
else
{
// accumulate local sums
for(IndexType jj = row_start + thread_lane; jj < row_end; jj += THREADS_PER_VECTOR)
sum += Ax[jj] * fetch_x<UseCache>(Aj[jj], x);
}
// store local sum in shared memory
sdata[threadIdx.x] = sum;
// reduce local sums to row sum
if (THREADS_PER_VECTOR > 16) sdata[threadIdx.x] = sum = sum + sdata[threadIdx.x + 16];
if (THREADS_PER_VECTOR > 8) sdata[threadIdx.x] = sum = sum + sdata[threadIdx.x + 8];
if (THREADS_PER_VECTOR > 4) sdata[threadIdx.x] = sum = sum + sdata[threadIdx.x + 4];
if (THREADS_PER_VECTOR > 2) sdata[threadIdx.x] = sum = sum + sdata[threadIdx.x + 2];
if (THREADS_PER_VECTOR > 1) sdata[threadIdx.x] = sum = sum + sdata[threadIdx.x + 1];
// first thread writes the result
if (thread_lane == 0)
y[row] = sdata[threadIdx.x];
}
}
template <bool UseCache, unsigned int THREADS_PER_VECTOR, typename Matrix, typename ValueType>
void __spmv_csr_vector(const Matrix& A,
const ValueType* x,
ValueType* y)
{
typedef typename Matrix::index_type IndexType;
const size_t THREADS_PER_BLOCK = 128;
const size_t VECTORS_PER_BLOCK = THREADS_PER_BLOCK / THREADS_PER_VECTOR;
const size_t MAX_BLOCKS = cusp::detail::device::arch::max_active_blocks(spmv_csr_vector_kernel<IndexType, ValueType, VECTORS_PER_BLOCK, THREADS_PER_VECTOR, UseCache>, THREADS_PER_BLOCK, (size_t) 0);
const size_t NUM_BLOCKS = std::min<size_t>(MAX_BLOCKS, DIVIDE_INTO(A.num_rows, VECTORS_PER_BLOCK));
if (UseCache)
bind_x(x);
spmv_csr_vector_kernel<IndexType, ValueType, VECTORS_PER_BLOCK, THREADS_PER_VECTOR, UseCache> <<<NUM_BLOCKS, THREADS_PER_BLOCK>>>
(A.num_rows,
thrust::raw_pointer_cast(&A.row_offsets[0]),
thrust::raw_pointer_cast(&A.column_indices[0]),
thrust::raw_pointer_cast(&A.values[0]),
x, y);
if (UseCache)
unbind_x(x);
}
template <typename Matrix,
typename ValueType>
void spmv_csr_vector(const Matrix& A,
const ValueType* x,
ValueType* y)
{
typedef typename Matrix::index_type IndexType;
const IndexType nnz_per_row = A.num_entries / A.num_rows;
if (nnz_per_row <= 2) { __spmv_csr_vector<false, 2>(A, x, y); return; }
if (nnz_per_row <= 4) { __spmv_csr_vector<false, 4>(A, x, y); return; }
if (nnz_per_row <= 8) { __spmv_csr_vector<false, 8>(A, x, y); return; }
if (nnz_per_row <= 16) { __spmv_csr_vector<false,16>(A, x, y); return; }
__spmv_csr_vector<false,32>(A, x, y);
}
template <typename Matrix,
typename ValueType>
void spmv_csr_vector_tex(const Matrix& A,
const ValueType* x,
ValueType* y)
{
typedef typename Matrix::index_type IndexType;
const IndexType nnz_per_row = A.num_entries / A.num_rows;
if (nnz_per_row <= 2) { __spmv_csr_vector<true, 2>(A, x, y); return; }
if (nnz_per_row <= 4) { __spmv_csr_vector<true, 4>(A, x, y); return; }
if (nnz_per_row <= 8) { __spmv_csr_vector<true, 8>(A, x, y); return; }
if (nnz_per_row <= 16) { __spmv_csr_vector<true,16>(A, x, y); return; }
__spmv_csr_vector<true,32>(A, x, y);
}
} // end namespace device
} // end namespace detail
} // end namespace cusp
|
