1 2 #ifndef __SELL_H 3 #define __SELL_H 4 5 #include <petsc/private/matimpl.h> 6 #include <petsc/private/hashmapi.h> 7 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 */ 42 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; 51 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 0; 65 } 66 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 postion 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 } 98 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 postion 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 } 170 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 208