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