1 2 #if !defined(__SELL_H) 3 #define __SELL_H 4 5 #include <petsc/private/matimpl.h> 6 #include <petscctable.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 *ibdiag; /* inverses of block diagonals */ 49 PetscBool ibdiagvalid; /* inverses of block diagonals are valid. */ 50 PetscScalar fshift,omega; /* last used omega and fshift */ 51 ISColoring coloring; /* set with MatADSetColoring() used by MatADSetValues() */ 52 } Mat_SeqSELL; 53 54 /* 55 Frees the arrays from the XSELLPACK matrix type 56 */ 57 PETSC_STATIC_INLINE PetscErrorCode MatSeqXSELLFreeSELL(Mat AA,MatScalar **val,PetscInt **colidx) 58 { 59 Mat_SeqSELL *A = (Mat_SeqSELL*) AA->data; 60 PetscErrorCode ierr; 61 if (A->singlemalloc) { 62 ierr = PetscFree2(*val,*colidx);CHKERRQ(ierr); 63 } else { 64 if (A->free_val) {ierr = PetscFree(*val);CHKERRQ(ierr);} 65 if (A->free_colidx) {ierr = PetscFree(*colidx);CHKERRQ(ierr);} 66 } 67 return 0; 68 } 69 70 #define MatSeqXSELLReallocateSELL(Amat,AM,BS2,WIDTH,SIDX,SID,ROW,COL,COLIDX,VAL,CP,VP,NONEW,datatype) \ 71 if (WIDTH >= (SIDX[SID+1]-SIDX[SID])/8) { \ 72 Mat_SeqSELL *Ain = (Mat_SeqSELL*)Amat->data; \ 73 /* there is no extra room in row, therefore enlarge 8 elements (1 slice column) */ \ 74 PetscInt new_size=Ain->maxallocmat+8,*new_colidx; \ 75 datatype *new_val; \ 76 \ 77 if (NONEW == -2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at (%D,%D) caused a malloc\nUse MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check",ROW,COL); \ 78 /* malloc new storage space */ \ 79 ierr = PetscMalloc2(BS2*new_size,&new_val,BS2*new_size,&new_colidx);CHKERRQ(ierr); \ 80 \ 81 /* copy over old data into new slots by two steps: one step for data before the current slice and the other for the rest */ \ 82 ierr = PetscMemcpy(new_val,VAL,SIDX[SID+1]*sizeof(datatype));CHKERRQ(ierr); \ 83 ierr = PetscMemcpy(new_colidx,COLIDX,SIDX[SID+1]*sizeof(PetscInt));CHKERRQ(ierr); \ 84 ierr = PetscMemcpy(new_val+SIDX[SID+1]+8,VAL+SIDX[SID+1],(SIDX[AM>>3]-SIDX[SID+1])*sizeof(datatype));CHKERRQ(ierr); \ 85 ierr = PetscMemcpy(new_colidx+SIDX[SID+1]+8,COLIDX+SIDX[SID+1],(SIDX[AM>>3]-SIDX[SID+1])*sizeof(PetscInt));CHKERRQ(ierr); \ 86 /* update slice_idx */ \ 87 for (ii=SID+1;ii<=AM>>3;ii++) { SIDX[ii] += 8; } \ 88 /* update pointers. Notice that they point to the FIRST postion of the row */ \ 89 CP = new_colidx+SIDX[SID]+(ROW & 0x07); \ 90 VP = new_val+SIDX[SID]+(ROW & 0x07); \ 91 /* free up old matrix storage */ \ 92 ierr = MatSeqXSELLFreeSELL(A,&Ain->val,&Ain->colidx);CHKERRQ(ierr); \ 93 Ain->val = (MatScalar*) new_val; \ 94 Ain->colidx = new_colidx; \ 95 Ain->singlemalloc = PETSC_TRUE; \ 96 Ain->maxallocmat = new_size; \ 97 Ain->reallocs++; \ 98 if (WIDTH>=Ain->maxallocrow) Ain->maxallocrow++; \ 99 if (WIDTH>=Ain->rlenmax) Ain->rlenmax++; \ 100 } \ 101 102 #define MatSetValue_SeqSELL_Private(A,row,col,value,addv,orow,ocol,cp,vp,lastcol,low,high) \ 103 { \ 104 Mat_SeqSELL *a=(Mat_SeqSELL*)A->data; \ 105 found=PETSC_FALSE; \ 106 if (col <= lastcol) low = 0; \ 107 else high = a->rlen[row]; \ 108 lastcol = col; \ 109 while (high-low > 5) { \ 110 t = (low+high)/2; \ 111 if (*(cp+8*t) > col) high = t; \ 112 else low = t; \ 113 } \ 114 for (_i=low; _i<high; _i++) { \ 115 if (*(cp+8*_i) > col) break; \ 116 if (*(cp+8*_i) == col) { \ 117 if (addv == ADD_VALUES)*(vp+8*_i) += value; \ 118 else *(vp+8*_i) = value; \ 119 found = PETSC_TRUE; \ 120 break; \ 121 } \ 122 } \ 123 if (!found) { \ 124 if (a->nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at global row/column (%D, %D) into matrix", orow, ocol); \ 125 if (a->nonew != 1 && !(value == 0.0 && a->ignorezeroentries) && a->rlen[row] >= (a->sliidx[row/8+1]-a->sliidx[row/8])/8) { \ 126 /* there is no extra room in row, therefore enlarge 8 elements (1 slice column) */ \ 127 if (a->maxallocmat < a->sliidx[a->totalslices]+8) { \ 128 /* allocates a larger array for the XSELL matrix types; only extend the current slice by one more column. */ \ 129 PetscInt new_size=a->maxallocmat+8,*new_colidx; \ 130 MatScalar *new_val; \ 131 if (a->nonew == -2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at (%D,%D) caused a malloc\nUse MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check",orow,ocol); \ 132 /* malloc new storage space */ \ 133 ierr = PetscMalloc2(new_size,&new_val,new_size,&new_colidx);CHKERRQ(ierr); \ 134 /* copy over old data into new slots by two steps: one step for data before the current slice and the other for the rest */ \ 135 ierr = PetscMemcpy(new_val,a->val,a->sliidx[row/8+1]*sizeof(MatScalar));CHKERRQ(ierr); \ 136 ierr = PetscMemcpy(new_colidx,a->colidx,a->sliidx[row/8+1]*sizeof(PetscInt));CHKERRQ(ierr); \ 137 ierr = PetscMemcpy(new_val+a->sliidx[row/8+1]+8,a->val+a->sliidx[row/8+1],(a->sliidx[a->totalslices]-a->sliidx[row/8+1])*sizeof(MatScalar));CHKERRQ(ierr); \ 138 ierr = PetscMemcpy(new_colidx+a->sliidx[row/8+1]+8,a->colidx+a->sliidx[row/8+1],(a->sliidx[a->totalslices]-a->sliidx[row/8+1])*sizeof(PetscInt));CHKERRQ(ierr); \ 139 /* update pointers. Notice that they point to the FIRST postion of the row */ \ 140 cp = new_colidx+a->sliidx[row/8]+(row & 0x07); \ 141 vp = new_val+a->sliidx[row/8]+(row & 0x07); \ 142 /* free up old matrix storage */ \ 143 ierr = MatSeqXSELLFreeSELL(A,&a->val,&a->colidx);CHKERRQ(ierr); \ 144 a->val = (MatScalar*)new_val; \ 145 a->colidx = new_colidx; \ 146 a->singlemalloc = PETSC_TRUE; \ 147 a->maxallocmat = new_size; \ 148 a->reallocs++; \ 149 } else { \ 150 /* no need to reallocate, just shift the following slices to create space for the added slice column */ \ 151 ierr = PetscMemmove(a->val+a->sliidx[row/8+1]+8,a->val+a->sliidx[row/8+1],(a->sliidx[a->totalslices]-a->sliidx[row/8+1])*sizeof(MatScalar));CHKERRQ(ierr); \ 152 ierr = PetscMemmove(a->colidx+a->sliidx[row/8+1]+8,a->colidx+a->sliidx[row/8+1],(a->sliidx[a->totalslices]-a->sliidx[row/8+1])*sizeof(PetscInt));CHKERRQ(ierr); \ 153 } \ 154 /* update slice_idx */ \ 155 for (ii=row/8+1;ii<=a->totalslices;ii++) a->sliidx[ii] += 8; \ 156 if (a->rlen[row]>=a->maxallocrow) a->maxallocrow++; \ 157 if (a->rlen[row]>=a->rlenmax) a->rlenmax++; \ 158 } \ 159 /* shift up all the later entries in this row */ \ 160 for (ii=a->rlen[row]-1; ii>=_i; ii--) { \ 161 *(cp+8*(ii+1)) = *(cp+8*ii); \ 162 *(vp+8*(ii+1)) = *(vp+8*ii); \ 163 } \ 164 *(cp+8*_i) = col; \ 165 *(vp+8*_i) = value; \ 166 a->nz++; a->rlen[row]++; A->nonzerostate++; \ 167 low = _i+1; 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