Actual source code: sell.c
1: /*
2: Defines the basic matrix operations for the SELL matrix storage format.
3: */
4: #include <../src/mat/impls/sell/seq/sell.h>
5: #include <petscblaslapack.h>
6: #include <petsc/private/kernels/blocktranspose.h>
8: static PetscBool cited = PETSC_FALSE;
9: static const char citation[] = "@inproceedings{ZhangELLPACK2018,\n"
10: " author = {Hong Zhang and Richard T. Mills and Karl Rupp and Barry F. Smith},\n"
11: " title = {Vectorized Parallel Sparse Matrix-Vector Multiplication in {PETSc} Using {AVX-512}},\n"
12: " booktitle = {Proceedings of the 47th International Conference on Parallel Processing},\n"
13: " year = 2018\n"
14: "}\n";
16: #if PetscDefined(HAVE_IMMINTRIN_H) && (defined(__AVX512F__) || (defined(__AVX2__) && defined(__FMA__)) || defined(__AVX__)) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
18: #include <immintrin.h>
20: #if !defined(_MM_SCALE_8)
21: #define _MM_SCALE_8 8
22: #endif
24: #if defined(__AVX512F__)
25: /* these do not work
26: vec_idx = _mm512_loadunpackhi_epi32(vec_idx,acolidx);
27: vec_vals = _mm512_loadunpackhi_pd(vec_vals,aval);
28: */
29: #define AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y) \
30: /* if the mask bit is set, copy from acolidx, otherwise from vec_idx */ \
31: vec_idx = _mm256_loadu_si256((__m256i const *)acolidx); \
32: vec_vals = _mm512_loadu_pd(aval); \
33: vec_x = _mm512_i32gather_pd(vec_idx, x, _MM_SCALE_8); \
34: vec_y = _mm512_fmadd_pd(vec_x, vec_vals, vec_y)
35: #elif defined(__AVX2__) && defined(__FMA__)
36: #define AVX2_Mult_Private(vec_idx, vec_x, vec_vals, vec_y) \
37: vec_vals = _mm256_loadu_pd(aval); \
38: vec_idx = _mm_loadu_si128((__m128i const *)acolidx); /* SSE2 */ \
39: vec_x = _mm256_i32gather_pd(x, vec_idx, _MM_SCALE_8); \
40: vec_y = _mm256_fmadd_pd(vec_x, vec_vals, vec_y)
41: #endif
42: #endif /* PETSC_HAVE_IMMINTRIN_H */
44: /*@
45: MatSeqSELLSetPreallocation - For good matrix assembly performance
46: the user should preallocate the matrix storage by setting the parameter `nz`
47: (or the array `nnz`).
49: Collective
51: Input Parameters:
52: + B - The `MATSEQSELL` matrix
53: . rlenmax - number of nonzeros per row (same for all rows), ignored if `rlen` is provided
54: - rlen - array containing the number of nonzeros in the various rows (possibly different for each row) or `NULL`
56: Level: intermediate
58: Notes:
59: Specify the preallocated storage with either `rlenmax` or `rlen` (not both).
60: Set `rlenmax` = `PETSC_DEFAULT` and `rlen` = `NULL` for PETSc to control dynamic memory
61: allocation.
63: You can call `MatGetInfo()` to get information on how effective the preallocation was;
64: for example the fields mallocs,nz_allocated,nz_used,nz_unneeded;
65: You can also run with the option `-info` and look for messages with the string
66: malloc in them to see if additional memory allocation was needed.
68: Developer Notes:
69: Use `rlenmax` of `MAT_SKIP_ALLOCATION` to not allocate any space for the matrix
70: entries or columns indices.
72: The maximum number of nonzeos in any row should be as accurate as possible.
73: If it is underestimated, you will get bad performance due to reallocation
74: (`MatSeqXSELLReallocateSELL()`).
76: .seealso: `Mat`, `MATSEQSELL`, `MATSELL`, `MatCreate()`, `MatCreateSELL()`, `MatSetValues()`, `MatGetInfo()`
77: @*/
78: PetscErrorCode MatSeqSELLSetPreallocation(Mat B, PetscInt rlenmax, const PetscInt rlen[])
79: {
80: PetscFunctionBegin;
83: PetscTryMethod(B, "MatSeqSELLSetPreallocation_C", (Mat, PetscInt, const PetscInt[]), (B, rlenmax, rlen));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: PetscErrorCode MatSeqSELLSetPreallocation_SeqSELL(Mat B, PetscInt maxallocrow, const PetscInt rlen[])
88: {
89: Mat_SeqSELL *b;
90: PetscInt i, j, totalslices;
91: #if PetscDefined(HAVE_CUPM)
92: PetscInt rlenmax = 0;
93: #endif
94: PetscBool skipallocation = PETSC_FALSE, realalloc = PETSC_FALSE;
96: PetscFunctionBegin;
97: if (maxallocrow >= 0 || rlen) realalloc = PETSC_TRUE;
98: if (maxallocrow == MAT_SKIP_ALLOCATION) {
99: skipallocation = PETSC_TRUE;
100: maxallocrow = 0;
101: }
103: PetscCall(PetscLayoutSetUp(B->rmap));
104: PetscCall(PetscLayoutSetUp(B->cmap));
106: /* FIXME: if one preallocates more space than needed, the matrix does not shrink automatically, but for best performance it should */
107: if (maxallocrow == PETSC_DEFAULT || maxallocrow == PETSC_DECIDE) maxallocrow = 5;
108: PetscCheck(maxallocrow >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "maxallocrow cannot be less than 0: value %" PetscInt_FMT, maxallocrow);
109: if (rlen) {
110: for (i = 0; i < B->rmap->n; i++) {
111: PetscCheck(rlen[i] >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "rlen cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT, i, rlen[i]);
112: PetscCheck(rlen[i] <= B->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "rlen cannot be greater than row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT, i, rlen[i], B->cmap->n);
113: }
114: }
116: B->preallocated = PETSC_TRUE;
118: b = (Mat_SeqSELL *)B->data;
120: if (!b->sliceheight) { /* not set yet */
121: #if PetscDefined(HAVE_CUPM)
122: b->sliceheight = 16;
123: #else
124: b->sliceheight = 8;
125: #endif
126: }
127: totalslices = PetscCeilInt(B->rmap->n, b->sliceheight);
128: b->totalslices = totalslices;
129: if (!skipallocation) {
130: if (B->rmap->n % b->sliceheight) PetscCall(PetscInfo(B, "Padding rows to the SEQSELL matrix because the number of rows is not the multiple of the slice height (value %" PetscInt_FMT ")\n", B->rmap->n));
132: if (!b->sliidx) { /* sliidx gives the starting index of each slice, the last element is the total space allocated */
133: PetscCall(PetscMalloc1(totalslices + 1, &b->sliidx));
134: }
135: if (!rlen) { /* if rlen is not provided, allocate same space for all the slices */
136: if (maxallocrow == PETSC_DEFAULT || maxallocrow == PETSC_DECIDE) maxallocrow = 10;
137: else if (maxallocrow < 0) maxallocrow = 1;
138: #if PetscDefined(HAVE_CUPM)
139: rlenmax = maxallocrow;
140: /* Pad the slice to DEVICE_MEM_ALIGN */
141: while (b->sliceheight * maxallocrow % DEVICE_MEM_ALIGN) maxallocrow++;
142: #endif
143: for (i = 0; i <= totalslices; i++) b->sliidx[i] = b->sliceheight * i * maxallocrow;
144: } else {
145: #if PetscDefined(HAVE_CUPM)
146: PetscInt mul = DEVICE_MEM_ALIGN / b->sliceheight;
147: #endif
148: maxallocrow = 0;
149: b->sliidx[0] = 0;
150: for (i = 1; i < totalslices; i++) {
151: b->sliidx[i] = 0;
152: for (j = 0; j < b->sliceheight; j++) b->sliidx[i] = PetscMax(b->sliidx[i], rlen[b->sliceheight * (i - 1) + j]);
153: #if PetscDefined(HAVE_CUPM)
154: if (mul != 0) { /* Pad the slice to DEVICE_MEM_ALIGN if sliceheight < DEVICE_MEM_ALIGN */
155: rlenmax = PetscMax(b->sliidx[i], rlenmax);
156: b->sliidx[i] = ((b->sliidx[i] - 1) / mul + 1) * mul;
157: }
158: #endif
159: maxallocrow = PetscMax(b->sliidx[i], maxallocrow);
160: PetscCall(PetscIntSumError(b->sliidx[i - 1], b->sliceheight * b->sliidx[i], &b->sliidx[i]));
161: }
162: /* last slice */
163: b->sliidx[totalslices] = 0;
164: for (j = b->sliceheight * (totalslices - 1); j < B->rmap->n; j++) b->sliidx[totalslices] = PetscMax(b->sliidx[totalslices], rlen[j]);
165: #if PetscDefined(HAVE_CUPM)
166: if (mul != 0) {
167: rlenmax = PetscMax(b->sliidx[i], rlenmax);
168: b->sliidx[totalslices] = ((b->sliidx[totalslices] - 1) / mul + 1) * mul;
169: }
170: #endif
171: maxallocrow = PetscMax(b->sliidx[totalslices], maxallocrow);
172: b->sliidx[totalslices] = b->sliidx[totalslices - 1] + b->sliceheight * b->sliidx[totalslices];
173: }
175: /* allocate space for val, colidx, rlen */
176: /* FIXME: should B's old memory be unlogged? */
177: PetscCall(MatSeqXSELLFreeSELL(B, &b->val, &b->colidx));
178: /* FIXME: assuming an element of the bit array takes 8 bits */
179: PetscCall(PetscMalloc2(b->sliidx[totalslices], &b->val, b->sliidx[totalslices], &b->colidx));
180: /* b->rlen will count nonzeros in each row so far. We dont copy rlen to b->rlen because the matrix has not been set. */
181: PetscCall(PetscCalloc1(b->sliceheight * totalslices, &b->rlen));
183: b->singlemalloc = PETSC_TRUE;
184: b->free_val = PETSC_TRUE;
185: b->free_colidx = PETSC_TRUE;
186: } else {
187: b->free_val = PETSC_FALSE;
188: b->free_colidx = PETSC_FALSE;
189: }
191: b->nz = 0;
192: b->maxallocrow = maxallocrow;
193: #if PetscDefined(HAVE_CUPM)
194: b->rlenmax = rlenmax;
195: #else
196: b->rlenmax = maxallocrow;
197: #endif
198: b->maxallocmat = b->sliidx[totalslices];
199: B->info.nz_unneeded = (double)b->maxallocmat;
200: if (realalloc) PetscCall(MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
201: PetscFunctionReturn(PETSC_SUCCESS);
202: }
204: static PetscErrorCode MatGetRow_SeqSELL(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
205: {
206: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
207: PetscInt shift;
209: PetscFunctionBegin;
210: PetscCheck(row >= 0 && row < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row %" PetscInt_FMT " out of range", row);
211: if (nz) *nz = a->rlen[row];
212: shift = a->sliidx[row / a->sliceheight] + (row % a->sliceheight);
213: if (!a->getrowcols) PetscCall(PetscMalloc2(a->rlenmax, &a->getrowcols, a->rlenmax, &a->getrowvals));
214: if (idx) {
215: PetscInt j;
216: for (j = 0; j < a->rlen[row]; j++) a->getrowcols[j] = a->colidx[shift + a->sliceheight * j];
217: *idx = a->getrowcols;
218: }
219: if (v) {
220: for (PetscInt j = 0; j < a->rlen[row]; j++) a->getrowvals[j] = a->val[shift + a->sliceheight * j];
221: *v = a->getrowvals;
222: }
223: PetscFunctionReturn(PETSC_SUCCESS);
224: }
226: static PetscErrorCode MatRestoreRow_SeqSELL(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
227: {
228: PetscFunctionBegin;
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: PetscErrorCode MatConvert_SeqSELL_SeqAIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
233: {
234: Mat B;
235: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
236: PetscInt i;
238: PetscFunctionBegin;
239: if (reuse == MAT_REUSE_MATRIX) {
240: B = *newmat;
241: PetscCall(MatZeroEntries(B));
242: } else {
243: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
244: PetscCall(MatSetSizes(B, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
245: PetscCall(MatSetType(B, MATSEQAIJ));
246: PetscCall(MatSeqAIJSetPreallocation(B, 0, a->rlen));
247: }
249: for (i = 0; i < A->rmap->n; i++) {
250: PetscInt nz = 0, *cols = NULL;
251: PetscScalar *vals = NULL;
253: PetscCall(MatGetRow_SeqSELL(A, i, &nz, &cols, &vals));
254: PetscCall(MatSetValues(B, 1, &i, nz, cols, vals, INSERT_VALUES));
255: PetscCall(MatRestoreRow_SeqSELL(A, i, &nz, &cols, &vals));
256: }
258: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
259: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
260: B->rmap->bs = A->rmap->bs;
262: if (reuse == MAT_INPLACE_MATRIX) {
263: PetscCall(MatHeaderReplace(A, &B));
264: } else {
265: *newmat = B;
266: }
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: #include <../src/mat/impls/aij/seq/aij.h>
272: PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
273: {
274: Mat B;
275: Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data;
276: PetscInt *ai = a->i, m = A->rmap->N, n = A->cmap->N, i, *rowlengths, row, ncols;
277: const PetscInt *cols;
278: const PetscScalar *vals;
280: PetscFunctionBegin;
281: if (reuse == MAT_REUSE_MATRIX) {
282: B = *newmat;
283: } else {
284: if (PetscDefined(USE_DEBUG) || !a->ilen) {
285: PetscCall(PetscMalloc1(m, &rowlengths));
286: for (i = 0; i < m; i++) rowlengths[i] = ai[i + 1] - ai[i];
287: }
288: if (PetscDefined(USE_DEBUG) && a->ilen) {
289: PetscBool eq;
290: PetscCall(PetscArraycmp(rowlengths, a->ilen, m, &eq));
291: PetscCheck(eq, PETSC_COMM_SELF, PETSC_ERR_PLIB, "SeqAIJ ilen array incorrect");
292: PetscCall(PetscFree(rowlengths));
293: rowlengths = a->ilen;
294: } else if (a->ilen) rowlengths = a->ilen;
295: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &B));
296: PetscCall(MatSetSizes(B, m, n, m, n));
297: PetscCall(MatSetType(B, MATSEQSELL));
298: PetscCall(MatSeqSELLSetPreallocation(B, 0, rowlengths));
299: if (rowlengths != a->ilen) PetscCall(PetscFree(rowlengths));
300: }
302: for (row = 0; row < m; row++) {
303: PetscCall(MatGetRow_SeqAIJ(A, row, &ncols, (PetscInt **)&cols, (PetscScalar **)&vals));
304: PetscCall(MatSetValues_SeqSELL(B, 1, &row, ncols, cols, vals, INSERT_VALUES));
305: PetscCall(MatRestoreRow_SeqAIJ(A, row, &ncols, (PetscInt **)&cols, (PetscScalar **)&vals));
306: }
307: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
308: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
309: B->rmap->bs = A->rmap->bs;
311: if (reuse == MAT_INPLACE_MATRIX) {
312: PetscCall(MatHeaderReplace(A, &B));
313: } else {
314: *newmat = B;
315: }
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: PetscErrorCode MatMult_SeqSELL(Mat A, Vec xx, Vec yy)
320: {
321: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
322: PetscScalar *y;
323: const PetscScalar *x;
324: const MatScalar *aval = a->val;
325: PetscInt totalslices = a->totalslices;
326: const PetscInt *acolidx = a->colidx;
327: PetscInt i, j;
328: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX512F__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
329: __m512d vec_x, vec_y, vec_vals;
330: __m256i vec_idx;
331: __mmask8 mask;
332: __m512d vec_x2, vec_y2, vec_vals2, vec_x3, vec_y3, vec_vals3, vec_x4, vec_y4, vec_vals4;
333: __m256i vec_idx2, vec_idx3, vec_idx4;
334: #elif PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
335: __m128i vec_idx;
336: __m256d vec_x, vec_y, vec_y2, vec_vals;
337: MatScalar yval;
338: PetscInt r, rows_left, row, nnz_in_row;
339: #elif PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
340: __m128d vec_x_tmp;
341: __m256d vec_x, vec_y, vec_y2, vec_vals;
342: MatScalar yval;
343: PetscInt r, rows_left, row, nnz_in_row;
344: #else
345: PetscInt k, sliceheight = a->sliceheight;
346: PetscScalar *sum;
347: #endif
349: #if PetscDefined(HAVE_PRAGMA_DISJOINT)
350: #pragma disjoint(*x, *y, *aval)
351: #endif
353: PetscFunctionBegin;
354: PetscCall(VecGetArrayRead(xx, &x));
355: PetscCall(VecGetArray(yy, &y));
356: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX512F__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
357: PetscCheck(a->sliceheight == 8, PETSC_COMM_SELF, PETSC_ERR_SUP, "The kernel requires a slice height of 8, but the input matrix has a slice height of %" PetscInt_FMT, a->sliceheight);
358: for (i = 0; i < totalslices; i++) { /* loop over slices */
359: PetscPrefetchBlock(acolidx, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
360: PetscPrefetchBlock(aval, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
362: vec_y = _mm512_setzero_pd();
363: vec_y2 = _mm512_setzero_pd();
364: vec_y3 = _mm512_setzero_pd();
365: vec_y4 = _mm512_setzero_pd();
367: j = a->sliidx[i] >> 3; /* 8 bytes are read at each time, corresponding to a slice column */
368: switch ((a->sliidx[i + 1] - a->sliidx[i]) / 8 & 3) {
369: case 3:
370: AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
371: acolidx += 8;
372: aval += 8;
373: AVX512_Mult_Private(vec_idx2, vec_x2, vec_vals2, vec_y2);
374: acolidx += 8;
375: aval += 8;
376: AVX512_Mult_Private(vec_idx3, vec_x3, vec_vals3, vec_y3);
377: acolidx += 8;
378: aval += 8;
379: j += 3;
380: break;
381: case 2:
382: AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
383: acolidx += 8;
384: aval += 8;
385: AVX512_Mult_Private(vec_idx2, vec_x2, vec_vals2, vec_y2);
386: acolidx += 8;
387: aval += 8;
388: j += 2;
389: break;
390: case 1:
391: AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
392: acolidx += 8;
393: aval += 8;
394: j += 1;
395: break;
396: }
397: #pragma novector
398: for (; j < (a->sliidx[i + 1] >> 3); j += 4) {
399: AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
400: acolidx += 8;
401: aval += 8;
402: AVX512_Mult_Private(vec_idx2, vec_x2, vec_vals2, vec_y2);
403: acolidx += 8;
404: aval += 8;
405: AVX512_Mult_Private(vec_idx3, vec_x3, vec_vals3, vec_y3);
406: acolidx += 8;
407: aval += 8;
408: AVX512_Mult_Private(vec_idx4, vec_x4, vec_vals4, vec_y4);
409: acolidx += 8;
410: aval += 8;
411: }
413: vec_y = _mm512_add_pd(vec_y, vec_y2);
414: vec_y = _mm512_add_pd(vec_y, vec_y3);
415: vec_y = _mm512_add_pd(vec_y, vec_y4);
416: if (i == totalslices - 1 && A->rmap->n & 0x07) { /* if last slice has padding rows */
417: mask = (__mmask8)(0xff >> (8 - (A->rmap->n & 0x07)));
418: _mm512_mask_storeu_pd(&y[8 * i], mask, vec_y);
419: } else {
420: _mm512_storeu_pd(&y[8 * i], vec_y);
421: }
422: }
423: #elif PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX2__) && defined(__FMA__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
424: PetscCheck(a->sliceheight == 8, PETSC_COMM_SELF, PETSC_ERR_SUP, "The kernel requires a slice height of 8, but the input matrix has a slice height of %" PetscInt_FMT, a->sliceheight);
425: for (i = 0; i < totalslices; i++) { /* loop over full slices */
426: PetscPrefetchBlock(acolidx, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
427: PetscPrefetchBlock(aval, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
429: /* last slice may have padding rows. Don't use vectorization. */
430: if (i == totalslices - 1 && (A->rmap->n & 0x07)) {
431: rows_left = A->rmap->n - 8 * i;
432: for (r = 0; r < rows_left; ++r) {
433: yval = (MatScalar)0;
434: row = 8 * i + r;
435: nnz_in_row = a->rlen[row];
436: for (j = 0; j < nnz_in_row; ++j) yval += aval[8 * j + r] * x[acolidx[8 * j + r]];
437: y[row] = yval;
438: }
439: break;
440: }
442: vec_y = _mm256_setzero_pd();
443: vec_y2 = _mm256_setzero_pd();
445: /* Process slice of height 8 (512 bits) via two subslices of height 4 (256 bits) via AVX */
446: #pragma novector
447: #pragma unroll(2)
448: for (j = a->sliidx[i]; j < a->sliidx[i + 1]; j += 8) {
449: AVX2_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
450: aval += 4;
451: acolidx += 4;
452: AVX2_Mult_Private(vec_idx, vec_x, vec_vals, vec_y2);
453: aval += 4;
454: acolidx += 4;
455: }
457: _mm256_storeu_pd(y + i * 8, vec_y);
458: _mm256_storeu_pd(y + i * 8 + 4, vec_y2);
459: }
460: #elif PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
461: PetscCheck(a->sliceheight == 8, PETSC_COMM_SELF, PETSC_ERR_SUP, "The kernel requires a slice height of 8, but the input matrix has a slice height of %" PetscInt_FMT, a->sliceheight);
462: for (i = 0; i < totalslices; i++) { /* loop over full slices */
463: PetscPrefetchBlock(acolidx, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
464: PetscPrefetchBlock(aval, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
466: vec_y = _mm256_setzero_pd();
467: vec_y2 = _mm256_setzero_pd();
469: /* last slice may have padding rows. Don't use vectorization. */
470: if (i == totalslices - 1 && (A->rmap->n & 0x07)) {
471: rows_left = A->rmap->n - 8 * i;
472: for (r = 0; r < rows_left; ++r) {
473: yval = (MatScalar)0;
474: row = 8 * i + r;
475: nnz_in_row = a->rlen[row];
476: for (j = 0; j < nnz_in_row; ++j) yval += aval[8 * j + r] * x[acolidx[8 * j + r]];
477: y[row] = yval;
478: }
479: break;
480: }
482: /* Process slice of height 8 (512 bits) via two subslices of height 4 (256 bits) via AVX */
483: #pragma novector
484: #pragma unroll(2)
485: for (j = a->sliidx[i]; j < a->sliidx[i + 1]; j += 8) {
486: vec_vals = _mm256_loadu_pd(aval);
487: vec_x_tmp = _mm_setzero_pd();
488: vec_x_tmp = _mm_loadl_pd(vec_x_tmp, x + *acolidx++);
489: vec_x_tmp = _mm_loadh_pd(vec_x_tmp, x + *acolidx++);
490: vec_x = _mm256_insertf128_pd(vec_x, vec_x_tmp, 0);
491: vec_x_tmp = _mm_loadl_pd(vec_x_tmp, x + *acolidx++);
492: vec_x_tmp = _mm_loadh_pd(vec_x_tmp, x + *acolidx++);
493: vec_x = _mm256_insertf128_pd(vec_x, vec_x_tmp, 1);
494: vec_y = _mm256_add_pd(_mm256_mul_pd(vec_x, vec_vals), vec_y);
495: aval += 4;
497: vec_vals = _mm256_loadu_pd(aval);
498: vec_x_tmp = _mm_loadl_pd(vec_x_tmp, x + *acolidx++);
499: vec_x_tmp = _mm_loadh_pd(vec_x_tmp, x + *acolidx++);
500: vec_x = _mm256_insertf128_pd(vec_x, vec_x_tmp, 0);
501: vec_x_tmp = _mm_loadl_pd(vec_x_tmp, x + *acolidx++);
502: vec_x_tmp = _mm_loadh_pd(vec_x_tmp, x + *acolidx++);
503: vec_x = _mm256_insertf128_pd(vec_x, vec_x_tmp, 1);
504: vec_y2 = _mm256_add_pd(_mm256_mul_pd(vec_x, vec_vals), vec_y2);
505: aval += 4;
506: }
508: _mm256_storeu_pd(y + i * 8, vec_y);
509: _mm256_storeu_pd(y + i * 8 + 4, vec_y2);
510: }
511: #else
512: PetscCall(PetscMalloc1(sliceheight, &sum));
513: for (i = 0; i < totalslices; i++) { /* loop over slices */
514: for (j = 0; j < sliceheight; j++) {
515: sum[j] = 0.0;
516: for (k = a->sliidx[i] + j; k < a->sliidx[i + 1]; k += sliceheight) sum[j] += aval[k] * x[acolidx[k]];
517: }
518: if (i == totalslices - 1 && (A->rmap->n % sliceheight)) { /* if last slice has padding rows */
519: for (j = 0; j < (A->rmap->n % sliceheight); j++) y[sliceheight * i + j] = sum[j];
520: } else {
521: for (j = 0; j < sliceheight; j++) y[sliceheight * i + j] = sum[j];
522: }
523: }
524: PetscCall(PetscFree(sum));
525: #endif
527: PetscCall(PetscLogFlops(2.0 * a->nz - a->nonzerorowcnt)); /* theoretical minimal FLOPs */
528: PetscCall(VecRestoreArrayRead(xx, &x));
529: PetscCall(VecRestoreArray(yy, &y));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h>
534: PetscErrorCode MatMultAdd_SeqSELL(Mat A, Vec xx, Vec yy, Vec zz)
535: {
536: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
537: PetscScalar *y, *z;
538: const PetscScalar *x;
539: const MatScalar *aval = a->val;
540: PetscInt totalslices = a->totalslices;
541: const PetscInt *acolidx = a->colidx;
542: PetscInt i, j;
543: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX512F__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
544: __m512d vec_x, vec_y, vec_vals;
545: __m256i vec_idx;
546: __mmask8 mask = 0;
547: __m512d vec_x2, vec_y2, vec_vals2, vec_x3, vec_y3, vec_vals3, vec_x4, vec_y4, vec_vals4;
548: __m256i vec_idx2, vec_idx3, vec_idx4;
549: #elif PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
550: __m128d vec_x_tmp;
551: __m256d vec_x, vec_y, vec_y2, vec_vals;
552: MatScalar yval;
553: PetscInt r, row, nnz_in_row;
554: #else
555: PetscInt k, sliceheight = a->sliceheight;
556: PetscScalar *sum;
557: #endif
559: #if PetscDefined(HAVE_PRAGMA_DISJOINT)
560: #pragma disjoint(*x, *y, *aval)
561: #endif
563: PetscFunctionBegin;
564: if (!a->nz) {
565: PetscCall(VecCopy(yy, zz));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
568: PetscCall(VecGetArrayRead(xx, &x));
569: PetscCall(VecGetArrayPair(yy, zz, &y, &z));
570: #if PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX512F__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
571: PetscCheck(a->sliceheight == 8, PETSC_COMM_SELF, PETSC_ERR_SUP, "The kernel requires a slice height of 8, but the input matrix has a slice height of %" PetscInt_FMT, a->sliceheight);
572: for (i = 0; i < totalslices; i++) { /* loop over slices */
573: PetscPrefetchBlock(acolidx, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
574: PetscPrefetchBlock(aval, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
576: if (i == totalslices - 1 && A->rmap->n & 0x07) { /* if last slice has padding rows */
577: mask = (__mmask8)(0xff >> (8 - (A->rmap->n & 0x07)));
578: vec_y = _mm512_maskz_loadu_pd(mask, &y[8 * i]);
579: } else {
580: vec_y = _mm512_loadu_pd(&y[8 * i]);
581: }
582: vec_y2 = _mm512_setzero_pd();
583: vec_y3 = _mm512_setzero_pd();
584: vec_y4 = _mm512_setzero_pd();
586: j = a->sliidx[i] >> 3; /* 8 bytes are read at each time, corresponding to a slice column */
587: switch ((a->sliidx[i + 1] - a->sliidx[i]) / 8 & 3) {
588: case 3:
589: AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
590: acolidx += 8;
591: aval += 8;
592: AVX512_Mult_Private(vec_idx2, vec_x2, vec_vals2, vec_y2);
593: acolidx += 8;
594: aval += 8;
595: AVX512_Mult_Private(vec_idx3, vec_x3, vec_vals3, vec_y3);
596: acolidx += 8;
597: aval += 8;
598: j += 3;
599: break;
600: case 2:
601: AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
602: acolidx += 8;
603: aval += 8;
604: AVX512_Mult_Private(vec_idx2, vec_x2, vec_vals2, vec_y2);
605: acolidx += 8;
606: aval += 8;
607: j += 2;
608: break;
609: case 1:
610: AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
611: acolidx += 8;
612: aval += 8;
613: j += 1;
614: break;
615: }
616: #pragma novector
617: for (; j < (a->sliidx[i + 1] >> 3); j += 4) {
618: AVX512_Mult_Private(vec_idx, vec_x, vec_vals, vec_y);
619: acolidx += 8;
620: aval += 8;
621: AVX512_Mult_Private(vec_idx2, vec_x2, vec_vals2, vec_y2);
622: acolidx += 8;
623: aval += 8;
624: AVX512_Mult_Private(vec_idx3, vec_x3, vec_vals3, vec_y3);
625: acolidx += 8;
626: aval += 8;
627: AVX512_Mult_Private(vec_idx4, vec_x4, vec_vals4, vec_y4);
628: acolidx += 8;
629: aval += 8;
630: }
632: vec_y = _mm512_add_pd(vec_y, vec_y2);
633: vec_y = _mm512_add_pd(vec_y, vec_y3);
634: vec_y = _mm512_add_pd(vec_y, vec_y4);
635: if (i == totalslices - 1 && A->rmap->n & 0x07) { /* if last slice has padding rows */
636: _mm512_mask_storeu_pd(&z[8 * i], mask, vec_y);
637: } else {
638: _mm512_storeu_pd(&z[8 * i], vec_y);
639: }
640: }
641: #elif PetscDefined(HAVE_IMMINTRIN_H) && defined(__AVX__) && PetscDefined(USE_REAL_DOUBLE) && !PetscDefined(USE_COMPLEX) && !PetscDefined(USE_64BIT_INDICES)
642: PetscCheck(a->sliceheight == 8, PETSC_COMM_SELF, PETSC_ERR_SUP, "The kernel requires a slice height of 8, but the input matrix has a slice height of %" PetscInt_FMT, a->sliceheight);
643: for (i = 0; i < totalslices; i++) { /* loop over full slices */
644: PetscPrefetchBlock(acolidx, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
645: PetscPrefetchBlock(aval, a->sliidx[i + 1] - a->sliidx[i], 0, PETSC_PREFETCH_HINT_T0);
647: /* last slice may have padding rows. Don't use vectorization. */
648: if (i == totalslices - 1 && (A->rmap->n & 0x07)) {
649: for (r = 0; r < (A->rmap->n & 0x07); ++r) {
650: row = 8 * i + r;
651: yval = (MatScalar)0.0;
652: nnz_in_row = a->rlen[row];
653: for (j = 0; j < nnz_in_row; ++j) yval += aval[8 * j + r] * x[acolidx[8 * j + r]];
654: z[row] = y[row] + yval;
655: }
656: break;
657: }
659: vec_y = _mm256_loadu_pd(y + 8 * i);
660: vec_y2 = _mm256_loadu_pd(y + 8 * i + 4);
662: /* Process slice of height 8 (512 bits) via two subslices of height 4 (256 bits) via AVX */
663: for (j = a->sliidx[i]; j < a->sliidx[i + 1]; j += 8) {
664: vec_vals = _mm256_loadu_pd(aval);
665: vec_x_tmp = _mm_setzero_pd();
666: vec_x_tmp = _mm_loadl_pd(vec_x_tmp, x + *acolidx++);
667: vec_x_tmp = _mm_loadh_pd(vec_x_tmp, x + *acolidx++);
668: vec_x = _mm256_setzero_pd();
669: vec_x = _mm256_insertf128_pd(vec_x, vec_x_tmp, 0);
670: vec_x_tmp = _mm_loadl_pd(vec_x_tmp, x + *acolidx++);
671: vec_x_tmp = _mm_loadh_pd(vec_x_tmp, x + *acolidx++);
672: vec_x = _mm256_insertf128_pd(vec_x, vec_x_tmp, 1);
673: vec_y = _mm256_add_pd(_mm256_mul_pd(vec_x, vec_vals), vec_y);
674: aval += 4;
676: vec_vals = _mm256_loadu_pd(aval);
677: vec_x_tmp = _mm_loadl_pd(vec_x_tmp, x + *acolidx++);
678: vec_x_tmp = _mm_loadh_pd(vec_x_tmp, x + *acolidx++);
679: vec_x = _mm256_insertf128_pd(vec_x, vec_x_tmp, 0);
680: vec_x_tmp = _mm_loadl_pd(vec_x_tmp, x + *acolidx++);
681: vec_x_tmp = _mm_loadh_pd(vec_x_tmp, x + *acolidx++);
682: vec_x = _mm256_insertf128_pd(vec_x, vec_x_tmp, 1);
683: vec_y2 = _mm256_add_pd(_mm256_mul_pd(vec_x, vec_vals), vec_y2);
684: aval += 4;
685: }
687: _mm256_storeu_pd(z + i * 8, vec_y);
688: _mm256_storeu_pd(z + i * 8 + 4, vec_y2);
689: }
690: #else
691: PetscCall(PetscMalloc1(sliceheight, &sum));
692: for (i = 0; i < totalslices; i++) { /* loop over slices */
693: for (j = 0; j < sliceheight; j++) {
694: sum[j] = 0.0;
695: for (k = a->sliidx[i] + j; k < a->sliidx[i + 1]; k += sliceheight) sum[j] += aval[k] * x[acolidx[k]];
696: }
697: if (i == totalslices - 1 && (A->rmap->n % sliceheight)) {
698: for (j = 0; j < (A->rmap->n % sliceheight); j++) z[sliceheight * i + j] = y[sliceheight * i + j] + sum[j];
699: } else {
700: for (j = 0; j < sliceheight; j++) z[sliceheight * i + j] = y[sliceheight * i + j] + sum[j];
701: }
702: }
703: PetscCall(PetscFree(sum));
704: #endif
706: PetscCall(PetscLogFlops(2.0 * a->nz));
707: PetscCall(VecRestoreArrayRead(xx, &x));
708: PetscCall(VecRestoreArrayPair(yy, zz, &y, &z));
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: PetscErrorCode MatMultTransposeAdd_SeqSELL(Mat A, Vec xx, Vec zz, Vec yy)
713: {
714: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
715: PetscScalar *y;
716: const PetscScalar *x;
717: const MatScalar *aval = a->val;
718: const PetscInt *acolidx = a->colidx;
719: PetscInt i, j, r, row, nnz_in_row, totalslices = a->totalslices, sliceheight = a->sliceheight;
721: #if PetscDefined(HAVE_PRAGMA_DISJOINT)
722: #pragma disjoint(*x, *y, *aval)
723: #endif
725: PetscFunctionBegin;
726: if (A->symmetric == PETSC_BOOL3_TRUE) {
727: PetscCall(MatMultAdd_SeqSELL(A, xx, zz, yy));
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
730: if (zz != yy) PetscCall(VecCopy(zz, yy));
732: if (a->nz) {
733: PetscCall(VecGetArrayRead(xx, &x));
734: PetscCall(VecGetArray(yy, &y));
735: for (i = 0; i < a->totalslices; i++) { /* loop over slices */
736: if (i == totalslices - 1 && (A->rmap->n % sliceheight)) {
737: for (r = 0; r < (A->rmap->n % sliceheight); ++r) {
738: row = sliceheight * i + r;
739: nnz_in_row = a->rlen[row];
740: for (j = 0; j < nnz_in_row; ++j) y[acolidx[sliceheight * j + r]] += aval[sliceheight * j + r] * x[row];
741: }
742: break;
743: }
744: for (r = 0; r < sliceheight; ++r)
745: for (j = a->sliidx[i] + r; j < a->sliidx[i + 1]; j += sliceheight) y[acolidx[j]] += aval[j] * x[sliceheight * i + r];
746: }
747: PetscCall(PetscLogFlops(2.0 * a->nz));
748: PetscCall(VecRestoreArrayRead(xx, &x));
749: PetscCall(VecRestoreArray(yy, &y));
750: }
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: PetscErrorCode MatMultTranspose_SeqSELL(Mat A, Vec xx, Vec yy)
755: {
756: PetscFunctionBegin;
757: if (A->symmetric == PETSC_BOOL3_TRUE) {
758: PetscCall(MatMult_SeqSELL(A, xx, yy));
759: } else {
760: PetscCall(VecSet(yy, 0.0));
761: PetscCall(MatMultTransposeAdd_SeqSELL(A, xx, yy, yy));
762: }
763: PetscFunctionReturn(PETSC_SUCCESS);
764: }
766: static PetscErrorCode MatGetDiagonalMarkers_SeqSELL(Mat A, const PetscInt **diag, PetscBool *diagDense)
767: {
768: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
770: PetscFunctionBegin;
771: if (A->factortype != MAT_FACTOR_NONE) {
772: PetscAssertPointer(diag, 2);
773: PetscCheck(!diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot check for dense diagonal with factored matrices");
774: *diag = a->diag;
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
777: PetscCheck(diag || diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "At least one of diag or diagDense must be requested");
778: if (a->diagNonzeroState != A->nonzerostate || (diag && !a->diag)) {
779: const PetscInt m = A->rmap->n;
780: PetscInt shift;
782: if (!diag && !a->diag) {
783: a->diagDense = PETSC_TRUE;
784: for (PetscInt i = 0; i < m; i++) {
785: PetscBool found = PETSC_FALSE;
787: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight; /* starting index of the row i */
788: for (PetscInt j = 0; j < a->rlen[i]; j++) {
789: if (a->colidx[shift + a->sliceheight * j] == i) {
790: a->diag[i] = shift + a->sliceheight * j;
791: found = PETSC_TRUE;
792: break;
793: }
794: }
795: if (!found) {
796: a->diagDense = PETSC_FALSE;
797: *diagDense = a->diagDense;
798: a->diagNonzeroState = A->nonzerostate;
799: PetscFunctionReturn(PETSC_SUCCESS);
800: }
801: }
802: } else {
803: if (!a->diag) PetscCall(PetscMalloc1(m, &a->diag));
804: a->diagDense = PETSC_TRUE;
805: for (PetscInt i = 0; i < m; i++) {
806: PetscBool found = PETSC_FALSE;
808: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight; /* starting index of the row i */
809: a->diag[i] = -1;
810: for (PetscInt j = 0; j < a->rlen[i]; j++) {
811: if (a->colidx[shift + a->sliceheight * j] == i) {
812: a->diag[i] = shift + a->sliceheight * j;
813: found = PETSC_TRUE;
814: break;
815: }
816: }
817: if (!found) a->diagDense = PETSC_FALSE;
818: }
819: }
820: a->diagNonzeroState = A->nonzerostate;
821: }
822: if (diag) *diag = a->diag;
823: if (diagDense) *diagDense = a->diagDense;
824: PetscFunctionReturn(PETSC_SUCCESS);
825: }
827: /*
828: Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways
829: */
830: static PetscErrorCode MatInvertDiagonalForSOR_SeqSELL(Mat A, PetscScalar omega, PetscScalar fshift)
831: {
832: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
833: PetscInt i, m = A->rmap->n;
834: MatScalar *val = a->val;
835: PetscScalar *idiag, *mdiag;
836: const PetscInt *diag;
837: PetscBool diagDense;
839: PetscFunctionBegin;
840: if (a->idiagState == ((PetscObject)A)->state && a->omega == omega && a->fshift == fshift) PetscFunctionReturn(PETSC_SUCCESS);
841: PetscCall(MatGetDiagonalMarkers_SeqSELL(A, &diag, &diagDense));
842: PetscCheck(diagDense, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix must have all diagonal locations to invert them");
844: if (!a->idiag) {
845: PetscCall(PetscMalloc3(m, &a->idiag, m, &a->mdiag, m, &a->ssor_work));
846: val = a->val;
847: }
848: mdiag = a->mdiag;
849: idiag = a->idiag;
851: if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) {
852: for (i = 0; i < m; i++) {
853: mdiag[i] = val[diag[i]];
854: if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */
855: PetscCheck(PetscRealPart(fshift), PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Zero diagonal on row %" PetscInt_FMT, i);
856: PetscCall(PetscInfo(A, "Zero diagonal on row %" PetscInt_FMT "\n", i));
857: A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT;
858: A->factorerror_zeropivot_value = 0.0;
859: A->factorerror_zeropivot_row = i;
860: }
861: idiag[i] = 1.0 / val[diag[i]];
862: }
863: PetscCall(PetscLogFlops(m));
864: } else {
865: for (i = 0; i < m; i++) {
866: mdiag[i] = val[diag[i]];
867: idiag[i] = omega / (fshift + val[diag[i]]);
868: }
869: PetscCall(PetscLogFlops(2.0 * m));
870: }
871: a->idiagState = ((PetscObject)A)->state;
872: a->omega = omega;
873: a->fshift = fshift;
874: PetscFunctionReturn(PETSC_SUCCESS);
875: }
877: PetscErrorCode MatZeroEntries_SeqSELL(Mat A)
878: {
879: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
881: PetscFunctionBegin;
882: PetscCall(PetscArrayzero(a->val, a->sliidx[a->totalslices]));
883: PetscFunctionReturn(PETSC_SUCCESS);
884: }
886: PetscErrorCode MatDestroy_SeqSELL(Mat A)
887: {
888: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
890: PetscFunctionBegin;
891: PetscCall(PetscLogObjectState((PetscObject)A, "Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT, A->rmap->n, A->cmap->n, a->nz));
892: PetscCall(MatSeqXSELLFreeSELL(A, &a->val, &a->colidx));
893: PetscCall(ISDestroy(&a->row));
894: PetscCall(ISDestroy(&a->col));
895: PetscCall(PetscFree(a->diag));
896: PetscCall(PetscFree(a->rlen));
897: PetscCall(PetscFree(a->sliidx));
898: PetscCall(PetscFree3(a->idiag, a->mdiag, a->ssor_work));
899: PetscCall(PetscFree(a->solve_work));
900: PetscCall(ISDestroy(&a->icol));
901: PetscCall(PetscFree(a->saved_values));
902: PetscCall(PetscFree2(a->getrowcols, a->getrowvals));
903: #if PetscDefined(HAVE_CUPM)
904: PetscCall(PetscFree(a->chunk_slice_map));
905: #endif
906: PetscCall(PetscFree(A->data));
908: PetscCall(PetscObjectChangeTypeName((PetscObject)A, NULL));
909: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatStoreValues_C", NULL));
910: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatRetrieveValues_C", NULL));
911: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSELLSetPreallocation_C", NULL));
912: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSELLGetArray_C", NULL));
913: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSELLRestoreArray_C", NULL));
914: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsell_seqaij_C", NULL));
915: #if PetscDefined(HAVE_CUDA)
916: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsell_seqsellcuda_C", NULL));
917: #endif
918: #if PetscDefined(HAVE_HIP)
919: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_seqsell_seqsellhip_C", NULL));
920: #endif
921: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSELLGetFillRatio_C", NULL));
922: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSELLGetMaxSliceWidth_C", NULL));
923: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSELLGetAvgSliceWidth_C", NULL));
924: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSELLGetVarSliceSize_C", NULL));
925: PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSeqSELLSetSliceHeight_C", NULL));
926: PetscFunctionReturn(PETSC_SUCCESS);
927: }
929: PetscErrorCode MatSetOption_SeqSELL(Mat A, MatOption op, PetscBool flg)
930: {
931: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
933: PetscFunctionBegin;
934: switch (op) {
935: case MAT_ROW_ORIENTED:
936: a->roworiented = flg;
937: break;
938: case MAT_KEEP_NONZERO_PATTERN:
939: a->keepnonzeropattern = flg;
940: break;
941: case MAT_NEW_NONZERO_LOCATIONS:
942: a->nonew = (flg ? 0 : 1);
943: break;
944: case MAT_NEW_NONZERO_LOCATION_ERR:
945: a->nonew = (flg ? -1 : 0);
946: break;
947: case MAT_NEW_NONZERO_ALLOCATION_ERR:
948: a->nonew = (flg ? -2 : 0);
949: break;
950: case MAT_UNUSED_NONZERO_LOCATION_ERR:
951: a->nounused = (flg ? -1 : 0);
952: break;
953: case MAT_IGNORE_ZERO_ENTRIES:
954: a->ignorezeroentries = flg;
955: break;
956: default:
957: break;
958: }
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: PetscErrorCode MatGetDiagonal_SeqSELL(Mat A, Vec v)
963: {
964: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
965: PetscInt i, j, n, shift;
966: PetscScalar *x, zero = 0.0;
968: PetscFunctionBegin;
969: PetscCall(VecGetLocalSize(v, &n));
970: PetscCheck(n == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Nonconforming matrix and vector");
972: if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) {
973: const PetscInt *diag;
975: PetscCall(MatGetDiagonalMarkers_SeqSELL(A, &diag, NULL));
976: PetscCall(VecGetArrayWrite(v, &x));
977: for (i = 0; i < n; i++) x[i] = 1.0 / a->val[diag[i]];
978: PetscCall(VecRestoreArrayWrite(v, &x));
979: PetscFunctionReturn(PETSC_SUCCESS);
980: }
982: PetscCall(VecSet(v, zero));
983: PetscCall(VecGetArray(v, &x));
984: for (i = 0; i < n; i++) { /* loop over rows */
985: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight; /* starting index of the row i */
986: x[i] = 0;
987: for (j = 0; j < a->rlen[i]; j++) {
988: if (a->colidx[shift + a->sliceheight * j] == i) {
989: x[i] = a->val[shift + a->sliceheight * j];
990: break;
991: }
992: }
993: }
994: PetscCall(VecRestoreArray(v, &x));
995: PetscFunctionReturn(PETSC_SUCCESS);
996: }
998: PetscErrorCode MatDiagonalScale_SeqSELL(Mat A, Vec ll, Vec rr)
999: {
1000: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1001: const PetscScalar *l, *r;
1002: PetscInt i, j, m, n, row;
1004: PetscFunctionBegin;
1005: if (ll) {
1006: /* The local size is used so that VecMPI can be passed to this routine
1007: by MatDiagonalScale_MPISELL */
1008: PetscCall(VecGetLocalSize(ll, &m));
1009: PetscCheck(m == A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Left scaling vector wrong length");
1010: PetscCall(VecGetArrayRead(ll, &l));
1011: for (i = 0; i < a->totalslices; i++) { /* loop over slices */
1012: if (i == a->totalslices - 1 && (A->rmap->n % a->sliceheight)) { /* if last slice has padding rows */
1013: for (j = a->sliidx[i], row = 0; j < a->sliidx[i + 1]; j++, row = (row + 1) % a->sliceheight) {
1014: if (row < (A->rmap->n % a->sliceheight)) a->val[j] *= l[a->sliceheight * i + row];
1015: }
1016: } else {
1017: for (j = a->sliidx[i], row = 0; j < a->sliidx[i + 1]; j++, row = (row + 1) % a->sliceheight) a->val[j] *= l[a->sliceheight * i + row];
1018: }
1019: }
1020: PetscCall(VecRestoreArrayRead(ll, &l));
1021: PetscCall(PetscLogFlops(a->nz));
1022: }
1023: if (rr) {
1024: PetscCall(VecGetLocalSize(rr, &n));
1025: PetscCheck(n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Right scaling vector wrong length");
1026: PetscCall(VecGetArrayRead(rr, &r));
1027: for (i = 0; i < a->totalslices; i++) { /* loop over slices */
1028: if (i == a->totalslices - 1 && (A->rmap->n % a->sliceheight)) { /* if last slice has padding rows */
1029: for (j = a->sliidx[i], row = 0; j < a->sliidx[i + 1]; j++, row = ((row + 1) % a->sliceheight)) {
1030: if (row < (A->rmap->n % a->sliceheight)) a->val[j] *= r[a->colidx[j]];
1031: }
1032: } else {
1033: for (j = a->sliidx[i]; j < a->sliidx[i + 1]; j++) a->val[j] *= r[a->colidx[j]];
1034: }
1035: }
1036: PetscCall(VecRestoreArrayRead(rr, &r));
1037: PetscCall(PetscLogFlops(a->nz));
1038: }
1039: #if PetscDefined(HAVE_CUPM)
1040: if (A->offloadmask != PETSC_OFFLOAD_UNALLOCATED) A->offloadmask = PETSC_OFFLOAD_CPU;
1041: #endif
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: PetscErrorCode MatGetValues_SeqSELL(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], PetscScalar v[])
1046: {
1047: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1048: PetscInt *cp, i, k, low, high, t, row, col, l;
1049: PetscInt shift;
1050: MatScalar *vp;
1051: PetscBool roworiented = a->roworiented;
1052: PetscScalar *value;
1054: PetscFunctionBegin;
1055: for (k = 0; k < m; k++) { /* loop over requested rows */
1056: row = im[k];
1057: if (row < 0) continue;
1058: PetscCheck(row < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, A->rmap->n - 1);
1059: shift = a->sliidx[row / a->sliceheight] + (row % a->sliceheight); /* starting index of the row */
1060: cp = a->colidx + shift; /* pointer to the row */
1061: vp = a->val + shift; /* pointer to the row */
1062: for (l = 0; l < n; l++) { /* loop over requested columns */
1063: col = in[l];
1064: if (col < 0) continue;
1065: PetscCheck(col < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Column too large: row %" PetscInt_FMT " max %" PetscInt_FMT, col, A->cmap->n - 1);
1066: value = roworiented ? &v[l + k * n] : &v[k + l * m];
1067: high = a->rlen[row];
1068: low = 0; /* assume unsorted */
1069: while (high - low > 5) {
1070: t = (low + high) / 2;
1071: if (*(cp + a->sliceheight * t) > col) high = t;
1072: else low = t;
1073: }
1074: for (i = low; i < high; i++) {
1075: if (*(cp + a->sliceheight * i) > col) break;
1076: if (*(cp + a->sliceheight * i) == col) {
1077: *value = *(vp + a->sliceheight * i);
1078: goto finished;
1079: }
1080: }
1081: *value = 0.0;
1082: finished:;
1083: }
1084: }
1085: PetscFunctionReturn(PETSC_SUCCESS);
1086: }
1088: static PetscErrorCode MatView_SeqSELL_ASCII(Mat A, PetscViewer viewer)
1089: {
1090: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1091: PetscInt i, j, m = A->rmap->n, shift;
1092: const char *name;
1093: PetscViewerFormat format;
1095: PetscFunctionBegin;
1096: PetscCall(PetscViewerGetFormat(viewer, &format));
1097: if (format == PETSC_VIEWER_ASCII_MATLAB) {
1098: PetscInt nofinalvalue = 0;
1099: /*
1100: if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) nofinalvalue = 1;
1101: */
1102: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1103: PetscCall(PetscViewerASCIIPrintf(viewer, "%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n", m, A->cmap->n));
1104: PetscCall(PetscViewerASCIIPrintf(viewer, "%% Nonzeros = %" PetscInt_FMT " \n", a->nz));
1105: PetscCall(PetscViewerASCIIPrintf(viewer, "zzz = zeros(%" PetscInt_FMT ",%d);\n", a->nz + nofinalvalue, PetscDefined(USE_COMPLEX) ? 4 : 3));
1106: PetscCall(PetscViewerASCIIPrintf(viewer, "zzz = [\n"));
1108: for (i = 0; i < m; i++) {
1109: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1110: for (j = 0; j < a->rlen[i]; j++) {
1111: #if PetscDefined(USE_COMPLEX)
1112: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n", i + 1, a->colidx[shift + a->sliceheight * j] + 1, (double)PetscRealPart(a->val[shift + a->sliceheight * j]), (double)PetscImaginaryPart(a->val[shift + a->sliceheight * j])));
1113: #else
1114: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n", i + 1, a->colidx[shift + a->sliceheight * j] + 1, (double)a->val[shift + a->sliceheight * j]));
1115: #endif
1116: }
1117: }
1118: /*
1119: if (nofinalvalue) {
1120: #if PetscDefined(USE_COMPLEX)
1121: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",m,A->cmap->n,0.,0.));
1122: #else
1123: PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",m,A->cmap->n,0.0));
1124: #endif
1125: }
1126: */
1127: PetscCall(PetscObjectGetName((PetscObject)A, &name));
1128: PetscCall(PetscViewerASCIIPrintf(viewer, "];\n %s = spconvert(zzz);\n", name));
1129: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1130: } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) {
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: } else if (format == PETSC_VIEWER_ASCII_COMMON) {
1133: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1134: for (i = 0; i < m; i++) {
1135: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1136: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1137: for (j = 0; j < a->rlen[i]; j++) {
1138: #if PetscDefined(USE_COMPLEX)
1139: if (PetscImaginaryPart(a->val[shift + a->sliceheight * j]) > 0.0 && PetscRealPart(a->val[shift + a->sliceheight * j]) != 0.0) {
1140: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->colidx[shift + a->sliceheight * j], (double)PetscRealPart(a->val[shift + a->sliceheight * j]), (double)PetscImaginaryPart(a->val[shift + a->sliceheight * j])));
1141: } else if (PetscImaginaryPart(a->val[shift + a->sliceheight * j]) < 0.0 && PetscRealPart(a->val[shift + a->sliceheight * j]) != 0.0) {
1142: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->colidx[shift + a->sliceheight * j], (double)PetscRealPart(a->val[shift + a->sliceheight * j]), (double)-PetscImaginaryPart(a->val[shift + a->sliceheight * j])));
1143: } else if (PetscRealPart(a->val[shift + a->sliceheight * j]) != 0.0) {
1144: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->colidx[shift + a->sliceheight * j], (double)PetscRealPart(a->val[shift + a->sliceheight * j])));
1145: }
1146: #else
1147: if (a->val[shift + a->sliceheight * j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->colidx[shift + a->sliceheight * j], (double)a->val[shift + a->sliceheight * j]));
1148: #endif
1149: }
1150: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1151: }
1152: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1153: } else if (format == PETSC_VIEWER_ASCII_DENSE) {
1154: PetscInt cnt = 0, jcnt;
1155: PetscScalar value;
1156: PetscBool realonly = PETSC_TRUE;
1158: if (PetscDefined(USE_COMPLEX)) {
1159: for (i = 0; i < a->sliidx[a->totalslices]; i++) {
1160: if (PetscImaginaryPart(a->val[i]) != 0.0) {
1161: realonly = PETSC_FALSE;
1162: break;
1163: }
1164: }
1165: }
1167: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1168: for (i = 0; i < m; i++) {
1169: jcnt = 0;
1170: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1171: for (j = 0; j < A->cmap->n; j++) {
1172: if (jcnt < a->rlen[i] && j == a->colidx[shift + a->sliceheight * j]) {
1173: value = a->val[cnt++];
1174: jcnt++;
1175: } else {
1176: value = 0.0;
1177: }
1178: if (!PetscDefined(USE_COMPLEX) || realonly) {
1179: PetscCall(PetscViewerASCIIPrintf(viewer, " %7.5e ", (double)PetscRealPart(value)));
1180: } else {
1181: PetscCall(PetscViewerASCIIPrintf(viewer, " %7.5e+%7.5e i ", (double)PetscRealPart(value), (double)PetscImaginaryPart(value)));
1182: }
1183: }
1184: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1185: }
1186: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1187: } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) {
1188: PetscInt fshift = 1;
1189: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1190: PetscCall(PetscViewerASCIIPrintf(viewer, "%%%%MatrixMarket matrix coordinate %s general\n", PetscDefined(USE_COMPLEX) ? "complex" : "real"));
1191: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, A->cmap->n, a->nz));
1192: for (i = 0; i < m; i++) {
1193: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1194: for (j = 0; j < a->rlen[i]; j++) {
1195: #if PetscDefined(USE_COMPLEX)
1196: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %g %g\n", i + fshift, a->colidx[shift + a->sliceheight * j] + fshift, (double)PetscRealPart(a->val[shift + a->sliceheight * j]), (double)PetscImaginaryPart(a->val[shift + a->sliceheight * j])));
1197: #else
1198: PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " %" PetscInt_FMT " %g\n", i + fshift, a->colidx[shift + a->sliceheight * j] + fshift, (double)a->val[shift + a->sliceheight * j]));
1199: #endif
1200: }
1201: }
1202: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1203: } else if (format == PETSC_VIEWER_NATIVE) {
1204: for (i = 0; i < a->totalslices; i++) { /* loop over slices */
1205: PetscInt row;
1206: PetscCall(PetscViewerASCIIPrintf(viewer, "slice %" PetscInt_FMT ": %" PetscInt_FMT " %" PetscInt_FMT "\n", i, a->sliidx[i], a->sliidx[i + 1]));
1207: for (j = a->sliidx[i], row = 0; j < a->sliidx[i + 1]; j++, row = (row + 1) % a->sliceheight) {
1208: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[j]) > 0.0) {
1209: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " %g + %g i\n", a->sliceheight * i + row, a->colidx[j], (double)PetscRealPart(a->val[j]), (double)PetscImaginaryPart(a->val[j])));
1210: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[j]) < 0.0) {
1211: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " %g - %g i\n", a->sliceheight * i + row, a->colidx[j], (double)PetscRealPart(a->val[j]), -(double)PetscImaginaryPart(a->val[j])));
1212: } else {
1213: PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT " %" PetscInt_FMT " %g\n", a->sliceheight * i + row, a->colidx[j], (double)PetscRealPart(a->val[j])));
1214: }
1215: }
1216: }
1217: } else {
1218: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
1219: if (A->factortype) {
1220: for (i = 0; i < m; i++) {
1221: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1222: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1223: /* L part */
1224: for (j = shift; j < a->diag[i]; j += a->sliceheight) {
1225: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[shift + a->sliceheight * j]) > 0.0) {
1226: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->colidx[j], (double)PetscRealPart(a->val[j]), (double)PetscImaginaryPart(a->val[j])));
1227: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[shift + a->sliceheight * j]) < 0.0) {
1228: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->colidx[j], (double)PetscRealPart(a->val[j]), (double)(-PetscImaginaryPart(a->val[j]))));
1229: } else {
1230: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->colidx[j], (double)PetscRealPart(a->val[j])));
1231: }
1232: }
1233: /* diagonal */
1234: j = a->diag[i];
1235: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[j]) > 0.0) {
1236: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->colidx[j], (double)PetscRealPart(1.0 / a->val[j]), (double)PetscImaginaryPart(1.0 / a->val[j])));
1237: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[j]) < 0.0) {
1238: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->colidx[j], (double)PetscRealPart(1.0 / a->val[j]), (double)(-PetscImaginaryPart(1.0 / a->val[j]))));
1239: } else {
1240: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->colidx[j], (double)PetscRealPart(1.0 / a->val[j])));
1241: }
1243: /* U part */
1244: for (j = a->diag[i] + 1; j < shift + a->sliceheight * a->rlen[i]; j += a->sliceheight) {
1245: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[j]) > 0.0) {
1246: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->colidx[j], (double)PetscRealPart(a->val[j]), (double)PetscImaginaryPart(a->val[j])));
1247: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[j]) < 0.0) {
1248: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->colidx[j], (double)PetscRealPart(a->val[j]), (double)(-PetscImaginaryPart(a->val[j]))));
1249: } else {
1250: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->colidx[j], (double)PetscRealPart(a->val[j])));
1251: }
1252: }
1253: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1254: }
1255: } else {
1256: for (i = 0; i < m; i++) {
1257: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1258: PetscCall(PetscViewerASCIIPrintf(viewer, "row %" PetscInt_FMT ":", i));
1259: for (j = 0; j < a->rlen[i]; j++) {
1260: if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[j]) > 0.0) {
1261: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g + %g i)", a->colidx[shift + a->sliceheight * j], (double)PetscRealPart(a->val[shift + a->sliceheight * j]), (double)PetscImaginaryPart(a->val[shift + a->sliceheight * j])));
1262: } else if (PetscDefined(USE_COMPLEX) && PetscImaginaryPart(a->val[j]) < 0.0) {
1263: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g - %g i)", a->colidx[shift + a->sliceheight * j], (double)PetscRealPart(a->val[shift + a->sliceheight * j]), (double)-PetscImaginaryPart(a->val[shift + a->sliceheight * j])));
1264: } else {
1265: PetscCall(PetscViewerASCIIPrintf(viewer, " (%" PetscInt_FMT ", %g) ", a->colidx[shift + a->sliceheight * j], (double)PetscRealPart(a->val[shift + a->sliceheight * j])));
1266: }
1267: }
1268: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
1269: }
1270: }
1271: PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
1272: }
1273: PetscCall(PetscViewerFlush(viewer));
1274: PetscFunctionReturn(PETSC_SUCCESS);
1275: }
1277: #include <petscdraw.h>
1278: #if defined(__GNUC__) && !defined(__clang__)
1279: #pragma GCC diagnostic push
1280: #pragma GCC diagnostic ignored "-Wclobbered"
1281: #endif
1282: static PetscErrorCode MatView_SeqSELL_Draw_Zoom(PetscDraw draw, void *Aa)
1283: {
1284: Mat A = (Mat)Aa;
1285: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1286: PetscInt i, j, m = A->rmap->n, shift;
1287: int color;
1288: PetscReal xl, yl, xr, yr, x_l, x_r, y_l, y_r;
1289: PetscViewer viewer;
1290: PetscViewerFormat format;
1292: PetscFunctionBegin;
1293: PetscCall(PetscObjectQuery((PetscObject)A, "Zoomviewer", (PetscObject *)&viewer));
1294: PetscCall(PetscViewerGetFormat(viewer, &format));
1295: PetscCall(PetscDrawGetCoordinates(draw, &xl, &yl, &xr, &yr));
1297: /* loop over matrix elements drawing boxes */
1299: if (format != PETSC_VIEWER_DRAW_CONTOUR) {
1300: PetscDrawCollectiveBegin(draw);
1301: /* Blue for negative, Cyan for zero and Red for positive */
1302: color = PETSC_DRAW_BLUE;
1303: for (i = 0; i < m; i++) {
1304: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight; /* starting index of the row i */
1305: y_l = m - i - 1.0;
1306: y_r = y_l + 1.0;
1307: for (j = 0; j < a->rlen[i]; j++) {
1308: x_l = a->colidx[shift + a->sliceheight * j];
1309: x_r = x_l + 1.0;
1310: if (PetscRealPart(a->val[shift + a->sliceheight * j]) >= 0.) continue;
1311: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1312: }
1313: }
1314: color = PETSC_DRAW_CYAN;
1315: for (i = 0; i < m; i++) {
1316: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1317: y_l = m - i - 1.0;
1318: y_r = y_l + 1.0;
1319: for (j = 0; j < a->rlen[i]; j++) {
1320: x_l = a->colidx[shift + a->sliceheight * j];
1321: x_r = x_l + 1.0;
1322: if (a->val[shift + a->sliceheight * j] != 0.) continue;
1323: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1324: }
1325: }
1326: color = PETSC_DRAW_RED;
1327: for (i = 0; i < m; i++) {
1328: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1329: y_l = m - i - 1.0;
1330: y_r = y_l + 1.0;
1331: for (j = 0; j < a->rlen[i]; j++) {
1332: x_l = a->colidx[shift + a->sliceheight * j];
1333: x_r = x_l + 1.0;
1334: if (PetscRealPart(a->val[shift + a->sliceheight * j]) <= 0.) continue;
1335: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1336: }
1337: }
1338: PetscDrawCollectiveEnd(draw);
1339: } else {
1340: /* use contour shading to indicate magnitude of values */
1341: /* first determine max of all nonzero values */
1342: PetscReal minv = 0.0, maxv = 0.0;
1343: PetscInt count = 0;
1344: PetscDraw popup;
1345: for (i = 0; i < a->sliidx[a->totalslices]; i++) {
1346: if (PetscAbsScalar(a->val[i]) > maxv) maxv = PetscAbsScalar(a->val[i]);
1347: }
1348: if (minv >= maxv) maxv = minv + PETSC_SMALL;
1349: PetscCall(PetscDrawGetPopup(draw, &popup));
1350: PetscCall(PetscDrawScalePopup(popup, minv, maxv));
1352: PetscDrawCollectiveBegin(draw);
1353: for (i = 0; i < m; i++) {
1354: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight;
1355: y_l = m - i - 1.0;
1356: y_r = y_l + 1.0;
1357: for (j = 0; j < a->rlen[i]; j++) {
1358: x_l = a->colidx[shift + a->sliceheight * j];
1359: x_r = x_l + 1.0;
1360: color = PetscDrawRealToColor(PetscAbsScalar(a->val[count]), minv, maxv);
1361: PetscCall(PetscDrawRectangle(draw, x_l, y_l, x_r, y_r, color, color, color, color));
1362: count++;
1363: }
1364: }
1365: PetscDrawCollectiveEnd(draw);
1366: }
1367: PetscFunctionReturn(PETSC_SUCCESS);
1368: }
1369: #if defined(__GNUC__) && !defined(__clang__)
1370: #pragma GCC diagnostic pop
1371: #endif
1373: #include <petscdraw.h>
1374: static PetscErrorCode MatView_SeqSELL_Draw(Mat A, PetscViewer viewer)
1375: {
1376: PetscDraw draw;
1377: PetscReal xr, yr, xl, yl, h, w;
1378: PetscBool isnull;
1380: PetscFunctionBegin;
1381: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
1382: PetscCall(PetscDrawIsNull(draw, &isnull));
1383: if (isnull) PetscFunctionReturn(PETSC_SUCCESS);
1385: xr = A->cmap->n;
1386: yr = A->rmap->n;
1387: h = yr / 10.0;
1388: w = xr / 10.0;
1389: xr += w;
1390: yr += h;
1391: xl = -w;
1392: yl = -h;
1393: PetscCall(PetscDrawSetCoordinates(draw, xl, yl, xr, yr));
1394: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", (PetscObject)viewer));
1395: PetscCall(PetscDrawZoom(draw, MatView_SeqSELL_Draw_Zoom, A));
1396: PetscCall(PetscObjectCompose((PetscObject)A, "Zoomviewer", NULL));
1397: PetscCall(PetscDrawSave(draw));
1398: PetscFunctionReturn(PETSC_SUCCESS);
1399: }
1401: PetscErrorCode MatView_SeqSELL(Mat A, PetscViewer viewer)
1402: {
1403: PetscBool isascii, isbinary, isdraw;
1405: PetscFunctionBegin;
1406: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1407: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary));
1408: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
1409: if (isascii) {
1410: PetscCall(MatView_SeqSELL_ASCII(A, viewer));
1411: } else if (isbinary) {
1412: /* PetscCall(MatView_SeqSELL_Binary(A,viewer)); */
1413: } else if (isdraw) PetscCall(MatView_SeqSELL_Draw(A, viewer));
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: PetscErrorCode MatAssemblyEnd_SeqSELL(Mat A, MatAssemblyType mode)
1418: {
1419: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1420: PetscInt i, shift, row_in_slice, row, nrow, *cp, lastcol, j, k;
1421: MatScalar *vp;
1422: #if PetscDefined(HAVE_CUPM)
1423: PetscInt totalchunks = 0;
1424: #endif
1426: PetscFunctionBegin;
1427: if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(PETSC_SUCCESS);
1428: /* To do: compress out the unused elements */
1429: PetscCall(PetscInfo(A, "Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; storage space: %" PetscInt_FMT " allocated %" PetscInt_FMT " used (%" PetscInt_FMT " nonzeros+%" PetscInt_FMT " paddedzeros)\n", A->rmap->n, A->cmap->n, a->maxallocmat, a->sliidx[a->totalslices], a->nz, a->sliidx[a->totalslices] - a->nz));
1430: PetscCall(PetscInfo(A, "Number of mallocs during MatSetValues() is %" PetscInt_FMT "\n", a->reallocs));
1431: PetscCall(PetscInfo(A, "Maximum nonzeros in any row is %" PetscInt_FMT "\n", a->rlenmax));
1432: a->nonzerorowcnt = 0;
1433: /* Set unused slots for column indices to last valid column index. Set unused slots for values to zero. This allows for a use of unmasked intrinsics -> higher performance */
1434: for (i = 0; i < a->totalslices; ++i) {
1435: shift = a->sliidx[i]; /* starting index of the slice */
1436: cp = PetscSafePointerPlusOffset(a->colidx, shift); /* pointer to the column indices of the slice */
1437: vp = PetscSafePointerPlusOffset(a->val, shift); /* pointer to the nonzero values of the slice */
1438: for (row_in_slice = 0; row_in_slice < a->sliceheight; ++row_in_slice) { /* loop over rows in the slice */
1439: row = a->sliceheight * i + row_in_slice;
1440: nrow = a->rlen[row]; /* number of nonzeros in row */
1441: /*
1442: Search for the nearest nonzero. Normally setting the index to zero may cause extra communication.
1443: But if the entire slice are empty, it is fine to use 0 since the index will not be loaded.
1444: */
1445: lastcol = 0;
1446: if (nrow > 0) { /* nonempty row */
1447: a->nonzerorowcnt++;
1448: lastcol = cp[a->sliceheight * (nrow - 1) + row_in_slice]; /* use the index from the last nonzero at current row */
1449: } else if (!row_in_slice) { /* first row of the correct slice is empty */
1450: for (j = 1; j < a->sliceheight; j++) {
1451: if (a->rlen[a->sliceheight * i + j]) {
1452: lastcol = cp[j];
1453: break;
1454: }
1455: }
1456: } else {
1457: if (a->sliidx[i + 1] != shift) lastcol = cp[row_in_slice - 1]; /* use the index from the previous row */
1458: }
1460: for (k = nrow; k < (a->sliidx[i + 1] - shift) / a->sliceheight; ++k) {
1461: cp[a->sliceheight * k + row_in_slice] = lastcol;
1462: vp[a->sliceheight * k + row_in_slice] = (MatScalar)0;
1463: }
1464: }
1465: }
1467: A->info.mallocs += a->reallocs;
1468: a->reallocs = 0;
1470: #if PetscDefined(HAVE_CUPM)
1471: if (a->totalslices) {
1472: if (!a->chunksize) {
1473: a->chunksize = 64;
1474: while (a->chunksize < 1024 && 2 * a->chunksize <= a->sliidx[a->totalslices] / a->totalslices) a->chunksize *= 2;
1475: }
1476: totalchunks = PetscCeilInt(a->sliidx[a->totalslices], a->chunksize);
1477: }
1478: if (totalchunks != a->totalchunks) {
1479: PetscCall(PetscFree(a->chunk_slice_map));
1480: PetscCall(PetscMalloc1(totalchunks, &a->chunk_slice_map));
1481: a->totalchunks = totalchunks;
1482: }
1483: j = 0;
1484: for (i = 0; i < totalchunks; i++) {
1485: while (a->sliidx[j + 1] <= i * a->chunksize && j < a->totalslices) j++;
1486: a->chunk_slice_map[i] = j;
1487: }
1488: #endif
1489: PetscFunctionReturn(PETSC_SUCCESS);
1490: }
1492: PetscErrorCode MatGetInfo_SeqSELL(Mat A, MatInfoType flag, MatInfo *info)
1493: {
1494: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1496: PetscFunctionBegin;
1497: info->block_size = 1.0;
1498: info->nz_allocated = a->maxallocmat;
1499: info->nz_used = a->sliidx[a->totalslices]; /* include padding zeros */
1500: info->nz_unneeded = (a->maxallocmat - a->sliidx[a->totalslices]);
1501: info->assemblies = A->num_ass;
1502: info->mallocs = A->info.mallocs;
1503: info->memory = 0; /* REVIEW ME */
1504: if (A->factortype) {
1505: info->fill_ratio_given = A->info.fill_ratio_given;
1506: info->fill_ratio_needed = A->info.fill_ratio_needed;
1507: info->factor_mallocs = A->info.factor_mallocs;
1508: } else {
1509: info->fill_ratio_given = 0;
1510: info->fill_ratio_needed = 0;
1511: info->factor_mallocs = 0;
1512: }
1513: PetscFunctionReturn(PETSC_SUCCESS);
1514: }
1516: PetscErrorCode MatSetValues_SeqSELL(Mat A, PetscInt m, const PetscInt im[], PetscInt n, const PetscInt in[], const PetscScalar v[], InsertMode is)
1517: {
1518: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1519: PetscInt shift, i, k, l, low, high, t, ii, row, col, nrow;
1520: PetscInt *cp, nonew = a->nonew, lastcol = -1;
1521: MatScalar *vp, value;
1522: #if PetscDefined(HAVE_CUPM)
1523: PetscBool inserted = PETSC_FALSE;
1524: PetscInt mul = DEVICE_MEM_ALIGN / a->sliceheight;
1525: #endif
1527: PetscFunctionBegin;
1528: for (k = 0; k < m; k++) { /* loop over added rows */
1529: row = im[k];
1530: if (row < 0) continue;
1531: PetscCheck(row < A->rmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT, row, A->rmap->n - 1);
1532: shift = a->sliidx[row / a->sliceheight] + row % a->sliceheight; /* starting index of the row */
1533: cp = a->colidx + shift; /* pointer to the row */
1534: vp = a->val + shift; /* pointer to the row */
1535: nrow = a->rlen[row];
1536: low = 0;
1537: high = nrow;
1539: for (l = 0; l < n; l++) { /* loop over added columns */
1540: col = in[l];
1541: if (col < 0) continue;
1542: PetscCheck(col < A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Col too large: row %" PetscInt_FMT " max %" PetscInt_FMT, col, A->cmap->n - 1);
1543: if (a->roworiented) {
1544: value = v[l + k * n];
1545: } else {
1546: value = v[k + l * m];
1547: }
1548: if (value == 0.0 && a->ignorezeroentries && is == ADD_VALUES && row != col) continue;
1550: /* search in this row for the specified column, i indicates the column to be set */
1551: if (col <= lastcol) low = 0;
1552: else high = nrow;
1553: lastcol = col;
1554: while (high - low > 5) {
1555: t = (low + high) / 2;
1556: if (*(cp + a->sliceheight * t) > col) high = t;
1557: else low = t;
1558: }
1559: for (i = low; i < high; i++) {
1560: if (*(cp + a->sliceheight * i) > col) break;
1561: if (*(cp + a->sliceheight * i) == col) {
1562: if (is == ADD_VALUES) *(vp + a->sliceheight * i) += value;
1563: else *(vp + a->sliceheight * i) = value;
1564: #if PetscDefined(HAVE_CUPM)
1565: inserted = PETSC_TRUE;
1566: #endif
1567: low = i + 1;
1568: goto noinsert;
1569: }
1570: }
1571: if (value == 0.0 && a->ignorezeroentries && row != col) goto noinsert;
1572: if (nonew == 1) goto noinsert;
1573: PetscCheck(nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero (%" PetscInt_FMT ", %" PetscInt_FMT ") in the matrix", row, col);
1574: #if PetscDefined(HAVE_CUPM)
1575: MatSeqXSELLReallocateSELL(A, A->rmap->n, 1, nrow, a->sliidx, a->sliceheight, row / a->sliceheight, row, col, a->colidx, a->val, cp, vp, nonew, MatScalar, mul);
1576: #else
1577: /* If the current row length exceeds the slice width (e.g. nrow==slice_width), allocate a new space, otherwise do nothing */
1578: MatSeqXSELLReallocateSELL(A, A->rmap->n, 1, nrow, a->sliidx, a->sliceheight, row / a->sliceheight, row, col, a->colidx, a->val, cp, vp, nonew, MatScalar, 1);
1579: #endif
1580: /* add the new nonzero to the high position, shift the remaining elements in current row to the right by one slot */
1581: for (ii = nrow - 1; ii >= i; ii--) {
1582: *(cp + a->sliceheight * (ii + 1)) = *(cp + a->sliceheight * ii);
1583: *(vp + a->sliceheight * (ii + 1)) = *(vp + a->sliceheight * ii);
1584: }
1585: a->rlen[row]++;
1586: *(cp + a->sliceheight * i) = col;
1587: *(vp + a->sliceheight * i) = value;
1588: a->nz++;
1589: #if PetscDefined(HAVE_CUPM)
1590: inserted = PETSC_TRUE;
1591: #endif
1592: low = i + 1;
1593: high++;
1594: nrow++;
1595: noinsert:;
1596: }
1597: a->rlen[row] = nrow;
1598: }
1599: #if PetscDefined(HAVE_CUPM)
1600: if (A->offloadmask != PETSC_OFFLOAD_UNALLOCATED && inserted) A->offloadmask = PETSC_OFFLOAD_CPU;
1601: #endif
1602: PetscFunctionReturn(PETSC_SUCCESS);
1603: }
1605: PetscErrorCode MatCopy_SeqSELL(Mat A, Mat B, MatStructure str)
1606: {
1607: PetscFunctionBegin;
1608: /* If the two matrices have the same copy implementation, use fast copy. */
1609: if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) {
1610: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1611: Mat_SeqSELL *b = (Mat_SeqSELL *)B->data;
1613: PetscCheck(a->sliidx[a->totalslices] == b->sliidx[b->totalslices], PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number of nonzeros in two matrices are different");
1614: PetscCall(PetscArraycpy(b->val, a->val, a->sliidx[a->totalslices]));
1615: #if PetscDefined(HAVE_CUPM)
1616: if (B->offloadmask != PETSC_OFFLOAD_UNALLOCATED) B->offloadmask = PETSC_OFFLOAD_CPU;
1617: #endif
1618: } else {
1619: PetscCall(MatCopy_Basic(A, B, str));
1620: }
1621: PetscFunctionReturn(PETSC_SUCCESS);
1622: }
1624: PetscErrorCode MatSetUp_SeqSELL(Mat A)
1625: {
1626: PetscFunctionBegin;
1627: PetscCall(MatSeqSELLSetPreallocation(A, PETSC_DEFAULT, NULL));
1628: PetscFunctionReturn(PETSC_SUCCESS);
1629: }
1631: PetscErrorCode MatSeqSELLGetArray_SeqSELL(Mat A, PetscScalar *array[])
1632: {
1633: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1635: PetscFunctionBegin;
1636: *array = a->val;
1637: PetscFunctionReturn(PETSC_SUCCESS);
1638: }
1640: PetscErrorCode MatSeqSELLRestoreArray_SeqSELL(Mat A, PetscScalar *array[])
1641: {
1642: PetscFunctionBegin;
1643: PetscFunctionReturn(PETSC_SUCCESS);
1644: }
1646: PetscErrorCode MatScale_SeqSELL(Mat inA, PetscScalar alpha)
1647: {
1648: Mat_SeqSELL *a = (Mat_SeqSELL *)inA->data;
1649: MatScalar *aval = a->val;
1650: PetscScalar oalpha = alpha;
1651: PetscBLASInt one = 1, size;
1653: PetscFunctionBegin;
1654: PetscCall(PetscBLASIntCast(a->sliidx[a->totalslices], &size));
1655: PetscCallBLAS("BLASscal", BLASscal_(&size, &oalpha, aval, &one));
1656: PetscCall(PetscLogFlops(a->nz));
1657: #if PetscDefined(HAVE_CUPM)
1658: if (inA->offloadmask != PETSC_OFFLOAD_UNALLOCATED) inA->offloadmask = PETSC_OFFLOAD_CPU;
1659: #endif
1660: PetscFunctionReturn(PETSC_SUCCESS);
1661: }
1663: PetscErrorCode MatShift_SeqSELL(Mat Y, PetscScalar a)
1664: {
1665: Mat_SeqSELL *y = (Mat_SeqSELL *)Y->data;
1667: PetscFunctionBegin;
1668: if (!Y->preallocated || !y->nz) PetscCall(MatSeqSELLSetPreallocation(Y, 1, NULL));
1669: PetscCall(MatShift_Basic(Y, a));
1670: PetscFunctionReturn(PETSC_SUCCESS);
1671: }
1673: PetscErrorCode MatSOR_SeqSELL(Mat A, Vec bb, PetscReal omega, MatSORType flag, PetscReal fshift, PetscInt its, PetscInt lits, Vec xx)
1674: {
1675: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
1676: PetscScalar *x, sum, *t;
1677: const MatScalar *idiag = NULL, *mdiag;
1678: const PetscScalar *b, *xb;
1679: PetscInt n, m = A->rmap->n, i, j, shift;
1680: const PetscInt *diag;
1682: PetscFunctionBegin;
1683: its = its * lits;
1685: PetscCall(MatInvertDiagonalForSOR_SeqSELL(A, omega, fshift));
1686: diag = a->diag;
1687: t = a->ssor_work;
1688: idiag = a->idiag;
1689: mdiag = a->mdiag;
1691: PetscCall(VecGetArray(xx, &x));
1692: PetscCall(VecGetArrayRead(bb, &b));
1693: /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */
1694: PetscCheck(flag != SOR_APPLY_UPPER, PETSC_COMM_SELF, PETSC_ERR_SUP, "SOR_APPLY_UPPER is not implemented");
1695: PetscCheck(flag != SOR_APPLY_LOWER, PETSC_COMM_SELF, PETSC_ERR_SUP, "SOR_APPLY_LOWER is not implemented");
1696: PetscCheck(!(flag & SOR_EISENSTAT), PETSC_COMM_SELF, PETSC_ERR_SUP, "No support yet for Eisenstat");
1698: if (flag & SOR_ZERO_INITIAL_GUESS) {
1699: if ((flag & SOR_FORWARD_SWEEP) || (flag & SOR_LOCAL_FORWARD_SWEEP)) {
1700: for (i = 0; i < m; i++) {
1701: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight; /* starting index of the row i */
1702: sum = b[i];
1703: n = (diag[i] - shift) / a->sliceheight;
1704: for (j = 0; j < n; j++) sum -= a->val[shift + a->sliceheight * j] * x[a->colidx[shift + a->sliceheight * j]];
1705: t[i] = sum;
1706: x[i] = sum * idiag[i];
1707: }
1708: xb = t;
1709: PetscCall(PetscLogFlops(a->nz));
1710: } else xb = b;
1711: if ((flag & SOR_BACKWARD_SWEEP) || (flag & SOR_LOCAL_BACKWARD_SWEEP)) {
1712: for (i = m - 1; i >= 0; i--) {
1713: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight; /* starting index of the row i */
1714: sum = xb[i];
1715: n = a->rlen[i] - (diag[i] - shift) / a->sliceheight - 1;
1716: for (j = 1; j <= n; j++) sum -= a->val[diag[i] + a->sliceheight * j] * x[a->colidx[diag[i] + a->sliceheight * j]];
1717: if (xb == b) {
1718: x[i] = sum * idiag[i];
1719: } else {
1720: x[i] = (1. - omega) * x[i] + sum * idiag[i]; /* omega in idiag */
1721: }
1722: }
1723: PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */
1724: }
1725: its--;
1726: }
1727: while (its--) {
1728: if ((flag & SOR_FORWARD_SWEEP) || (flag & SOR_LOCAL_FORWARD_SWEEP)) {
1729: for (i = 0; i < m; i++) {
1730: /* lower */
1731: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight; /* starting index of the row i */
1732: sum = b[i];
1733: n = (diag[i] - shift) / a->sliceheight;
1734: for (j = 0; j < n; j++) sum -= a->val[shift + a->sliceheight * j] * x[a->colidx[shift + a->sliceheight * j]];
1735: t[i] = sum; /* save application of the lower-triangular part */
1736: /* upper */
1737: n = a->rlen[i] - (diag[i] - shift) / a->sliceheight - 1;
1738: for (j = 1; j <= n; j++) sum -= a->val[diag[i] + a->sliceheight * j] * x[a->colidx[diag[i] + a->sliceheight * j]];
1739: x[i] = (1. - omega) * x[i] + sum * idiag[i]; /* omega in idiag */
1740: }
1741: xb = t;
1742: PetscCall(PetscLogFlops(2.0 * a->nz));
1743: } else xb = b;
1744: if ((flag & SOR_BACKWARD_SWEEP) || (flag & SOR_LOCAL_BACKWARD_SWEEP)) {
1745: for (i = m - 1; i >= 0; i--) {
1746: shift = a->sliidx[i / a->sliceheight] + i % a->sliceheight; /* starting index of the row i */
1747: sum = xb[i];
1748: if (xb == b) {
1749: /* whole matrix (no checkpointing available) */
1750: n = a->rlen[i];
1751: for (j = 0; j < n; j++) sum -= a->val[shift + a->sliceheight * j] * x[a->colidx[shift + a->sliceheight * j]];
1752: x[i] = (1. - omega) * x[i] + (sum + mdiag[i] * x[i]) * idiag[i];
1753: } else { /* lower-triangular part has been saved, so only apply upper-triangular */
1754: n = a->rlen[i] - (diag[i] - shift) / a->sliceheight - 1;
1755: for (j = 1; j <= n; j++) sum -= a->val[diag[i] + a->sliceheight * j] * x[a->colidx[diag[i] + a->sliceheight * j]];
1756: x[i] = (1. - omega) * x[i] + sum * idiag[i]; /* omega in idiag */
1757: }
1758: }
1759: if (xb == b) PetscCall(PetscLogFlops(2.0 * a->nz));
1760: else PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */
1761: }
1762: }
1763: PetscCall(VecRestoreArray(xx, &x));
1764: PetscCall(VecRestoreArrayRead(bb, &b));
1765: PetscFunctionReturn(PETSC_SUCCESS);
1766: }
1768: static struct _MatOps MatOps_Values = {MatSetValues_SeqSELL,
1769: MatGetRow_SeqSELL,
1770: MatRestoreRow_SeqSELL,
1771: MatMult_SeqSELL,
1772: /* 4*/ MatMultAdd_SeqSELL,
1773: MatMultTranspose_SeqSELL,
1774: MatMultTransposeAdd_SeqSELL,
1775: NULL,
1776: NULL,
1777: NULL,
1778: /* 10*/ NULL,
1779: NULL,
1780: NULL,
1781: MatSOR_SeqSELL,
1782: NULL,
1783: /* 15*/ MatGetInfo_SeqSELL,
1784: MatEqual_SeqSELL,
1785: MatGetDiagonal_SeqSELL,
1786: MatDiagonalScale_SeqSELL,
1787: NULL,
1788: /* 20*/ NULL,
1789: MatAssemblyEnd_SeqSELL,
1790: MatSetOption_SeqSELL,
1791: MatZeroEntries_SeqSELL,
1792: /* 24*/ NULL,
1793: NULL,
1794: NULL,
1795: NULL,
1796: NULL,
1797: /* 29*/ MatSetUp_SeqSELL,
1798: NULL,
1799: NULL,
1800: NULL,
1801: NULL,
1802: /* 34*/ MatDuplicate_SeqSELL,
1803: NULL,
1804: NULL,
1805: NULL,
1806: NULL,
1807: /* 39*/ NULL,
1808: NULL,
1809: NULL,
1810: MatGetValues_SeqSELL,
1811: MatCopy_SeqSELL,
1812: /* 44*/ NULL,
1813: MatScale_SeqSELL,
1814: MatShift_SeqSELL,
1815: NULL,
1816: NULL,
1817: /* 49*/ NULL,
1818: NULL,
1819: NULL,
1820: NULL,
1821: NULL,
1822: /* 54*/ MatFDColoringCreate_SeqXAIJ,
1823: NULL,
1824: NULL,
1825: NULL,
1826: NULL,
1827: /* 59*/ NULL,
1828: MatDestroy_SeqSELL,
1829: MatView_SeqSELL,
1830: NULL,
1831: NULL,
1832: /* 64*/ NULL,
1833: NULL,
1834: NULL,
1835: NULL,
1836: NULL,
1837: /* 69*/ NULL,
1838: NULL,
1839: NULL,
1840: MatFDColoringApply_AIJ, /* reuse the FDColoring function for AIJ */
1841: NULL,
1842: /* 74*/ NULL,
1843: NULL,
1844: NULL,
1845: NULL,
1846: NULL,
1847: /* 79*/ NULL,
1848: NULL,
1849: NULL,
1850: NULL,
1851: NULL,
1852: /* 84*/ NULL,
1853: NULL,
1854: NULL,
1855: NULL,
1856: NULL,
1857: /* 89*/ NULL,
1858: NULL,
1859: NULL,
1860: NULL,
1861: MatConjugate_SeqSELL,
1862: /* 94*/ NULL,
1863: NULL,
1864: NULL,
1865: NULL,
1866: NULL,
1867: /* 99*/ NULL,
1868: NULL,
1869: NULL,
1870: NULL,
1871: NULL,
1872: /*104*/ NULL,
1873: NULL,
1874: NULL,
1875: NULL,
1876: NULL,
1877: /*109*/ NULL,
1878: NULL,
1879: NULL,
1880: NULL,
1881: NULL,
1882: /*114*/ NULL,
1883: NULL,
1884: NULL,
1885: NULL,
1886: NULL,
1887: /*119*/ NULL,
1888: NULL,
1889: NULL,
1890: NULL,
1891: NULL,
1892: /*124*/ NULL,
1893: NULL,
1894: NULL,
1895: MatFDColoringSetUp_SeqXAIJ,
1896: NULL,
1897: /*129*/ NULL,
1898: NULL,
1899: NULL,
1900: NULL,
1901: NULL,
1902: /*134*/ NULL,
1903: NULL,
1904: NULL,
1905: NULL,
1906: NULL,
1907: /*139*/ NULL,
1908: NULL,
1909: NULL,
1910: NULL,
1911: NULL,
1912: /*144*/ NULL,
1913: NULL,
1914: NULL,
1915: NULL};
1917: static PetscErrorCode MatStoreValues_SeqSELL(Mat mat)
1918: {
1919: Mat_SeqSELL *a = (Mat_SeqSELL *)mat->data;
1921: PetscFunctionBegin;
1922: PetscCheck(a->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1924: /* allocate space for values if not already there */
1925: if (!a->saved_values) PetscCall(PetscMalloc1(a->sliidx[a->totalslices] + 1, &a->saved_values));
1927: /* copy values over */
1928: PetscCall(PetscArraycpy(a->saved_values, a->val, a->sliidx[a->totalslices]));
1929: PetscFunctionReturn(PETSC_SUCCESS);
1930: }
1932: static PetscErrorCode MatRetrieveValues_SeqSELL(Mat mat)
1933: {
1934: Mat_SeqSELL *a = (Mat_SeqSELL *)mat->data;
1936: PetscFunctionBegin;
1937: PetscCheck(a->nonew, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first");
1938: PetscCheck(a->saved_values, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Must call MatStoreValues(A);first");
1939: PetscCall(PetscArraycpy(a->val, a->saved_values, a->sliidx[a->totalslices]));
1940: #if PetscDefined(HAVE_CUPM)
1941: if (mat->offloadmask != PETSC_OFFLOAD_UNALLOCATED) mat->offloadmask = PETSC_OFFLOAD_CPU;
1942: #endif
1943: PetscFunctionReturn(PETSC_SUCCESS);
1944: }
1946: static PetscErrorCode MatSeqSELLGetFillRatio_SeqSELL(Mat mat, PetscReal *ratio)
1947: {
1948: Mat_SeqSELL *a = (Mat_SeqSELL *)mat->data;
1950: PetscFunctionBegin;
1951: if (a->totalslices && a->sliidx[a->totalslices]) {
1952: *ratio = (PetscReal)(a->sliidx[a->totalslices] - a->nz) / a->sliidx[a->totalslices];
1953: } else {
1954: *ratio = 0.0;
1955: }
1956: PetscFunctionReturn(PETSC_SUCCESS);
1957: }
1959: static PetscErrorCode MatSeqSELLGetMaxSliceWidth_SeqSELL(Mat mat, PetscInt *slicewidth)
1960: {
1961: Mat_SeqSELL *a = (Mat_SeqSELL *)mat->data;
1962: PetscInt i, current_slicewidth;
1964: PetscFunctionBegin;
1965: *slicewidth = 0;
1966: for (i = 0; i < a->totalslices; i++) {
1967: current_slicewidth = (a->sliidx[i + 1] - a->sliidx[i]) / a->sliceheight;
1968: if (current_slicewidth > *slicewidth) *slicewidth = current_slicewidth;
1969: }
1970: PetscFunctionReturn(PETSC_SUCCESS);
1971: }
1973: static PetscErrorCode MatSeqSELLGetAvgSliceWidth_SeqSELL(Mat mat, PetscReal *slicewidth)
1974: {
1975: Mat_SeqSELL *a = (Mat_SeqSELL *)mat->data;
1977: PetscFunctionBegin;
1978: *slicewidth = 0;
1979: if (a->totalslices) *slicewidth = (PetscReal)a->sliidx[a->totalslices] / a->sliceheight / a->totalslices;
1980: PetscFunctionReturn(PETSC_SUCCESS);
1981: }
1983: static PetscErrorCode MatSeqSELLGetVarSliceSize_SeqSELL(Mat mat, PetscReal *variance)
1984: {
1985: Mat_SeqSELL *a = (Mat_SeqSELL *)mat->data;
1986: PetscReal mean;
1987: PetscInt i, totalslices = a->totalslices, *sliidx = a->sliidx;
1989: PetscFunctionBegin;
1990: *variance = 0;
1991: if (totalslices) {
1992: mean = (PetscReal)sliidx[totalslices] / totalslices;
1993: for (i = 1; i <= totalslices; i++) *variance += ((PetscReal)(sliidx[i] - sliidx[i - 1]) - mean) * ((PetscReal)(sliidx[i] - sliidx[i - 1]) - mean) / totalslices;
1994: }
1995: PetscFunctionReturn(PETSC_SUCCESS);
1996: }
1998: static PetscErrorCode MatSeqSELLSetSliceHeight_SeqSELL(Mat A, PetscInt sliceheight)
1999: {
2000: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
2002: PetscFunctionBegin;
2003: if (A->preallocated) PetscFunctionReturn(PETSC_SUCCESS);
2004: PetscCheck(a->sliceheight <= 0 || a->sliceheight == sliceheight, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot change slice height %" PetscInt_FMT " to %" PetscInt_FMT, a->sliceheight, sliceheight);
2005: a->sliceheight = sliceheight;
2006: #if PetscDefined(HAVE_CUPM)
2007: PetscCheck(PetscMax(DEVICE_MEM_ALIGN, sliceheight) % PetscMin(DEVICE_MEM_ALIGN, sliceheight) == 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "The slice height is not compatible with DEVICE_MEM_ALIGN (one must be divisible by the other) %" PetscInt_FMT, sliceheight);
2008: #endif
2009: PetscFunctionReturn(PETSC_SUCCESS);
2010: }
2012: /*@
2013: MatSeqSELLGetFillRatio - returns a ratio that indicates the irregularity of the matrix.
2015: Not Collective
2017: Input Parameter:
2018: . A - a MATSEQSELL matrix
2020: Output Parameter:
2021: . ratio - ratio of number of padded zeros to number of allocated elements
2023: Level: intermediate
2025: .seealso: `MATSEQSELL`, `MatSeqSELLGetAvgSliceWidth()`
2026: @*/
2027: PetscErrorCode MatSeqSELLGetFillRatio(Mat A, PetscReal *ratio)
2028: {
2029: PetscFunctionBegin;
2030: PetscUseMethod(A, "MatSeqSELLGetFillRatio_C", (Mat, PetscReal *), (A, ratio));
2031: PetscFunctionReturn(PETSC_SUCCESS);
2032: }
2034: /*@
2035: MatSeqSELLGetMaxSliceWidth - returns the maximum slice width.
2037: Not Collective
2039: Input Parameter:
2040: . A - a MATSEQSELL matrix
2042: Output Parameter:
2043: . slicewidth - maximum slice width
2045: Level: intermediate
2047: .seealso: `MATSEQSELL`, `MatSeqSELLGetAvgSliceWidth()`
2048: @*/
2049: PetscErrorCode MatSeqSELLGetMaxSliceWidth(Mat A, PetscInt *slicewidth)
2050: {
2051: PetscFunctionBegin;
2052: PetscUseMethod(A, "MatSeqSELLGetMaxSliceWidth_C", (Mat, PetscInt *), (A, slicewidth));
2053: PetscFunctionReturn(PETSC_SUCCESS);
2054: }
2056: /*@
2057: MatSeqSELLGetAvgSliceWidth - returns the average slice width.
2059: Not Collective
2061: Input Parameter:
2062: . A - a MATSEQSELL matrix
2064: Output Parameter:
2065: . slicewidth - average slice width
2067: Level: intermediate
2069: .seealso: `MATSEQSELL`, `MatSeqSELLGetMaxSliceWidth()`
2070: @*/
2071: PetscErrorCode MatSeqSELLGetAvgSliceWidth(Mat A, PetscReal *slicewidth)
2072: {
2073: PetscFunctionBegin;
2074: PetscUseMethod(A, "MatSeqSELLGetAvgSliceWidth_C", (Mat, PetscReal *), (A, slicewidth));
2075: PetscFunctionReturn(PETSC_SUCCESS);
2076: }
2078: /*@
2079: MatSeqSELLSetSliceHeight - sets the slice height.
2081: Not Collective
2083: Input Parameters:
2084: + A - a MATSEQSELL matrix
2085: - sliceheight - slice height
2087: Notes:
2088: You cannot change the slice height once it have been set.
2090: The slice height must be set before MatSetUp() or MatXXXSetPreallocation() is called.
2092: Level: intermediate
2094: .seealso: `MATSEQSELL`, `MatSeqSELLGetVarSliceSize()`
2095: @*/
2096: PetscErrorCode MatSeqSELLSetSliceHeight(Mat A, PetscInt sliceheight)
2097: {
2098: PetscFunctionBegin;
2099: PetscUseMethod(A, "MatSeqSELLSetSliceHeight_C", (Mat, PetscInt), (A, sliceheight));
2100: PetscFunctionReturn(PETSC_SUCCESS);
2101: }
2103: /*@
2104: MatSeqSELLGetVarSliceSize - returns the variance of the slice size.
2106: Not Collective
2108: Input Parameter:
2109: . A - a MATSEQSELL matrix
2111: Output Parameter:
2112: . variance - variance of the slice size
2114: Level: intermediate
2116: .seealso: `MATSEQSELL`, `MatSeqSELLSetSliceHeight()`
2117: @*/
2118: PetscErrorCode MatSeqSELLGetVarSliceSize(Mat A, PetscReal *variance)
2119: {
2120: PetscFunctionBegin;
2121: PetscUseMethod(A, "MatSeqSELLGetVarSliceSize_C", (Mat, PetscReal *), (A, variance));
2122: PetscFunctionReturn(PETSC_SUCCESS);
2123: }
2125: #if PetscDefined(HAVE_CUDA)
2126: PETSC_EXTERN PetscErrorCode MatConvert_SeqSELL_SeqSELLCUDA(Mat);
2127: #endif
2128: #if PetscDefined(HAVE_HIP)
2129: PETSC_EXTERN PetscErrorCode MatConvert_SeqSELL_SeqSELLHIP(Mat);
2130: #endif
2132: PETSC_EXTERN PetscErrorCode MatCreate_SeqSELL(Mat B)
2133: {
2134: Mat_SeqSELL *b;
2135: PetscMPIInt size;
2137: PetscFunctionBegin;
2138: PetscCall(PetscCitationsRegister(citation, &cited));
2139: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B), &size));
2140: PetscCheck(size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Comm must be of size 1");
2142: PetscCall(PetscNew(&b));
2144: B->data = (void *)b;
2145: B->ops[0] = MatOps_Values;
2147: b->row = NULL;
2148: b->col = NULL;
2149: b->icol = NULL;
2150: b->reallocs = 0;
2151: b->ignorezeroentries = PETSC_FALSE;
2152: b->roworiented = PETSC_TRUE;
2153: b->nonew = 0;
2154: b->diag = NULL;
2155: b->solve_work = NULL;
2156: B->spptr = NULL;
2157: b->saved_values = NULL;
2158: b->idiag = NULL;
2159: b->mdiag = NULL;
2160: b->ssor_work = NULL;
2161: b->omega = 1.0;
2162: b->fshift = 0.0;
2163: b->keepnonzeropattern = PETSC_FALSE;
2164: b->sliceheight = 0;
2166: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATSEQSELL));
2167: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSELLGetArray_C", MatSeqSELLGetArray_SeqSELL));
2168: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSELLRestoreArray_C", MatSeqSELLRestoreArray_SeqSELL));
2169: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatStoreValues_C", MatStoreValues_SeqSELL));
2170: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatRetrieveValues_C", MatRetrieveValues_SeqSELL));
2171: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSELLSetPreallocation_C", MatSeqSELLSetPreallocation_SeqSELL));
2172: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsell_seqaij_C", MatConvert_SeqSELL_SeqAIJ));
2173: #if PetscDefined(HAVE_CUDA)
2174: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsell_seqsellcuda_C", MatConvert_SeqSELL_SeqSELLCUDA));
2175: #endif
2176: #if PetscDefined(HAVE_HIP)
2177: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_seqsell_seqsellhip_C", MatConvert_SeqSELL_SeqSELLHIP));
2178: #endif
2179: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSELLGetFillRatio_C", MatSeqSELLGetFillRatio_SeqSELL));
2180: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSELLGetMaxSliceWidth_C", MatSeqSELLGetMaxSliceWidth_SeqSELL));
2181: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSELLGetAvgSliceWidth_C", MatSeqSELLGetAvgSliceWidth_SeqSELL));
2182: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSELLGetVarSliceSize_C", MatSeqSELLGetVarSliceSize_SeqSELL));
2183: PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSeqSELLSetSliceHeight_C", MatSeqSELLSetSliceHeight_SeqSELL));
2185: PetscObjectOptionsBegin((PetscObject)B);
2186: {
2187: PetscInt newsh = -1;
2188: PetscBool flg;
2189: #if PetscDefined(HAVE_CUPM)
2190: PetscInt chunksize = 0;
2191: #endif
2193: PetscCall(PetscOptionsInt("-mat_sell_slice_height", "Set the slice height used to store SELL matrix", "MatSELLSetSliceHeight", newsh, &newsh, &flg));
2194: if (flg) PetscCall(MatSeqSELLSetSliceHeight(B, newsh));
2195: #if PetscDefined(HAVE_CUPM)
2196: PetscCall(PetscOptionsInt("-mat_sell_chunk_size", "Set the chunksize for load-balanced CUDA/HIP kernels. Choices include 64,128,256,512,1024", NULL, chunksize, &chunksize, &flg));
2197: if (flg) {
2198: PetscCheck(chunksize >= 64 && chunksize <= 1024 && chunksize % 64 == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "chunksize must be a number in {64,128,256,512,1024}: value %" PetscInt_FMT, chunksize);
2199: b->chunksize = chunksize;
2200: }
2201: #endif
2202: }
2203: PetscOptionsEnd();
2204: PetscFunctionReturn(PETSC_SUCCESS);
2205: }
2207: /*
2208: Given a matrix generated with MatGetFactor() duplicates all the information in A into B
2209: */
2210: static PetscErrorCode MatDuplicateNoCreate_SeqSELL(Mat C, Mat A, MatDuplicateOption cpvalues, PetscBool mallocmatspace)
2211: {
2212: Mat_SeqSELL *c = (Mat_SeqSELL *)C->data, *a = (Mat_SeqSELL *)A->data;
2213: PetscInt i, m = A->rmap->n;
2214: PetscInt totalslices = a->totalslices;
2216: PetscFunctionBegin;
2217: C->factortype = A->factortype;
2218: c->row = NULL;
2219: c->col = NULL;
2220: c->icol = NULL;
2221: c->reallocs = 0;
2222: C->assembled = PETSC_TRUE;
2224: PetscCall(PetscLayoutReference(A->rmap, &C->rmap));
2225: PetscCall(PetscLayoutReference(A->cmap, &C->cmap));
2227: c->sliceheight = a->sliceheight;
2228: c->totalslices = a->totalslices;
2229: #if PetscDefined(HAVE_CUPM)
2230: c->chunksize = a->chunksize;
2231: c->totalchunks = a->totalchunks;
2232: PetscCall(PetscMalloc1(c->totalchunks, &c->chunk_slice_map));
2233: PetscCall(PetscArraycpy(c->chunk_slice_map, a->chunk_slice_map, c->totalchunks));
2234: #endif
2235: PetscCall(PetscCalloc1(c->sliceheight * totalslices, &c->rlen));
2236: PetscCall(PetscMalloc1(totalslices + 1, &c->sliidx));
2238: for (i = 0; i < m; i++) c->rlen[i] = a->rlen[i];
2239: for (i = 0; i < totalslices + 1; i++) c->sliidx[i] = a->sliidx[i];
2241: /* allocate the matrix space */
2242: if (mallocmatspace) {
2243: PetscCall(PetscMalloc2(a->maxallocmat, &c->val, a->maxallocmat, &c->colidx));
2245: c->singlemalloc = PETSC_TRUE;
2247: if (m > 0) {
2248: PetscCall(PetscArraycpy(c->colidx, a->colidx, a->maxallocmat));
2249: if (cpvalues == MAT_COPY_VALUES) {
2250: PetscCall(PetscArraycpy(c->val, a->val, a->maxallocmat));
2251: } else {
2252: PetscCall(PetscArrayzero(c->val, a->maxallocmat));
2253: }
2254: }
2255: }
2257: c->ignorezeroentries = a->ignorezeroentries;
2258: c->roworiented = a->roworiented;
2259: c->nonew = a->nonew;
2260: c->solve_work = NULL;
2261: c->saved_values = NULL;
2262: c->idiag = NULL;
2263: c->ssor_work = NULL;
2264: c->keepnonzeropattern = a->keepnonzeropattern;
2265: c->free_val = PETSC_TRUE;
2266: c->free_colidx = PETSC_TRUE;
2268: c->maxallocmat = a->maxallocmat;
2269: c->maxallocrow = a->maxallocrow;
2270: c->rlenmax = a->rlenmax;
2271: c->nz = a->nz;
2272: C->preallocated = PETSC_TRUE;
2274: c->nonzerorowcnt = a->nonzerorowcnt;
2275: C->nonzerostate = A->nonzerostate;
2277: PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist, &((PetscObject)C)->qlist));
2278: PetscFunctionReturn(PETSC_SUCCESS);
2279: }
2281: PetscErrorCode MatDuplicate_SeqSELL(Mat A, MatDuplicateOption cpvalues, Mat *B)
2282: {
2283: PetscFunctionBegin;
2284: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), B));
2285: PetscCall(MatSetSizes(*B, A->rmap->n, A->cmap->n, A->rmap->n, A->cmap->n));
2286: if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) PetscCall(MatSetBlockSizesFromMats(*B, A, A));
2287: PetscCall(MatSetType(*B, ((PetscObject)A)->type_name));
2288: PetscCall(MatDuplicateNoCreate_SeqSELL(*B, A, cpvalues, PETSC_TRUE));
2289: PetscFunctionReturn(PETSC_SUCCESS);
2290: }
2292: /*MC
2293: MATSEQSELL - MATSEQSELL = "seqsell" - A matrix type to be used for sequential sparse matrices,
2294: based on the sliced Ellpack format, {cite}`zhangellpack2018`
2296: Options Database Key:
2297: . -mat_type seqsell - sets the matrix type to "`MATSEQELL` during a call to `MatSetFromOptions()`
2299: Level: beginner
2301: .seealso: `Mat`, `MatCreateSeqSELL()`, `MATSELL`, `MATMPISELL`, `MATSEQAIJ`, `MATAIJ`, `MATMPIAIJ`
2302: M*/
2304: /*MC
2305: MATSELL - MATSELL = "sell" - A matrix type to be used for sparse matrices, {cite}`zhangellpack2018`
2307: This matrix type is identical to `MATSEQSELL` when constructed with a single process communicator,
2308: and `MATMPISELL` otherwise. As a result, for single process communicators,
2309: `MatSeqSELLSetPreallocation()` is supported, and similarly `MatMPISELLSetPreallocation()` is supported
2310: for communicators controlling multiple processes. It is recommended that you call both of
2311: the above preallocation routines for simplicity.
2313: Options Database Key:
2314: . -mat_type sell - sets the matrix type to "sell" during a call to MatSetFromOptions()
2316: Level: beginner
2318: Notes:
2319: This format is only supported for real scalars, double precision, and 32-bit indices (the defaults).
2321: It can provide better performance on Intel and AMD processes with AVX2 or AVX512 support for matrices that have a similar number of
2322: non-zeros in contiguous groups of rows. However if the computation is memory bandwidth limited it may not provide much improvement.
2324: Developer Notes:
2325: On Intel (and AMD) systems some of the matrix operations use SIMD (AVX) instructions to achieve higher performance.
2327: The sparse matrix format is as follows. For simplicity we assume a slice size of 2, it is actually 8
2328: .vb
2329: (2 0 3 4)
2330: Consider the matrix A = (5 0 6 0)
2331: (0 0 7 8)
2332: (0 0 9 9)
2334: symbolically the Ellpack format can be written as
2336: (2 3 4 |) (0 2 3 |)
2337: v = (5 6 0 |) colidx = (0 2 2 |)
2338: -------- ---------
2339: (7 8 |) (2 3 |)
2340: (9 9 |) (2 3 |)
2342: The data for 2 contiguous rows of the matrix are stored together (in column-major format) (with any left-over rows handled as a special case).
2343: Any of the rows in a slice fewer columns than the rest of the slice (row 1 above) are padded with a previous valid column in their "extra" colidx[] locations and
2344: zeros in their "extra" v locations so that the matrix operations do not need special code to handle different length rows within the 2 rows in a slice.
2346: The one-dimensional representation of v used in the code is (2 5 3 6 4 0 7 9 8 9) and for colidx is (0 0 2 2 3 2 2 2 3 3)
2348: .ve
2350: See `MatMult_SeqSELL()` for how this format is used with the SIMD operations to achieve high performance.
2352: .seealso: `Mat`, `MatCreateSeqSELL()`, `MatCreateSeqAIJ()`, `MatCreateSELL()`, `MATSEQSELL`, `MATMPISELL`, `MATSEQAIJ`, `MATMPIAIJ`, `MATAIJ`
2353: M*/
2355: /*@
2356: MatCreateSeqSELL - Creates a sparse matrix in `MATSEQSELL` format.
2358: Collective
2360: Input Parameters:
2361: + comm - MPI communicator, set to `PETSC_COMM_SELF`
2362: . m - number of rows
2363: . n - number of columns
2364: . rlenmax - maximum number of nonzeros in a row, ignored if `rlen` is provided
2365: - rlen - array containing the number of nonzeros in the various rows (possibly different for each row) or NULL
2367: Output Parameter:
2368: . A - the matrix
2370: Level: intermediate
2372: Notes:
2373: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`,
2374: MatXXXXSetPreallocation() paradigm instead of this routine directly.
2375: [MatXXXXSetPreallocation() is, for example, `MatSeqSELLSetPreallocation()`]
2377: Specify the preallocated storage with either `rlenmax` or `rlen` (not both).
2378: Set `rlenmax` = `PETSC_DEFAULT` and `rlen` = `NULL` for PETSc to control dynamic memory
2379: allocation.
2381: .seealso: `Mat`, `MATSEQSELL`, `MatCreate()`, `MatCreateSELL()`, `MatSetValues()`, `MatSeqSELLSetPreallocation()`, `MATSELL`, `MATMPISELL`
2382: @*/
2383: PetscErrorCode MatCreateSeqSELL(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt rlenmax, const PetscInt rlen[], Mat *A)
2384: {
2385: PetscFunctionBegin;
2386: PetscCall(MatCreate(comm, A));
2387: PetscCall(MatSetSizes(*A, m, n, m, n));
2388: PetscCall(MatSetType(*A, MATSEQSELL));
2389: PetscCall(MatSeqSELLSetPreallocation_SeqSELL(*A, rlenmax, rlen));
2390: PetscFunctionReturn(PETSC_SUCCESS);
2391: }
2393: PetscErrorCode MatEqual_SeqSELL(Mat A, Mat B, PetscBool *flg)
2394: {
2395: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data, *b = (Mat_SeqSELL *)B->data;
2396: PetscInt totalslices = a->totalslices;
2398: PetscFunctionBegin;
2399: /* If the matrix dimensions are not equal,or no of nonzeros */
2400: if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) || (a->nz != b->nz) || (a->rlenmax != b->rlenmax)) {
2401: *flg = PETSC_FALSE;
2402: PetscFunctionReturn(PETSC_SUCCESS);
2403: }
2404: /* if the a->colidx are the same */
2405: PetscCall(PetscArraycmp(a->colidx, b->colidx, a->sliidx[totalslices], flg));
2406: if (!*flg) PetscFunctionReturn(PETSC_SUCCESS);
2407: /* if a->val are the same */
2408: PetscCall(PetscArraycmp(a->val, b->val, a->sliidx[totalslices], flg));
2409: PetscFunctionReturn(PETSC_SUCCESS);
2410: }
2412: PetscErrorCode MatConjugate_SeqSELL(Mat A)
2413: {
2414: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data;
2415: PetscScalar *val = a->val;
2417: PetscFunctionBegin;
2418: for (PetscInt i = 0; i < a->sliidx[a->totalslices]; i++) val[i] = PetscConj(val[i]);
2419: #if PetscDefined(HAVE_CUPM)
2420: if (A->offloadmask != PETSC_OFFLOAD_UNALLOCATED) A->offloadmask = PETSC_OFFLOAD_CPU;
2421: #endif
2422: PetscFunctionReturn(PETSC_SUCCESS);
2423: }