Actual source code: sellcuda.cu
1: #include <petscdevice_cuda.h>
2: #include <petsc/private/cupmatomics.hpp>
3: #include <../src/mat/impls/sell/seq/sell.h>
5: #define SLICE_HEIGHT 16
7: typedef struct {
8: PetscInt maxallocmat;
9: PetscInt totalentries;
10: PetscInt *colidx; /* column index array, device pointer */
11: MatScalar *val; /* value array, device pointer */
12: PetscInt totalslices;
13: PetscInt *sliidx; /* slice index array, device pointer */
14: PetscObjectState nonzerostate;
15: PetscInt kernelchoice;
16: PetscInt blocky;
17: PetscInt chunksperblock;
18: PetscInt totalchunks;
19: PetscInt *chunk_slice_map; /* starting slice for each chunk, device pointer */
20: } Mat_SeqSELLCUDA;
22: static PetscErrorCode MatSeqSELLCUDA_Destroy(Mat_SeqSELLCUDA **cudastruct)
23: {
24: PetscFunctionBegin;
25: if (*cudastruct) {
26: PetscCallCUDA(cudaFree((*cudastruct)->colidx));
27: PetscCallCUDA(cudaFree((*cudastruct)->val));
28: PetscCallCUDA(cudaFree((*cudastruct)->sliidx));
29: PetscCallCUDA(cudaFree((*cudastruct)->chunk_slice_map));
30: PetscCall(PetscFree(*cudastruct));
31: }
32: PetscFunctionReturn(PETSC_SUCCESS);
33: }
35: static PetscErrorCode MatSeqSELLCUDACopyToGPU(Mat A)
36: {
37: Mat_SeqSELLCUDA *cudastruct = (Mat_SeqSELLCUDA *)A->spptr;
38: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
40: PetscFunctionBegin;
41: if (A->offloadmask == PETSC_OFFLOAD_UNALLOCATED || A->offloadmask == PETSC_OFFLOAD_CPU) {
42: PetscCall(PetscLogEventBegin(MAT_CUDACopyToGPU, A, 0, 0, 0));
43: if (cudastruct->sliidx && A->nonzerostate == cudastruct->nonzerostate) {
44: /* copy values only */
45: PetscCallCUDA(cudaMemcpy(cudastruct->val, a->val, a->sliidx[a->totalslices] * sizeof(MatScalar), cudaMemcpyHostToDevice));
46: PetscCall(PetscLogCpuToGpu(a->sliidx[a->totalslices] * (sizeof(MatScalar))));
47: } else {
48: PetscCallCUDA(cudaFree(cudastruct->colidx));
49: PetscCallCUDA(cudaFree(cudastruct->val));
50: PetscCallCUDA(cudaFree(cudastruct->sliidx));
51: PetscCallCUDA(cudaFree(cudastruct->chunk_slice_map));
52: cudastruct->maxallocmat = a->maxallocmat;
53: cudastruct->totalentries = a->sliidx[a->totalslices];
54: cudastruct->totalslices = a->totalslices;
55: cudastruct->totalchunks = a->totalchunks;
56: PetscCallCUDA(cudaMalloc((void **)&cudastruct->colidx, a->maxallocmat * sizeof(*cudastruct->colidx)));
57: PetscCallCUDA(cudaMalloc((void **)&cudastruct->val, a->maxallocmat * sizeof(*cudastruct->val)));
58: /* copy values, nz or maxallocmat? */
59: PetscCallCUDA(cudaMemcpy(cudastruct->colidx, a->colidx, a->sliidx[a->totalslices] * sizeof(*a->colidx), cudaMemcpyHostToDevice));
60: PetscCallCUDA(cudaMemcpy(cudastruct->val, a->val, a->sliidx[a->totalslices] * sizeof(*a->val), cudaMemcpyHostToDevice));
62: PetscCallCUDA(cudaMalloc((void **)&cudastruct->sliidx, (a->totalslices + 1) * sizeof(*cudastruct->sliidx)));
63: PetscCallCUDA(cudaMemcpy(cudastruct->sliidx, a->sliidx, (a->totalslices + 1) * sizeof(*a->sliidx), cudaMemcpyHostToDevice));
64: PetscCallCUDA(cudaMalloc((void **)&cudastruct->chunk_slice_map, a->totalchunks * sizeof(*cudastruct->chunk_slice_map)));
65: PetscCallCUDA(cudaMemcpy(cudastruct->chunk_slice_map, a->chunk_slice_map, a->totalchunks * sizeof(*a->chunk_slice_map), cudaMemcpyHostToDevice));
66: PetscCall(PetscLogCpuToGpu(a->sliidx[a->totalslices] * (sizeof(MatScalar) + sizeof(PetscInt)) + (a->totalslices + 1 + a->totalchunks) * sizeof(PetscInt)));
67: }
68: PetscCallCUDA(WaitForCUDA());
69: PetscCall(PetscLogEventEnd(MAT_CUDACopyToGPU, A, 0, 0, 0));
70: cudastruct->nonzerostate = A->nonzerostate;
71: A->offloadmask = PETSC_OFFLOAD_BOTH;
72: }
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: static __global__ void matmult_seqsell_basic_kernel(PetscInt nrows, PetscInt sliceheight, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
77: {
78: PetscInt i, row, slice_id, row_in_slice;
79: MatScalar sum;
80: /* one thread per row. */
81: row = blockIdx.x * blockDim.x + threadIdx.x;
82: if (row < nrows) {
83: slice_id = row / sliceheight;
84: row_in_slice = row % sliceheight;
85: sum = 0.0;
86: for (i = sliidx[slice_id] + row_in_slice; i < sliidx[slice_id + 1]; i += sliceheight) sum += aval[i] * x[acolidx[i]];
87: y[row] = sum;
88: }
89: }
91: static __global__ void matmultadd_seqsell_basic_kernel(PetscInt nrows, PetscInt sliceheight, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
92: {
93: PetscInt i, row, slice_id, row_in_slice;
94: MatScalar sum;
95: /* one thread per row. */
96: row = blockIdx.x * blockDim.x + threadIdx.x;
97: if (row < nrows) {
98: slice_id = row / sliceheight;
99: row_in_slice = row % sliceheight;
100: sum = 0.0;
101: for (i = sliidx[slice_id] + row_in_slice; i < sliidx[slice_id + 1]; i += sliceheight) sum += aval[i] * x[acolidx[i]];
102: z[row] = y[row] + sum;
103: }
104: }
106: #if !PetscDefined(USE_COMPLEX)
107: /* use 1 block per slice, suitable for large slice width */
108: template <int BLOCKY>
109: __global__ void matmult_seqsell_tiled_kernel9(PetscInt nrows, PetscInt sliceheight, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
110: {
111: __shared__ MatScalar shared[32][BLOCKY];
112: PetscInt i, row, slice_id = blockIdx.x;
113: int tid = threadIdx.x + threadIdx.y * 32;
114: /* transposed index */
115: int tidx = tid % BLOCKY;
116: int tidy = tid / BLOCKY;
117: PetscScalar t = 0.0;
119: row = slice_id * sliceheight + threadIdx.x % sliceheight;
120: if (row < nrows) {
121: for (i = sliidx[slice_id] + threadIdx.x + 32 * threadIdx.y; i < sliidx[slice_id + 1]; i += 32 * BLOCKY) t += aval[i] * x[acolidx[i]];
122: }
123: #pragma unroll
124: for (int offset = 16; offset >= sliceheight; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset);
125: /* transpose layout to reduce each row using warp shfl */
126: if (threadIdx.x < sliceheight) shared[threadIdx.x][threadIdx.y] = t;
127: __syncthreads();
128: if (tidy < sliceheight) t = shared[tidy][tidx];
129: #pragma unroll
130: for (int offset = BLOCKY / 2; offset > 0; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset, BLOCKY);
131: if (tidx == 0 && tidy < sliceheight) shared[0][tidy] = t;
132: __syncthreads();
133: if (row < nrows && threadIdx.y == 0 && threadIdx.x < sliceheight) y[row] = shared[0][threadIdx.x];
134: }
136: /* use 1 block per slice, suitable for large slice width */
137: template <int BLOCKY>
138: __global__ void matmultadd_seqsell_tiled_kernel9(PetscInt nrows, PetscInt sliceheight, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
139: {
140: __shared__ MatScalar shared[32][BLOCKY];
141: PetscInt i, row, slice_id = blockIdx.x;
142: int tid = threadIdx.x + threadIdx.y * 32;
143: /* transposed index */
144: int tidx = tid % BLOCKY;
145: int tidy = tid / BLOCKY;
146: PetscScalar t = 0.0;
148: row = slice_id * sliceheight + threadIdx.x % sliceheight;
149: if (row < nrows) {
150: for (i = sliidx[slice_id] + threadIdx.x + 32 * threadIdx.y; i < sliidx[slice_id + 1]; i += 32 * BLOCKY) t += aval[i] * x[acolidx[i]];
151: }
152: #pragma unroll
153: for (int offset = 16; offset >= sliceheight; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset);
154: /* transpose layout to reduce each row using warp shfl */
155: if (threadIdx.x < sliceheight) shared[threadIdx.x][threadIdx.y] = t;
156: __syncthreads();
157: if (tidy < sliceheight) t = shared[tidy][tidx];
158: #pragma unroll
159: for (int offset = BLOCKY / 2; offset > 0; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset, BLOCKY);
160: if (tidx == 0 && tidy < sliceheight) shared[0][tidy] = t;
161: __syncthreads();
162: if (row < nrows && threadIdx.y == 0 && threadIdx.x < sliceheight) z[row] = y[row] + shared[0][threadIdx.x];
163: }
165: template <int BLOCKY>
166: __device__ __forceinline__ static bool segment_scan(PetscInt flag[], MatScalar shared[], PetscScalar *val)
167: {
168: bool head = true;
169: #pragma unroll
170: for (int i = 1; i < BLOCKY * 2; i <<= 1) {
171: int halfwarpid = threadIdx.y * 2 + threadIdx.x / 16;
172: shared[threadIdx.x + threadIdx.y * 32] = 0;
173: if (halfwarpid >= i && flag[halfwarpid - i] == flag[halfwarpid]) {
174: shared[threadIdx.x + threadIdx.y * 32] = *val;
175: if (i == 1) head = false;
176: }
177: __syncthreads();
178: if (halfwarpid < BLOCKY * 2 - i) *val += shared[threadIdx.x + threadIdx.y * 32 + i * 16];
179: __syncthreads();
180: }
181: return head;
182: }
184: /* load-balancing version. Chunksize is equal to the number of threads per block */
185: template <int BLOCKY>
186: __global__ void matmult_seqsell_tiled_kernel8(PetscInt nrows, PetscInt sliceheight, PetscInt chunksperblock, PetscInt totalchunks, const PetscInt *chunk_slice_map, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
187: {
188: __shared__ MatScalar shared[BLOCKY * 32];
189: PetscInt gid, row, start_slice, cid;
190: PetscScalar t = 0.0;
191: AtomicAdd<MatScalar> atomAdd;
192: /* zero out y */
193: for (int iter = 0; iter < 1 + (nrows - 1) / (gridDim.x * 32 * BLOCKY); iter++) {
194: gid = gridDim.x * 32 * BLOCKY * iter + blockIdx.x * BLOCKY * 32 + threadIdx.y * 32 + threadIdx.x;
195: if (gid < nrows) y[gid] = 0.0;
196: }
197: for (int iter = 0; iter < chunksperblock; iter++) {
198: cid = blockIdx.x * chunksperblock + iter; /* chunk id */
199: if (cid < totalchunks) {
200: start_slice = chunk_slice_map[cid]; /* starting slice at each iteration */
201: gid = cid * BLOCKY * 32 + threadIdx.y * 32 + threadIdx.x;
202: if ((cid + 1) * BLOCKY * 32 > sliidx[start_slice + 1]) { /* this iteration covers more than one slice */
203: __shared__ PetscInt flag[BLOCKY * 2];
204: bool write;
205: PetscInt slice_id = start_slice, totalslices = PetscCeilIntMacro(nrows, sliceheight), totalentries = sliidx[totalslices];
206: /* find out the slice that this element belongs to */
207: while (gid < totalentries && gid >= sliidx[slice_id + 1]) slice_id++;
208: if (threadIdx.x % 16 == 0) flag[threadIdx.y * 2 + threadIdx.x / 16] = slice_id;
209: row = slice_id * sliceheight + threadIdx.x % sliceheight;
210: if (row < nrows && gid < totalentries) t = aval[gid] * x[acolidx[gid]];
211: __syncthreads();
212: write = segment_scan<BLOCKY>(flag, shared, &t);
213: if (row < nrows && gid < totalentries && write) atomAdd(y[row], t);
214: t = 0.0;
215: } else { /* this iteration covers only one slice */
216: row = start_slice * sliceheight + threadIdx.x % sliceheight;
217: if (row < nrows) t += aval[gid] * x[acolidx[gid]];
218: if (iter == chunksperblock - 1 || (cid + 2) * BLOCKY * 32 > sliidx[start_slice + 1]) { /* last iteration or next iteration covers more than one slice */
219: int tid = threadIdx.x + threadIdx.y * 32, tidx = tid % BLOCKY, tidy = tid / BLOCKY;
220: /* reduction and write to output vector */
221: #pragma unroll
222: for (int offset = 16; offset >= sliceheight; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset);
223: /* transpose layout to reduce each row using warp shfl */
224: if (threadIdx.x < sliceheight) shared[threadIdx.x * BLOCKY + threadIdx.y] = t; /* shared[threadIdx.x][threadIdx.y] = t */
225: __syncthreads();
226: if (tidy < sliceheight) t = shared[tidy * BLOCKY + tidx]; /* shared[tidy][tidx] */
227: #pragma unroll
228: for (int offset = BLOCKY / 2; offset > 0; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset, BLOCKY);
229: if (tidx == 0 && tidy < sliceheight) shared[tidy] = t; /* shared[0][tidy] = t */
230: __syncthreads();
231: if (row < nrows && threadIdx.y == 0 && threadIdx.x < sliceheight) atomAdd(y[row], shared[threadIdx.x]); /* shared[0][threadIdx.x] */
232: t = 0.0;
233: }
234: }
235: }
236: }
237: }
239: /* load-balancing version. Chunksize is equal to the number of threads per block */
240: template <int BLOCKY>
241: __global__ void matmultadd_seqsell_tiled_kernel8(PetscInt nrows, PetscInt sliceheight, PetscInt chunksperblock, PetscInt totalchunks, const PetscInt *chunk_slice_map, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
242: {
243: __shared__ MatScalar shared[BLOCKY * 32];
244: PetscInt gid, row, start_slice, cid;
245: PetscScalar t = 0.0;
246: AtomicAdd<MatScalar> atomAdd;
247: /* copy y to z */
248: for (int iter = 0; iter < 1 + (nrows - 1) / (gridDim.x * 32 * BLOCKY); iter++) {
249: gid = gridDim.x * 32 * BLOCKY * iter + blockIdx.x * BLOCKY * 32 + threadIdx.y * 32 + threadIdx.x;
250: if (gid < nrows) z[gid] = y[gid];
251: }
252: for (int iter = 0; iter < chunksperblock; iter++) {
253: cid = blockIdx.x * chunksperblock + iter; /* chunk id */
254: if (cid < totalchunks) {
255: start_slice = chunk_slice_map[cid]; /* starting slice at each iteration */
256: gid = cid * BLOCKY * 32 + threadIdx.y * 32 + threadIdx.x;
257: if ((cid + 1) * BLOCKY * 32 > sliidx[start_slice + 1]) { /* this iteration covers more than one slice */
258: __shared__ PetscInt flag[BLOCKY * 2];
259: bool write;
260: PetscInt slice_id = start_slice, totalslices = PetscCeilIntMacro(nrows, sliceheight), totalentries = sliidx[totalslices];
261: /* find out the slice that this element belongs to */
262: while (gid < totalentries && gid >= sliidx[slice_id + 1]) slice_id++;
263: if (threadIdx.x % 16 == 0) flag[threadIdx.y * 2 + threadIdx.x / 16] = slice_id;
264: row = slice_id * sliceheight + threadIdx.x % sliceheight;
265: if (row < nrows && gid < totalentries) t = aval[gid] * x[acolidx[gid]];
266: __syncthreads();
267: write = segment_scan<BLOCKY>(flag, shared, &t);
268: if (row < nrows && gid < totalentries && write) atomAdd(z[row], t);
269: t = 0.0;
270: } else { /* this iteration covers only one slice */
271: row = start_slice * sliceheight + threadIdx.x % sliceheight;
272: if (row < nrows) t += aval[gid] * x[acolidx[gid]];
273: if (iter == chunksperblock - 1 || (cid + 2) * BLOCKY * 32 > sliidx[start_slice + 1]) { /* last iteration or next iteration covers more than one slice */
274: int tid = threadIdx.x + threadIdx.y * 32, tidx = tid % BLOCKY, tidy = tid / BLOCKY;
275: /* reduction and write to output vector */
276: #pragma unroll
277: for (int offset = 16; offset >= sliceheight; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset);
278: /* transpose layout to reduce each row using warp shfl */
279: if (threadIdx.x < sliceheight) shared[threadIdx.x * BLOCKY + threadIdx.y] = t; /* shared[threadIdx.x][threadIdx.y] = t */
280: __syncthreads();
281: if (tidy < sliceheight) t = shared[tidy * BLOCKY + tidx]; /* shared[tidy][tidx] */
282: #pragma unroll
283: for (int offset = BLOCKY / 2; offset > 0; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset, BLOCKY);
284: if (tidx == 0 && tidy < sliceheight) shared[tidy] = t; /* shared[0][tidy] = t */
285: __syncthreads();
286: if (row < nrows && threadIdx.y == 0 && threadIdx.x < sliceheight) atomAdd(z[row], shared[threadIdx.x]); /* shared[0][threadIdx.x] */
287: t = 0.0;
288: }
289: }
290: }
291: }
292: }
294: /* use 1 warp per slice, suitable for small slice width */
295: static __global__ void matmult_seqsell_tiled_kernel7(PetscInt nrows, PetscInt sliceheight, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
296: {
297: PetscInt i, row, slice_id;
298: slice_id = blockIdx.x * blockDim.y + threadIdx.y;
299: row = slice_id * sliceheight + threadIdx.x % sliceheight;
300: double t = 0.0;
301: if (row < nrows) {
302: for (i = sliidx[slice_id] + threadIdx.x; i < sliidx[slice_id + 1]; i += 32) t += aval[i] * x[acolidx[i]];
303: }
304: #pragma unroll
305: for (int offset = 16; offset >= sliceheight; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset);
306: if (row < nrows && threadIdx.x < sliceheight) y[row] = t;
307: }
309: /* use 1 warp per slice, suitable for small slice width */
310: static __global__ void matmultadd_seqsell_tiled_kernel7(PetscInt nrows, PetscInt sliceheight, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
311: {
312: PetscInt i, row, slice_id;
313: slice_id = blockIdx.x * blockDim.y + threadIdx.y;
314: row = slice_id * sliceheight + threadIdx.x % sliceheight;
315: double t = 0.0;
316: if (row < nrows) {
317: for (i = sliidx[slice_id] + threadIdx.x; i < sliidx[slice_id + 1]; i += 32) t += aval[i] * x[acolidx[i]];
318: }
319: #pragma unroll
320: for (int offset = 16; offset >= sliceheight; offset /= 2) t += __shfl_down_sync(0xffffffff, t, offset);
321: if (row < nrows && threadIdx.x < sliceheight) z[row] = y[row] + t;
322: }
323: #endif
325: /*********** Kernel 2-6 are tied to slice height 16. They are kept only for performance comparison **********/
327: static __global__ void matmult_seqsell_tiled_kernel6(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
328: {
329: __shared__ MatScalar shared[512];
330: PetscInt i, row, slice_id, row_in_slice;
331: /* multiple threads per row. */
332: row = blockIdx.x * blockDim.x + threadIdx.x;
333: if (row < nrows) {
334: slice_id = row / SLICE_HEIGHT;
335: row_in_slice = row % SLICE_HEIGHT;
337: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
338: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
339: __syncthreads();
340: if (threadIdx.y < 16) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 16) * blockDim.x + threadIdx.x];
341: __syncthreads();
342: if (threadIdx.y < 8) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 8) * blockDim.x + threadIdx.x];
343: __syncthreads();
344: if (threadIdx.y < 4) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 4) * blockDim.x + threadIdx.x];
345: __syncthreads();
346: if (threadIdx.y < 2) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 2) * blockDim.x + threadIdx.x];
347: __syncthreads();
348: if (threadIdx.y < 1) {
349: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
350: y[row] = shared[threadIdx.x];
351: }
352: }
353: }
355: static __global__ void matmult_seqsell_tiled_kernel5(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
356: {
357: __shared__ MatScalar shared[512];
358: PetscInt i, row, slice_id, row_in_slice;
359: /* multiple threads per row. */
360: row = blockIdx.x * blockDim.x + threadIdx.x;
361: if (row < nrows) {
362: slice_id = row / SLICE_HEIGHT;
363: row_in_slice = row % SLICE_HEIGHT;
365: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
366: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
367: __syncthreads();
368: if (threadIdx.y < 8) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 8) * blockDim.x + threadIdx.x];
369: __syncthreads();
370: if (threadIdx.y < 4) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 4) * blockDim.x + threadIdx.x];
371: __syncthreads();
372: if (threadIdx.y < 2) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 2) * blockDim.x + threadIdx.x];
373: __syncthreads();
374: if (threadIdx.y < 1) {
375: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
376: y[row] = shared[threadIdx.x];
377: }
378: }
379: }
381: static __global__ void matmult_seqsell_tiled_kernel4(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
382: {
383: __shared__ MatScalar shared[512];
384: PetscInt i, row, slice_id, row_in_slice;
385: /* multiple threads per row. */
386: row = blockIdx.x * blockDim.x + threadIdx.x;
387: if (row < nrows) {
388: slice_id = row / SLICE_HEIGHT;
389: row_in_slice = row % SLICE_HEIGHT;
391: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
392: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
393: __syncthreads();
394: if (threadIdx.y < 4) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 4) * blockDim.x + threadIdx.x];
395: __syncthreads();
396: if (threadIdx.y < 2) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 2) * blockDim.x + threadIdx.x];
397: __syncthreads();
398: if (threadIdx.y < 1) {
399: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
400: y[row] = shared[threadIdx.x];
401: }
402: }
403: }
405: static __global__ void matmult_seqsell_tiled_kernel3(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
406: {
407: __shared__ MatScalar shared[512];
408: PetscInt i, row, slice_id, row_in_slice;
409: /* multiple threads per row. */
410: row = blockIdx.x * blockDim.x + threadIdx.x;
411: if (row < nrows) {
412: slice_id = row / SLICE_HEIGHT;
413: row_in_slice = row % SLICE_HEIGHT;
415: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
416: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
417: __syncthreads();
418: if (threadIdx.y < 2) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 2) * blockDim.x + threadIdx.x];
419: __syncthreads();
420: if (threadIdx.y < 1) {
421: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
422: y[row] = shared[threadIdx.x];
423: }
424: }
425: }
427: static __global__ void matmult_seqsell_tiled_kernel2(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, PetscScalar *y)
428: {
429: __shared__ MatScalar shared[512];
430: PetscInt i, row, slice_id, row_in_slice;
431: /* multiple threads per row. */
432: row = blockIdx.x * blockDim.x + threadIdx.x;
433: if (row < nrows) {
434: slice_id = row / SLICE_HEIGHT;
435: row_in_slice = row % SLICE_HEIGHT;
437: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
438: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
439: __syncthreads();
440: if (threadIdx.y < 1) {
441: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
442: y[row] = shared[threadIdx.x];
443: }
444: }
445: }
447: static __global__ void matmultadd_seqsell_tiled_kernel6(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
448: {
449: __shared__ MatScalar shared[512];
450: PetscInt i, row, slice_id, row_in_slice;
451: /* multiple threads per row. */
452: row = blockIdx.x * blockDim.x + threadIdx.x;
453: if (row < nrows) {
454: slice_id = row / SLICE_HEIGHT;
455: row_in_slice = row % SLICE_HEIGHT;
457: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
458: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
459: __syncthreads();
460: if (threadIdx.y < 16) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 16) * blockDim.x + threadIdx.x];
461: __syncthreads();
462: if (threadIdx.y < 8) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 8) * blockDim.x + threadIdx.x];
463: __syncthreads();
464: if (threadIdx.y < 4) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 4) * blockDim.x + threadIdx.x];
465: __syncthreads();
466: if (threadIdx.y < 2) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 2) * blockDim.x + threadIdx.x];
467: __syncthreads();
468: if (threadIdx.y < 1) {
469: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
470: z[row] = y[row] + shared[threadIdx.x];
471: }
472: }
473: }
475: static __global__ void matmultadd_seqsell_tiled_kernel5(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
476: {
477: __shared__ MatScalar shared[512];
478: PetscInt i, row, slice_id, row_in_slice;
479: /* multiple threads per row. */
480: row = blockIdx.x * blockDim.x + threadIdx.x;
481: if (row < nrows) {
482: slice_id = row / SLICE_HEIGHT;
483: row_in_slice = row % SLICE_HEIGHT;
485: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
486: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
487: __syncthreads();
488: if (threadIdx.y < 8) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 8) * blockDim.x + threadIdx.x];
489: __syncthreads();
490: if (threadIdx.y < 4) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 4) * blockDim.x + threadIdx.x];
491: __syncthreads();
492: if (threadIdx.y < 2) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 2) * blockDim.x + threadIdx.x];
493: __syncthreads();
494: if (threadIdx.y < 1) {
495: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
496: z[row] = y[row] + shared[threadIdx.x];
497: }
498: }
499: }
501: static __global__ void matmultadd_seqsell_tiled_kernel4(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
502: {
503: __shared__ MatScalar shared[512];
504: PetscInt i, row, slice_id, row_in_slice;
505: /* multiple threads per row. */
506: row = blockIdx.x * blockDim.x + threadIdx.x;
507: if (row < nrows) {
508: slice_id = row / SLICE_HEIGHT;
509: row_in_slice = row % SLICE_HEIGHT;
511: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
512: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
513: __syncthreads();
514: if (threadIdx.y < 4) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 4) * blockDim.x + threadIdx.x];
515: __syncthreads();
516: if (threadIdx.y < 2) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 2) * blockDim.x + threadIdx.x];
517: __syncthreads();
518: if (threadIdx.y < 1) {
519: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
520: z[row] = y[row] + shared[threadIdx.x];
521: }
522: }
523: }
525: static __global__ void matmultadd_seqsell_tiled_kernel3(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
526: {
527: __shared__ MatScalar shared[512];
528: PetscInt i, row, slice_id, row_in_slice;
529: /* multiple threads per row. */
530: row = blockIdx.x * blockDim.x + threadIdx.x;
531: if (row < nrows) {
532: slice_id = row / SLICE_HEIGHT;
533: row_in_slice = row % SLICE_HEIGHT;
535: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
536: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
537: __syncthreads();
538: if (threadIdx.y < 2) shared[threadIdx.y * blockDim.x + threadIdx.x] += shared[(threadIdx.y + 2) * blockDim.x + threadIdx.x];
539: __syncthreads();
540: if (threadIdx.y < 1) {
541: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
542: z[row] = y[row] + shared[threadIdx.x];
543: }
544: }
545: }
547: static __global__ void matmultadd_seqsell_tiled_kernel2(PetscInt nrows, const PetscInt *acolidx, const MatScalar *aval, const PetscInt *sliidx, const PetscScalar *x, const PetscScalar *y, PetscScalar *z)
548: {
549: __shared__ MatScalar shared[512];
550: PetscInt i, row, slice_id, row_in_slice;
551: /* multiple threads per row. */
552: row = blockIdx.x * blockDim.x + threadIdx.x;
553: if (row < nrows) {
554: slice_id = row / SLICE_HEIGHT;
555: row_in_slice = row % SLICE_HEIGHT;
557: shared[threadIdx.y * blockDim.x + threadIdx.x] = 0.0;
558: for (i = sliidx[slice_id] + row_in_slice + SLICE_HEIGHT * threadIdx.y; i < sliidx[slice_id + 1]; i += SLICE_HEIGHT * blockDim.y) shared[threadIdx.y * blockDim.x + threadIdx.x] += aval[i] * x[acolidx[i]];
559: __syncthreads();
560: if (threadIdx.y < 1) {
561: shared[threadIdx.x] += shared[blockDim.x + threadIdx.x];
562: z[row] = y[row] + shared[threadIdx.x];
563: }
564: }
565: }
567: static PetscErrorCode MatMult_SeqSELLCUDA(Mat A, Vec xx, Vec yy)
568: {
569: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
570: Mat_SeqSELLCUDA *cudastruct = (Mat_SeqSELLCUDA *)A->spptr;
571: PetscScalar *y;
572: const PetscScalar *x;
573: PetscInt nrows = A->rmap->n, sliceheight = a->sliceheight;
574: MatScalar *aval;
575: PetscInt *acolidx;
576: PetscInt *sliidx;
577: PetscInt nblocks, blocksize = 512; /* blocksize must be multiple of SLICE_HEIGHT*32 */
578: dim3 block2(256, 2), block4(128, 4), block8(64, 8), block16(32, 16), block32(16, 32);
579: #if !PetscDefined(USE_COMPLEX)
580: PetscInt chunksperblock, nchunks, *chunk_slice_map;
581: PetscReal maxoveravg;
582: #endif
584: PetscFunctionBegin;
585: PetscCheck(32 % sliceheight == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "The kernel requires a slice height be a divisor of 32, but the input matrix has a slice height of %" PetscInt_FMT, sliceheight);
586: PetscCheck(!(cudastruct->kernelchoice >= 2 && cudastruct->kernelchoice <= 6 && sliceheight != SLICE_HEIGHT), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Kernel choices {2-6} requires the slice height of the matrix be 16, but the current slice height is %" PetscInt_FMT, sliceheight);
587: PetscCall(MatSeqSELLCUDACopyToGPU(A));
588: /* cudastruct may not be available until MatSeqSELLCUDACopyToGPU() is called */
589: aval = cudastruct->val;
590: acolidx = cudastruct->colidx;
591: sliidx = cudastruct->sliidx;
593: PetscCall(VecCUDAGetArrayRead(xx, &x));
594: PetscCall(VecCUDAGetArrayWrite(yy, &y));
595: PetscCall(PetscLogGpuTimeBegin());
597: switch (cudastruct->kernelchoice) {
598: #if !PetscDefined(USE_COMPLEX)
599: case 9:
600: nblocks = 1 + (nrows - 1) / sliceheight;
601: if (cudastruct->blocky == 2) {
602: matmult_seqsell_tiled_kernel9<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
603: } else if (cudastruct->blocky == 4) {
604: matmult_seqsell_tiled_kernel9<4><<<nblocks, dim3(32, 4)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
605: } else if (cudastruct->blocky == 8) {
606: matmult_seqsell_tiled_kernel9<8><<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
607: } else if (cudastruct->blocky == 16) {
608: matmult_seqsell_tiled_kernel9<16><<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
609: } else if (cudastruct->blocky == 32) {
610: matmult_seqsell_tiled_kernel9<32><<<nblocks, dim3(32, 32)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
611: } else {
612: matmult_seqsell_tiled_kernel9<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
613: }
614: break;
615: case 7:
616: nblocks = 1 + (nrows - 1) / (2 * sliceheight);
617: if (cudastruct->blocky == 2) {
618: matmult_seqsell_tiled_kernel7<<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
619: } else if (cudastruct->blocky == 4) {
620: matmult_seqsell_tiled_kernel7<<<nblocks, dim3(32, 4)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
621: } else if (cudastruct->blocky == 8) {
622: matmult_seqsell_tiled_kernel7<<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
623: } else if (cudastruct->blocky == 16) {
624: matmult_seqsell_tiled_kernel7<<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
625: } else if (cudastruct->blocky == 32) {
626: matmult_seqsell_tiled_kernel7<<<nblocks, dim3(32, 32)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
627: } else {
628: matmult_seqsell_tiled_kernel7<<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
629: }
630: break;
631: #endif
632: case 6:
633: nblocks = 1 + (nrows - 1) / (blocksize / 32); /* 1 slice per block if blocksize=512 */
634: matmult_seqsell_tiled_kernel6<<<nblocks, block32>>>(nrows, acolidx, aval, sliidx, x, y);
635: break;
636: case 5:
637: nblocks = 1 + (nrows - 1) / (blocksize / 16); /* 2 slices per block if blocksize=512*/
638: matmult_seqsell_tiled_kernel5<<<nblocks, block16>>>(nrows, acolidx, aval, sliidx, x, y);
639: break;
640: case 4:
641: nblocks = 1 + (nrows - 1) / (blocksize / 8); /* 4 slices per block if blocksize=512 */
642: matmult_seqsell_tiled_kernel4<<<nblocks, block8>>>(nrows, acolidx, aval, sliidx, x, y);
643: break;
644: case 3:
645: nblocks = 1 + (nrows - 1) / (blocksize / 4); /* 8 slices per block if blocksize=512 */
646: matmult_seqsell_tiled_kernel3<<<nblocks, block4>>>(nrows, acolidx, aval, sliidx, x, y);
647: break;
648: case 2: /* 16 slices per block if blocksize=512 */
649: nblocks = 1 + (nrows - 1) / (blocksize / 2);
650: matmult_seqsell_tiled_kernel2<<<nblocks, block2>>>(nrows, acolidx, aval, sliidx, x, y);
651: break;
652: case 1: /* 32 slices per block if blocksize=512 */
653: nblocks = 1 + (nrows - 1) / blocksize;
654: matmult_seqsell_basic_kernel<<<nblocks, blocksize>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
655: break;
656: #if !PetscDefined(USE_COMPLEX)
657: case 0:
658: maxoveravg = a->maxslicewidth / a->avgslicewidth;
659: if (maxoveravg > 12.0 && maxoveravg / nrows > 0.001) { /* important threshold */
660: /* each block handles approximately one slice */
661: PetscInt blocky = a->chunksize / 32;
662: nchunks = cudastruct->totalchunks;
663: chunksperblock = cudastruct->chunksperblock ? cudastruct->chunksperblock : 1 + (cudastruct->totalentries / cudastruct->totalslices - 1) / a->chunksize;
664: nblocks = 1 + (nchunks - 1) / chunksperblock;
665: chunk_slice_map = cudastruct->chunk_slice_map;
666: if (blocky == 2) {
667: matmult_seqsell_tiled_kernel8<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y);
668: } else if (blocky == 4) {
669: matmult_seqsell_tiled_kernel8<4><<<nblocks, dim3(32, 4)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y);
670: } else if (blocky == 8) {
671: matmult_seqsell_tiled_kernel8<8><<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y);
672: } else if (blocky == 16) {
673: matmult_seqsell_tiled_kernel8<16><<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y);
674: } else if (blocky == 32) {
675: matmult_seqsell_tiled_kernel8<32><<<nblocks, dim3(32, 32)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y);
676: } else {
677: matmult_seqsell_tiled_kernel8<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y);
678: }
679: } else {
680: PetscInt avgslicesize = sliceheight * a->avgslicewidth;
681: if (avgslicesize <= 432) {
682: if (sliceheight * a->maxslicewidth < 2048 && nrows > 100000) {
683: nblocks = 1 + (nrows - 1) / (2 * sliceheight); /* two slices per block */
684: matmult_seqsell_tiled_kernel7<<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
685: } else {
686: nblocks = 1 + (nrows - 1) / sliceheight;
687: matmult_seqsell_tiled_kernel9<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
688: }
689: } else if (avgslicesize <= 2400) {
690: nblocks = 1 + (nrows - 1) / sliceheight;
691: matmult_seqsell_tiled_kernel9<8><<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
692: } else {
693: nblocks = 1 + (nrows - 1) / sliceheight;
694: matmult_seqsell_tiled_kernel9<16><<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y);
695: }
696: }
697: break;
698: #endif
699: default:
700: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unsupported kernel choice %" PetscInt_FMT " for MatMult_SeqSELLCUDA.", cudastruct->kernelchoice);
701: }
702: PetscCall(PetscLogGpuTimeEnd());
703: PetscCall(VecCUDARestoreArrayRead(xx, &x));
704: PetscCall(VecCUDARestoreArrayWrite(yy, &y));
705: PetscCall(PetscLogGpuFlops(2.0 * a->nz - a->nonzerorowcnt));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: static PetscErrorCode MatMultAdd_SeqSELLCUDA(Mat A, Vec xx, Vec yy, Vec zz)
710: {
711: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
712: Mat_SeqSELLCUDA *cudastruct = (Mat_SeqSELLCUDA *)A->spptr;
713: PetscScalar *z;
714: const PetscScalar *y, *x;
715: PetscInt nrows = A->rmap->n, sliceheight = a->sliceheight;
716: MatScalar *aval;
717: PetscInt *acolidx;
718: PetscInt *sliidx;
719: #if !PetscDefined(USE_COMPLEX)
720: PetscReal maxoveravg;
721: PetscInt chunksperblock, nchunks, *chunk_slice_map;
722: PetscInt blocky = cudastruct->blocky;
723: #endif
725: PetscFunctionBegin;
726: PetscCheck(32 % sliceheight == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "The kernel requires a slice height be a divisor of 32, but the input matrix has a slice height of %" PetscInt_FMT, sliceheight);
727: PetscCheck(!(cudastruct->kernelchoice >= 2 && cudastruct->kernelchoice <= 6 && sliceheight != SLICE_HEIGHT), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Kernel choices {2-6} requires the slice height of the matrix be 16, but the current slice height is %" PetscInt_FMT, sliceheight);
728: PetscCall(MatSeqSELLCUDACopyToGPU(A));
729: aval = cudastruct->val;
730: acolidx = cudastruct->colidx;
731: sliidx = cudastruct->sliidx;
732: if (a->nz) {
733: PetscInt nblocks, blocksize = 512;
734: dim3 block2(256, 2), block4(128, 4), block8(64, 8), block16(32, 16), block32(16, 32);
735: PetscCall(VecCUDAGetArrayRead(xx, &x));
736: PetscCall(VecCUDAGetArrayRead(yy, &y));
737: PetscCall(VecCUDAGetArrayWrite(zz, &z));
738: PetscCall(PetscLogGpuTimeBegin());
740: switch (cudastruct->kernelchoice) {
741: #if !PetscDefined(USE_COMPLEX)
742: case 9:
743: nblocks = 1 + (nrows - 1) / sliceheight;
744: if (blocky == 2) {
745: matmultadd_seqsell_tiled_kernel9<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
746: } else if (blocky == 4) {
747: matmultadd_seqsell_tiled_kernel9<4><<<nblocks, dim3(32, 4)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
748: } else if (blocky == 8) {
749: matmultadd_seqsell_tiled_kernel9<8><<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
750: } else if (blocky == 16) {
751: matmultadd_seqsell_tiled_kernel9<16><<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
752: } else if (blocky == 32) {
753: matmultadd_seqsell_tiled_kernel9<32><<<nblocks, dim3(32, 32)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
754: } else {
755: matmultadd_seqsell_tiled_kernel9<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
756: }
757: break;
758: case 8:
759: /* each block handles approximately one slice */
760: nchunks = cudastruct->totalchunks;
761: blocky = a->chunksize / 32;
762: chunksperblock = cudastruct->chunksperblock ? cudastruct->chunksperblock : 1 + (cudastruct->totalentries / cudastruct->totalslices - 1) / a->chunksize;
763: nblocks = 1 + (nchunks - 1) / chunksperblock;
764: chunk_slice_map = cudastruct->chunk_slice_map;
765: if (blocky == 2) {
766: matmultadd_seqsell_tiled_kernel8<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
767: } else if (blocky == 4) {
768: matmultadd_seqsell_tiled_kernel8<4><<<nblocks, dim3(32, 4)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
769: } else if (blocky == 8) {
770: matmultadd_seqsell_tiled_kernel8<8><<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
771: } else if (blocky == 16) {
772: matmultadd_seqsell_tiled_kernel8<16><<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
773: } else if (blocky == 32) {
774: matmultadd_seqsell_tiled_kernel8<32><<<nblocks, dim3(32, 32)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
775: } else {
776: matmultadd_seqsell_tiled_kernel8<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
777: }
778: break;
779: case 7:
780: nblocks = 1 + (nrows - 1) / (2 * sliceheight);
781: if (blocky == 2) {
782: matmultadd_seqsell_tiled_kernel7<<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
783: } else if (blocky == 4) {
784: matmultadd_seqsell_tiled_kernel7<<<nblocks, dim3(32, 4)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
785: } else if (blocky == 8) {
786: matmultadd_seqsell_tiled_kernel7<<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
787: } else if (blocky == 16) {
788: matmultadd_seqsell_tiled_kernel7<<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
789: } else if (blocky == 32) {
790: matmultadd_seqsell_tiled_kernel7<<<nblocks, dim3(32, 32)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
791: } else {
792: matmultadd_seqsell_tiled_kernel7<<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
793: }
794: break;
795: #endif
796: case 6:
797: nblocks = 1 + (nrows - 1) / (blocksize / 32);
798: matmultadd_seqsell_tiled_kernel6<<<nblocks, block32>>>(nrows, acolidx, aval, sliidx, x, y, z);
799: break;
800: case 5:
801: nblocks = 1 + (nrows - 1) / (blocksize / 16);
802: matmultadd_seqsell_tiled_kernel5<<<nblocks, block16>>>(nrows, acolidx, aval, sliidx, x, y, z);
803: break;
804: case 4:
805: nblocks = 1 + (nrows - 1) / (blocksize / 8);
806: matmultadd_seqsell_tiled_kernel4<<<nblocks, block8>>>(nrows, acolidx, aval, sliidx, x, y, z);
807: break;
808: case 3:
809: nblocks = 1 + (nrows - 1) / (blocksize / 4);
810: matmultadd_seqsell_tiled_kernel3<<<nblocks, block4>>>(nrows, acolidx, aval, sliidx, x, y, z);
811: break;
812: case 2:
813: nblocks = 1 + (nrows - 1) / (blocksize / 2);
814: matmultadd_seqsell_tiled_kernel2<<<nblocks, block2>>>(nrows, acolidx, aval, sliidx, x, y, z);
815: break;
816: case 1:
817: nblocks = 1 + (nrows - 1) / blocksize;
818: matmultadd_seqsell_basic_kernel<<<nblocks, blocksize>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
819: break;
820: #if !PetscDefined(USE_COMPLEX)
821: case 0:
822: maxoveravg = a->maxslicewidth / a->avgslicewidth;
823: if (maxoveravg > 12.0 && maxoveravg / nrows > 0.001) { /* important threshold */
824: /* each block handles approximately one slice */
825: nchunks = cudastruct->totalchunks;
826: blocky = a->chunksize / 32;
827: chunksperblock = cudastruct->chunksperblock ? cudastruct->chunksperblock : 1 + (cudastruct->totalentries / cudastruct->totalslices - 1) / a->chunksize;
828: nblocks = 1 + (nchunks - 1) / chunksperblock;
829: chunk_slice_map = cudastruct->chunk_slice_map;
830: if (blocky == 2) {
831: matmultadd_seqsell_tiled_kernel8<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
832: } else if (blocky == 4) {
833: matmultadd_seqsell_tiled_kernel8<4><<<nblocks, dim3(32, 4)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
834: } else if (blocky == 8) {
835: matmultadd_seqsell_tiled_kernel8<8><<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
836: } else if (blocky == 16) {
837: matmultadd_seqsell_tiled_kernel8<16><<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
838: } else if (blocky == 32) {
839: matmultadd_seqsell_tiled_kernel8<32><<<nblocks, dim3(32, 32)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
840: } else {
841: matmultadd_seqsell_tiled_kernel8<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, chunksperblock, nchunks, chunk_slice_map, acolidx, aval, sliidx, x, y, z);
842: }
843: } else {
844: PetscInt avgslicesize = sliceheight * a->avgslicewidth;
845: if (avgslicesize <= 432) {
846: if (sliceheight * a->maxslicewidth < 2048 && nrows > 100000) {
847: nblocks = 1 + (nrows - 1) / (2 * sliceheight); /* two slices per block */
848: matmultadd_seqsell_tiled_kernel7<<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
849: } else {
850: nblocks = 1 + (nrows - 1) / sliceheight;
851: matmultadd_seqsell_tiled_kernel9<2><<<nblocks, dim3(32, 2)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
852: }
853: } else if (avgslicesize <= 2400) {
854: nblocks = 1 + (nrows - 1) / sliceheight;
855: matmultadd_seqsell_tiled_kernel9<8><<<nblocks, dim3(32, 8)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
856: } else {
857: nblocks = 1 + (nrows - 1) / sliceheight;
858: matmultadd_seqsell_tiled_kernel9<16><<<nblocks, dim3(32, 16)>>>(nrows, sliceheight, acolidx, aval, sliidx, x, y, z);
859: }
860: }
861: break;
862: #endif
863: default:
864: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "unsupported kernel choice %" PetscInt_FMT " for MatMult_SeqSELLCUDA.", cudastruct->kernelchoice);
865: }
866: PetscCall(PetscLogGpuTimeEnd());
867: PetscCall(VecCUDARestoreArrayRead(xx, &x));
868: PetscCall(VecCUDARestoreArrayRead(yy, &y));
869: PetscCall(VecCUDARestoreArrayWrite(zz, &z));
870: PetscCall(PetscLogGpuFlops(2.0 * a->nz));
871: } else {
872: PetscCall(VecCopy(yy, zz));
873: }
874: PetscFunctionReturn(PETSC_SUCCESS);
875: }
877: static PetscErrorCode MatSetFromOptions_SeqSELLCUDA(Mat A, PetscOptionItems PetscOptionsObject)
878: {
879: Mat_SeqSELLCUDA *cudastruct = (Mat_SeqSELLCUDA *)A->spptr;
880: PetscInt kernel, blocky;
881: PetscBool flg;
883: PetscFunctionBegin;
884: PetscOptionsHeadBegin(PetscOptionsObject, "SeqSELLCUDA options");
885: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_sell_spmv_cuda_blocky", &blocky, &flg));
886: if (flg) {
887: PetscCheck(blocky == 2 || blocky == 4 || blocky == 8 || blocky == 16 || blocky == 32, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unsupported blocky: %" PetscInt_FMT " it should be in {2,4,8,16,32}", blocky);
888: cudastruct->blocky = blocky;
889: }
890: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_sell_spmv_cuda_kernel", &kernel, &flg));
891: if (flg) {
892: PetscCheck(kernel >= 0 && kernel <= 9, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Wrong kernel choice: %" PetscInt_FMT " it should be in [0,9]", kernel);
893: cudastruct->kernelchoice = kernel;
894: if (kernel == 8) PetscCall(PetscOptionsGetInt(NULL, NULL, "-mat_sell_spmv_cuda_chunksperblock", &cudastruct->chunksperblock, &flg));
895: }
896: PetscOptionsHeadEnd();
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: PETSC_INTERN PetscErrorCode MatAssemblyEnd_SpMV_Preprocessing_Private(Mat A)
901: {
902: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
904: PetscFunctionBegin;
905: PetscCall(MatSeqSELLGetAvgSliceWidth(A, &a->avgslicewidth));
906: PetscCall(MatSeqSELLGetMaxSliceWidth(A, &a->maxslicewidth));
907: PetscCall(MatSeqSELLGetFillRatio(A, &a->fillratio));
908: PetscFunctionReturn(PETSC_SUCCESS);
909: }
911: static PetscErrorCode MatAssemblyEnd_SeqSELLCUDA(Mat A, MatAssemblyType mode)
912: {
913: PetscFunctionBegin;
914: PetscCall(MatAssemblyEnd_SeqSELL(A, mode));
915: PetscCall(MatAssemblyEnd_SpMV_Preprocessing_Private(A));
916: if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(PETSC_SUCCESS);
917: if (A->factortype == MAT_FACTOR_NONE) PetscCall(MatSeqSELLCUDACopyToGPU(A));
918: A->ops->mult = MatMult_SeqSELLCUDA;
919: A->ops->multadd = MatMultAdd_SeqSELLCUDA;
920: PetscFunctionReturn(PETSC_SUCCESS);
921: }
923: static PetscErrorCode MatZeroEntries_SeqSELLCUDA(Mat A)
924: {
925: PetscBool both = PETSC_FALSE;
926: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
928: PetscFunctionBegin;
929: if (A->factortype == MAT_FACTOR_NONE) {
930: Mat_SeqSELLCUDA *cudastruct = (Mat_SeqSELLCUDA *)A->spptr;
931: if (cudastruct->val) {
932: both = PETSC_TRUE;
933: PetscCallCUDA(cudaMemset(cudastruct->val, 0, a->sliidx[a->totalslices] * sizeof(*cudastruct->val)));
934: }
935: }
936: PetscCall(PetscArrayzero(a->val, a->sliidx[a->totalslices]));
937: if (both) A->offloadmask = PETSC_OFFLOAD_BOTH;
938: else A->offloadmask = PETSC_OFFLOAD_CPU;
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }
942: static PetscErrorCode MatDestroy_SeqSELLCUDA(Mat A)
943: {
944: PetscFunctionBegin;
945: if (A->factortype == MAT_FACTOR_NONE) PetscCall(MatSeqSELLCUDA_Destroy((Mat_SeqSELLCUDA **)&A->spptr));
946: PetscCall(MatDestroy_SeqSELL(A));
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: PETSC_INTERN PetscErrorCode MatConvert_SeqSELL_SeqSELLCUDA(Mat);
951: static PetscErrorCode MatDuplicate_SeqSELLCUDA(Mat A, MatDuplicateOption cpvalues, Mat *B)
952: {
953: PetscFunctionBegin;
954: PetscCall(MatDuplicate_SeqSELL(A, cpvalues, B));
955: PetscCall(MatConvert_SeqSELL_SeqSELLCUDA(*B));
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
959: PETSC_INTERN PetscErrorCode MatConvert_SeqSELL_SeqSELLCUDA(Mat B)
960: {
961: Mat_SeqSELLCUDA *cudastruct;
963: PetscFunctionBegin;
964: PetscCall(PetscFree(B->defaultvectype));
965: PetscCall(PetscStrallocpy(VECCUDA, &B->defaultvectype));
967: if (!B->spptr) {
968: if (B->factortype == MAT_FACTOR_NONE) {
969: PetscCall(PetscNew(&cudastruct));
970: B->spptr = cudastruct;
971: }
972: }
974: B->ops->assemblyend = MatAssemblyEnd_SeqSELLCUDA;
975: B->ops->destroy = MatDestroy_SeqSELLCUDA;
976: B->ops->setfromoptions = MatSetFromOptions_SeqSELLCUDA;
977: B->ops->mult = MatMult_SeqSELLCUDA;
978: B->ops->multadd = MatMultAdd_SeqSELLCUDA;
979: B->ops->duplicate = MatDuplicate_SeqSELLCUDA;
980: B->ops->zeroentries = MatZeroEntries_SeqSELLCUDA;
982: /* No need to assemble SeqSELL, but need to do the preprocessing for SpMV */
983: PetscCall(MatAssemblyEnd_SpMV_Preprocessing_Private(B));
985: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQSELLCUDA));
986: B->offloadmask = PETSC_OFFLOAD_UNALLOCATED;
987: PetscFunctionReturn(PETSC_SUCCESS);
988: }
990: /*MC
991: MATSEQSELLCUDA - MATSELLCUDA = "(seq)sellcuda" - A matrix type to be used for sparse matrices on NVIDIA GPUs.
993: Options Database Keys:
994: + -mat_type seqsellcuda - sets the matrix type to "seqsellcuda" during a call to `MatSetFromOptions()`
995: . -mat_sell_spmv_cuda_kernel - selects a spmv kernel for MatSELLCUDA
996: - -mat_sell_spmv_cuda_blocky - sets the y dimension of the block size of the spmv kernels. These kernels use a 2D block with the x dimension being 32
998: Level: beginner
1000: .seealso: [](ch_matrices), `Mat`, `MATSELLCUDA`
1001: M*/
1003: PETSC_EXTERN PetscErrorCode MatCreate_SeqSELLCUDA(Mat B)
1004: {
1005: PetscFunctionBegin;
1006: PetscCall(MatCreate_SeqSELL(B));
1007: PetscCall(MatConvert_SeqSELL_SeqSELLCUDA(B));
1008: PetscCall(MatSetFromOptions(B));
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }