Actual source code: sell.h
1: #pragma once
3: #include <petsc/private/matimpl.h>
4: #include <petsc/private/hashmapi.h>
6: /*
7: For NVIDIA GPUs each slice should be padded to the boundary of 16 elements for best performance.
8: The optimal memory alignment in device memory is 128 bytes, 64 bytes, 32 bytes for double precision, single precision and half precision.
9: */
10: #if defined(PETSC_HAVE_DEVICE)
11: #define DEVICE_MEM_ALIGN 16
12: #endif
14: /*
15: Struct header for SeqSELL matrix format
16: */
17: #define SEQSELLHEADER(datatype) \
18: PetscBool roworiented; /* if true, row-oriented input, default */ \
19: PetscInt nonew; /* 1 don't add new nonzeros, -1 generate error on new */ \
20: PetscInt nounused; /* -1 generate error on unused space */ \
21: PetscBool singlemalloc; /* if true a, i, and j have been obtained with one big malloc */ \
22: PetscInt maxallocmat; /* max allocated space for the matrix */ \
23: PetscInt maxallocrow; /* max allocated space for each row */ \
24: PetscInt nz; /* actual nonzeros */ \
25: PetscInt rlenmax; /* max actual row length, rmax cannot exceed maxallocrow */ \
26: PetscInt *rlen; /* actual length of each row (padding zeros excluded) */ \
27: PetscBool free_rlen; /* free rlen array ? */ \
28: PetscInt reallocs; /* number of mallocs done during MatSetValues() \
29: as more values are set than were prealloced */ \
30: PetscBool keepnonzeropattern; /* keeps matrix nonzero structure the same in calls to MatZeroRows()*/ \
31: PetscBool ignorezeroentries; \
32: PetscBool free_colidx; /* free the column indices colidx when the matrix is destroyed */ \
33: PetscBool free_val; /* free the numerical values when matrix is destroy */ \
34: PetscInt *colidx; /* column index */ \
35: PetscInt *diag; /* pointers to diagonal elements */ \
36: PetscObjectState diagNonzeroState; /* nonzero state of the matrix when diag was obtained */ \
37: PetscBool diagDense; /* matrix contains all the diagonal entries */ \
38: PetscInt nonzerorowcnt; /* how many rows have nonzero entries */ \
39: datatype *val; /* elements including nonzeros and padding zeros */ \
40: PetscScalar *solve_work; /* work space used in MatSolve */ \
41: IS row, col, icol; /* index sets, used for reorderings */ \
42: PetscBool pivotinblocks; /* pivot inside factorization of each diagonal block */ \
43: Mat parent; /* set if this matrix was formed with MatDuplicate(...,MAT_SHARE_NONZERO_PATTERN,....);
44: means that this shares some data structures with the parent including diag, ilen, imax, i, j */ \
45: PetscInt *sliidx; /* slice index */ \
46: PetscInt totalslices; /* total number of slices */ \
47: PetscInt sliceheight; /* slice height */ \
48: PetscReal fillratio; /* ratio of number of padded zeros over total number of elements */ \
49: PetscReal avgslicewidth; /* average slice width */ \
50: PetscInt maxslicewidth; /* maximum slice width */ \
51: PetscReal varslicesize; /* variance of slice size */ \
52: PetscInt *sliperm; /* slice permutation array, CUDA only */ \
53: PetscInt totalblocks; /* total number of blocks, CUDA only */ \
54: PetscInt *blockidx; /* block index, CUDA only */ \
55: PetscInt *block_row_map; /* starting row of the current block, CUDA only */ \
56: PetscInt chunksize; /* chunk size, CUDA only */ \
57: PetscInt totalchunks; /* total number of chunks, CUDA only */ \
58: PetscInt *chunk_slice_map; /* starting slice of the current chunk, CUDA only */ \
59: PetscInt *getrowcols; /* workarray for MatGetRow_SeqSELL */ \
60: PetscScalar *getrowvals /* workarray for MatGetRow_SeqSELL */
62: typedef struct {
63: SEQSELLHEADER(MatScalar);
64: MatScalar *saved_values; /* location for stashing nonzero values of matrix */
66: /* data needed for MatSOR_SeqAIJ() */
67: PetscScalar *mdiag, *idiag; /* diagonal values, inverse of diagonal entries */
68: PetscScalar *ssor_work; /* workspace for Eisenstat trick */
69: PetscObjectState idiagState; /* state of the matrix when mdiag and idiag was obtained */
70: PetscScalar fshift, omega; /* last used omega and fshift */
72: ISColoring coloring; /* set with MatADSetColoring() used by MatADSetValues() */
73: } Mat_SeqSELL;
75: /*
76: Frees the arrays from the XSELLPACK matrix type
77: */
78: static inline PetscErrorCode MatSeqXSELLFreeSELL(Mat AA, MatScalar **val, PetscInt **colidx)
79: {
80: Mat_SeqSELL *A = (Mat_SeqSELL *)AA->data;
81: if (A->singlemalloc) {
82: PetscCall(PetscFree2(*val, *colidx));
83: } else {
84: if (A->free_val) PetscCall(PetscFree(*val));
85: if (A->free_colidx) PetscCall(PetscFree(*colidx));
86: }
87: return PETSC_SUCCESS;
88: }
90: #define MatSeqXSELLReallocateSELL(Amat, AM, BS2, WIDTH, SIDX, SH, SID, ROW, COL, COLIDX, VAL, CP, VP, NONEW, datatype, MUL) \
91: do { \
92: if (WIDTH >= (SIDX[SID + 1] - SIDX[SID]) / SH) { \
93: Mat_SeqSELL *Ain = (Mat_SeqSELL *)Amat->data; \
94: /* there is no extra room in row, therefore enlarge 1 slice column */ \
95: PetscInt new_size = Ain->maxallocmat + SH * MUL, *new_colidx; \
96: datatype *new_val; \
97: \
98: PetscCheck(NONEW != -2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "New nonzero at (%" PetscInt_FMT ",%" PetscInt_FMT ") caused a malloc. Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check", ROW, COL); \
99: /* malloc new storage space */ \
100: PetscCall(PetscMalloc2(BS2 * new_size, &new_val, BS2 * new_size, &new_colidx)); \
101: \
102: /* copy over old data into new slots by two steps: one step for data before the current slice and the other for the rest */ \
103: PetscCall(PetscArraycpy(new_val, VAL, SIDX[SID + 1])); \
104: PetscCall(PetscArraycpy(new_colidx, COLIDX, SIDX[SID + 1])); \
105: PetscCall(PetscArraycpy(new_val + SIDX[SID + 1] + SH * MUL, VAL + SIDX[SID + 1], SIDX[Ain->totalslices] - SIDX[SID + 1])); \
106: PetscCall(PetscArraycpy(new_colidx + SIDX[SID + 1] + SH * MUL, COLIDX + SIDX[SID + 1], SIDX[Ain->totalslices] - SIDX[SID + 1])); \
107: /* update slice_idx */ \
108: for (ii = SID + 1; ii <= Ain->totalslices; ii++) SIDX[ii] += SH * MUL; \
109: /* update pointers. Notice that they point to the FIRST position of the row */ \
110: CP = new_colidx + SIDX[SID] + (ROW % SH); \
111: VP = new_val + SIDX[SID] + (ROW % SH); \
112: /* free up old matrix storage */ \
113: PetscCall(MatSeqXSELLFreeSELL(A, &Ain->val, &Ain->colidx)); \
114: Ain->val = new_val; \
115: Ain->colidx = new_colidx; \
116: Ain->singlemalloc = PETSC_TRUE; \
117: Ain->maxallocmat = new_size; \
118: Ain->reallocs++; \
119: A->nonzerostate++; \
120: if (WIDTH >= Ain->maxallocrow) Ain->maxallocrow += MUL; \
121: if (WIDTH >= Ain->rlenmax) Ain->rlenmax++; \
122: } \
123: } while (0)
125: #define MatSetValue_SeqSELL_Private(A, row, col, value, addv, orow, ocol, cp, vp, lastcol, low, high) \
126: do { \
127: Mat_SeqSELL *a = (Mat_SeqSELL *)A->data; \
128: found = PETSC_FALSE; \
129: if (col <= lastcol) low = 0; \
130: else high = a->rlen[row]; \
131: lastcol = col; \
132: while (high - low > 5) { \
133: t = (low + high) / 2; \
134: if (*(cp + a->sliceheight * t) > col) high = t; \
135: else low = t; \
136: } \
137: for (_i = low; _i < high; _i++) { \
138: if (*(cp + a->sliceheight * _i) > col) break; \
139: if (*(cp + a->sliceheight * _i) == col) { \
140: if (addv == ADD_VALUES) *(vp + a->sliceheight * _i) += value; \
141: else *(vp + a->sliceheight * _i) = value; \
142: found = PETSC_TRUE; \
143: break; \
144: } \
145: } \
146: if (!found) { \
147: PetscCheck(a->nonew != -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Inserting a new nonzero at global row/column (%" PetscInt_FMT ", %" PetscInt_FMT ") into matrix", orow, ocol); \
148: if (a->nonew != 1 && !(value == 0.0 && a->ignorezeroentries) && a->rlen[row] >= (a->sliidx[row / a->sliceheight + 1] - a->sliidx[row / a->sliceheight]) / a->sliceheight) { \
149: /* there is no extra room in row, therefore enlarge 1 slice column */ \
150: if (a->maxallocmat < a->sliidx[a->totalslices] + a->sliceheight) { \
151: /* allocates a larger array for the XSELL matrix types; only extend the current slice by one more column. */ \
152: PetscInt new_size = a->maxallocmat + a->sliceheight, *new_colidx; \
153: MatScalar *new_val; \
154: PetscCheck(a->nonew != -2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "New nonzero at (%" PetscInt_FMT ",%" PetscInt_FMT ") caused a malloc. Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check", orow, ocol); \
155: /* malloc new storage space */ \
156: PetscCall(PetscMalloc2(new_size, &new_val, new_size, &new_colidx)); \
157: /* copy over old data into new slots by two steps: one step for data before the current slice and the other for the rest */ \
158: PetscCall(PetscArraycpy(new_val, a->val, a->sliidx[row / a->sliceheight + 1])); \
159: PetscCall(PetscArraycpy(new_colidx, a->colidx, a->sliidx[row / a->sliceheight + 1])); \
160: PetscCall(PetscArraycpy(new_val + a->sliidx[row / a->sliceheight + 1] + a->sliceheight, PetscSafePointerPlusOffset(a->val, a->sliidx[row / a->sliceheight + 1]), a->sliidx[a->totalslices] - a->sliidx[row / a->sliceheight + 1])); \
161: PetscCall(PetscArraycpy(new_colidx + a->sliidx[row / a->sliceheight + 1] + a->sliceheight, PetscSafePointerPlusOffset(a->colidx, a->sliidx[row / a->sliceheight + 1]), a->sliidx[a->totalslices] - a->sliidx[row / a->sliceheight + 1])); \
162: /* update pointers. Notice that they point to the FIRST position of the row */ \
163: cp = new_colidx + a->sliidx[row / a->sliceheight] + (row % a->sliceheight); \
164: vp = new_val + a->sliidx[row / a->sliceheight] + (row % a->sliceheight); \
165: /* free up old matrix storage */ \
166: PetscCall(MatSeqXSELLFreeSELL(A, &a->val, &a->colidx)); \
167: a->val = new_val; \
168: a->colidx = new_colidx; \
169: a->singlemalloc = PETSC_TRUE; \
170: a->maxallocmat = new_size; \
171: a->reallocs++; \
172: } else { \
173: /* no need to reallocate, just shift the following slices to create space for the added slice column */ \
174: PetscCall(PetscArraymove(a->val + a->sliidx[row / a->sliceheight + 1] + a->sliceheight, a->val + a->sliidx[row / a->sliceheight + 1], a->sliidx[a->totalslices] - a->sliidx[row / a->sliceheight + 1])); \
175: PetscCall(PetscArraymove(a->colidx + a->sliidx[row / a->sliceheight + 1] + a->sliceheight, a->colidx + a->sliidx[row / a->sliceheight + 1], a->sliidx[a->totalslices] - a->sliidx[row / a->sliceheight + 1])); \
176: } \
177: /* update slice_idx */ \
178: for (ii = row / a->sliceheight + 1; ii <= a->totalslices; ii++) a->sliidx[ii] += a->sliceheight; \
179: if (a->rlen[row] >= a->maxallocrow) a->maxallocrow++; \
180: if (a->rlen[row] >= a->rlenmax) a->rlenmax++; \
181: } \
182: /* shift up all the later entries in this row */ \
183: for (ii = a->rlen[row] - 1; ii >= _i; ii--) { \
184: *(cp + a->sliceheight * (ii + 1)) = *(cp + a->sliceheight * ii); \
185: *(vp + a->sliceheight * (ii + 1)) = *(vp + a->sliceheight * ii); \
186: } \
187: *(cp + a->sliceheight * _i) = col; \
188: *(vp + a->sliceheight * _i) = value; \
189: a->nz++; \
190: a->rlen[row]++; \
191: A->nonzerostate++; \
192: low = _i + 1; \
193: high++; \
194: } \
195: } while (0)
197: #define PetscCeilIntMacro(x, y) ((((PetscInt)(x)) / ((PetscInt)(y))) + ((((PetscInt)(x)) % ((PetscInt)(y))) ? 1 : 0))
199: PETSC_INTERN PetscErrorCode MatSeqSELLSetPreallocation_SeqSELL(Mat, PetscInt, const PetscInt[]);
200: PETSC_INTERN PetscErrorCode MatMult_SeqSELL(Mat, Vec, Vec);
201: PETSC_INTERN PetscErrorCode MatMultAdd_SeqSELL(Mat, Vec, Vec, Vec);
202: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqSELL(Mat, Vec, Vec);
203: PETSC_INTERN PetscErrorCode MatMultTransposeAdd_SeqSELL(Mat, Vec, Vec, Vec);
204: PETSC_INTERN PetscErrorCode MatZeroEntries_SeqSELL(Mat);
205: PETSC_INTERN PetscErrorCode MatDestroy_SeqSELL(Mat);
206: PETSC_INTERN PetscErrorCode MatSetOption_SeqSELL(Mat, MatOption, PetscBool);
207: PETSC_INTERN PetscErrorCode MatGetDiagonal_SeqSELL(Mat, Vec v);
208: PETSC_INTERN PetscErrorCode MatGetValues_SeqSELL(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], PetscScalar[]);
209: PETSC_INTERN PetscErrorCode MatView_SeqSELL(Mat, PetscViewer);
210: PETSC_INTERN PetscErrorCode MatAssemblyEnd_SeqSELL(Mat, MatAssemblyType);
211: PETSC_INTERN PetscErrorCode MatGetInfo_SeqSELL(Mat, MatInfoType, MatInfo *);
212: PETSC_INTERN PetscErrorCode MatSetValues_SeqSELL(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode);
213: PETSC_INTERN PetscErrorCode MatCopy_SeqSELL(Mat, Mat, MatStructure);
214: PETSC_INTERN PetscErrorCode MatSetUp_SeqSELL(Mat);
215: PETSC_INTERN PetscErrorCode MatSeqSELLGetArray_SeqSELL(Mat, PetscScalar *[]);
216: PETSC_INTERN PetscErrorCode MatSeqSELLRestoreArray_SeqSELL(Mat, PetscScalar *[]);
217: PETSC_INTERN PetscErrorCode MatShift_SeqSELL(Mat, PetscScalar);
218: PETSC_INTERN PetscErrorCode MatSOR_SeqSELL(Mat, Vec, PetscReal, MatSORType, PetscReal, PetscInt, PetscInt, Vec);
219: PETSC_EXTERN PetscErrorCode MatCreate_SeqSELL(Mat);
220: PETSC_INTERN PetscErrorCode MatDuplicate_SeqSELL(Mat, MatDuplicateOption, Mat *);
221: PETSC_INTERN PetscErrorCode MatEqual_SeqSELL(Mat, Mat, PetscBool *);
222: PETSC_INTERN PetscErrorCode MatConvert_SeqSELL_SeqAIJ(Mat, MatType, MatReuse, Mat *);
223: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat, MatType, MatReuse, Mat *);
224: PETSC_INTERN PetscErrorCode MatFDColoringCreate_SeqSELL(Mat, ISColoring, MatFDColoring);
225: PETSC_INTERN PetscErrorCode MatFDColoringSetUp_SeqSELL(Mat, ISColoring, MatFDColoring);
226: PETSC_INTERN PetscErrorCode MatGetColumnIJ_SeqSELL_Color(Mat, PetscInt, PetscBool, PetscBool, PetscInt *, const PetscInt *[], const PetscInt *[], PetscInt *[], PetscBool *);
227: PETSC_INTERN PetscErrorCode MatRestoreColumnIJ_SeqSELL_Color(Mat, PetscInt, PetscBool, PetscBool, PetscInt *, const PetscInt *[], const PetscInt *[], PetscInt *[], PetscBool *);
228: PETSC_INTERN PetscErrorCode MatConjugate_SeqSELL(Mat A);
229: PETSC_INTERN PetscErrorCode MatScale_SeqSELL(Mat, PetscScalar);
230: PETSC_INTERN PetscErrorCode MatDiagonalScale_SeqSELL(Mat, Vec, Vec);