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