Actual source code: sell.h
petsc-3.9.4 2018-09-11
5: #include <petsc/private/matimpl.h>
6: #include <petscctable.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: PETSC_STATIC_INLINE PetscErrorCode MatSeqXSELLFreeSELL(Mat AA,MatScalar **val,PetscInt **colidx)
56: {
57: Mat_SeqSELL *A = (Mat_SeqSELL*) AA->data;
59: if (A->singlemalloc) {
60: PetscFree2(*val,*colidx);
61: } else {
62: if (A->free_val) {PetscFree(*val);}
63: if (A->free_colidx) {PetscFree(*colidx);}
64: }
65: return 0;
66: }
68: #define MatSeqXSELLReallocateSELL(Amat,AM,BS2,WIDTH,SIDX,SID,ROW,COL,COLIDX,VAL,CP,VP,NONEW,datatype) \
69: if (WIDTH >= (SIDX[SID+1]-SIDX[SID])/8) { \
70: Mat_SeqSELL *Ain = (Mat_SeqSELL*)Amat->data; \
71: /* there is no extra room in row, therefore enlarge 8 elements (1 slice column) */ \
72: PetscInt new_size=Ain->maxallocmat+8,*new_colidx; \
73: datatype *new_val; \
74: \
75: 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); \
76: /* malloc new storage space */ \
77: PetscMalloc2(BS2*new_size,&new_val,BS2*new_size,&new_colidx); \
78: \
79: /* copy over old data into new slots by two steps: one step for data before the current slice and the other for the rest */ \
80: PetscMemcpy(new_val,VAL,SIDX[SID+1]*sizeof(datatype)); \
81: PetscMemcpy(new_colidx,COLIDX,SIDX[SID+1]*sizeof(PetscInt)); \
82: PetscMemcpy(new_val+SIDX[SID+1]+8,VAL+SIDX[SID+1],(SIDX[AM>>3]-SIDX[SID+1])*sizeof(datatype)); \
83: PetscMemcpy(new_colidx+SIDX[SID+1]+8,COLIDX+SIDX[SID+1],(SIDX[AM>>3]-SIDX[SID+1])*sizeof(PetscInt)); \
84: /* update slice_idx */ \
85: for (ii=SID+1;ii<=AM>>3;ii++) { SIDX[ii] += 8; } \
86: /* update pointers. Notice that they point to the FIRST postion of the row */ \
87: CP = new_colidx+SIDX[SID]+(ROW & 0x07); \
88: VP = new_val+SIDX[SID]+(ROW & 0x07); \
89: /* free up old matrix storage */ \
90: MatSeqXSELLFreeSELL(A,&Ain->val,&Ain->colidx); \
91: Ain->val = (MatScalar*) new_val; \
92: Ain->colidx = new_colidx; \
93: Ain->singlemalloc = PETSC_TRUE; \
94: Ain->maxallocmat = new_size; \
95: Ain->reallocs++; \
96: if (WIDTH>=Ain->maxallocrow) Ain->maxallocrow++; \
97: if (WIDTH>=Ain->rlenmax) Ain->rlenmax++; \
98: } \
100: #define MatSetValue_SeqSELL_Private(A,row,col,value,addv,orow,ocol,cp,vp,lastcol,low,high) \
101: { \
102: Mat_SeqSELL *a=(Mat_SeqSELL*)A->data; \
103: found=PETSC_FALSE; \
104: if (col <= lastcol) low = 0; \
105: else high = a->rlen[row]; \
106: lastcol = col; \
107: while (high-low > 5) { \
108: t = (low+high)/2; \
109: if (*(cp+8*t) > col) high = t; \
110: else low = t; \
111: } \
112: for (_i=low; _i<high; _i++) { \
113: if (*(cp+8*_i) > col) break; \
114: if (*(cp+8*_i) == col) { \
115: if (addv == ADD_VALUES)*(vp+8*_i) += value; \
116: else *(vp+8*_i) = value; \
117: found = PETSC_TRUE; \
118: break; \
119: } \
120: } \
121: if (!found) { \
122: 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); \
123: if (a->nonew != 1 && !(value == 0.0 && a->ignorezeroentries) && a->rlen[row] >= (a->sliidx[row/8+1]-a->sliidx[row/8])/8) { \
124: /* there is no extra room in row, therefore enlarge 8 elements (1 slice column) */ \
125: if (a->maxallocmat < a->sliidx[a->totalslices]+8) { \
126: /* allocates a larger array for the XSELL matrix types; only extend the current slice by one more column. */ \
127: PetscInt new_size=a->maxallocmat+8,*new_colidx; \
128: MatScalar *new_val; \
129: 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); \
130: /* malloc new storage space */ \
131: PetscMalloc2(new_size,&new_val,new_size,&new_colidx); \
132: /* copy over old data into new slots by two steps: one step for data before the current slice and the other for the rest */ \
133: PetscMemcpy(new_val,a->val,a->sliidx[row/8+1]*sizeof(MatScalar)); \
134: PetscMemcpy(new_colidx,a->colidx,a->sliidx[row/8+1]*sizeof(PetscInt)); \
135: 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)); \
136: 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)); \
137: /* update pointers. Notice that they point to the FIRST postion of the row */ \
138: cp = new_colidx+a->sliidx[row/8]+(row & 0x07); \
139: vp = new_val+a->sliidx[row/8]+(row & 0x07); \
140: /* free up old matrix storage */ \
141: MatSeqXSELLFreeSELL(A,&a->val,&a->colidx); \
142: a->val = (MatScalar*)new_val; \
143: a->colidx = new_colidx; \
144: a->singlemalloc = PETSC_TRUE; \
145: a->maxallocmat = new_size; \
146: a->reallocs++; \
147: } else { \
148: /* no need to reallocate, just shift the following slices to create space for the added slice column */ \
149: 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)); \
150: 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)); \
151: } \
152: /* update slice_idx */ \
153: for (ii=row/8+1;ii<=a->totalslices;ii++) a->sliidx[ii] += 8; \
154: if (a->rlen[row]>=a->maxallocrow) a->maxallocrow++; \
155: if (a->rlen[row]>=a->rlenmax) a->rlenmax++; \
156: } \
157: /* shift up all the later entries in this row */ \
158: for (ii=a->rlen[row]-1; ii>=_i; ii--) { \
159: *(cp+8*(ii+1)) = *(cp+8*ii); \
160: *(vp+8*(ii+1)) = *(vp+8*ii); \
161: } \
162: *(cp+8*_i) = col; \
163: *(vp+8*_i) = value; \
164: a->nz++; a->rlen[row]++; A->nonzerostate++; \
165: low = _i+1; high++; \
166: } \
167: } \
169: PETSC_INTERN PetscErrorCode MatSeqSELLSetPreallocation_SeqSELL(Mat,PetscInt,const PetscInt[]);
170: PETSC_INTERN PetscErrorCode MatMult_SeqSELL(Mat,Vec,Vec);
171: PETSC_INTERN PetscErrorCode MatMultAdd_SeqSELL(Mat,Vec,Vec,Vec);
172: PETSC_INTERN PetscErrorCode MatMultTranspose_SeqSELL(Mat,Vec,Vec);
173: PETSC_INTERN PetscErrorCode MatMultTransposeAdd_SeqSELL(Mat,Vec,Vec,Vec);
174: PETSC_INTERN PetscErrorCode MatMissingDiagonal_SeqSELL(Mat,PetscBool*,PetscInt*);
175: PETSC_INTERN PetscErrorCode MatMarkDiagonal_SeqSELL(Mat);
176: PETSC_INTERN PetscErrorCode MatInvertDiagonal_SeqSELL(Mat,PetscScalar,PetscScalar);
177: PETSC_INTERN PetscErrorCode MatZeroEntries_SeqSELL(Mat);
178: PETSC_INTERN PetscErrorCode MatDestroy_SeqSELL(Mat);
179: PETSC_INTERN PetscErrorCode MatSetOption_SeqSELL(Mat,MatOption,PetscBool);
180: PETSC_INTERN PetscErrorCode MatGetDiagonal_SeqSELL(Mat,Vec v);
181: PETSC_INTERN PetscErrorCode MatGetValues_SeqSELL(Mat,PetscInt,const PetscInt [],PetscInt,const PetscInt[],PetscScalar[]);
182: PETSC_INTERN PetscErrorCode MatView_SeqSELL(Mat,PetscViewer);
183: PETSC_INTERN PetscErrorCode MatAssemblyEnd_SeqSELL(Mat,MatAssemblyType);
184: PETSC_INTERN PetscErrorCode MatGetInfo_SeqSELL(Mat,MatInfoType,MatInfo*);
185: PETSC_INTERN PetscErrorCode MatSetValues_SeqSELL(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode);
186: PETSC_INTERN PetscErrorCode MatCopy_SeqSELL(Mat,Mat,MatStructure);
187: PETSC_INTERN PetscErrorCode MatSetUp_SeqSELL(Mat);
188: PETSC_INTERN PetscErrorCode MatSeqSELLGetArray_SeqSELL(Mat,PetscScalar *[]);
189: PETSC_INTERN PetscErrorCode MatSeqSELLRestoreArray_SeqSELL(Mat,PetscScalar *[]);
190: PETSC_INTERN PetscErrorCode MatShift_SeqSELL(Mat,PetscScalar);
191: PETSC_INTERN PetscErrorCode MatSOR_SeqSELL(Mat,Vec,PetscReal,MatSORType,PetscReal,PetscInt,PetscInt,Vec);
192: PETSC_EXTERN PetscErrorCode MatCreate_SeqSELL(Mat);
193: PETSC_INTERN PetscErrorCode MatDuplicate_SeqSELL(Mat,MatDuplicateOption,Mat*);
194: PETSC_INTERN PetscErrorCode MatEqual_SeqSELL(Mat,Mat,PetscBool*);
195: PETSC_INTERN PetscErrorCode MatSeqSELLInvalidateDiagonal(Mat);
196: PETSC_INTERN PetscErrorCode MatConvert_SeqSELL_SeqAIJ(Mat,MatType,MatReuse,Mat*);
197: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*);
198: PETSC_INTERN PetscErrorCode MatFDColoringCreate_SeqSELL(Mat,ISColoring,MatFDColoring);
199: PETSC_INTERN PetscErrorCode MatFDColoringSetUp_SeqSELL(Mat,ISColoring,MatFDColoring);
200: PETSC_INTERN PetscErrorCode MatGetColumnIJ_SeqSELL_Color(Mat,PetscInt,PetscBool,PetscBool,PetscInt*,const PetscInt *[],const PetscInt *[],PetscInt *[],PetscBool*);
201: PETSC_INTERN PetscErrorCode MatRestoreColumnIJ_SeqSELL_Color(Mat,PetscInt,PetscBool,PetscBool,PetscInt*,const PetscInt *[],const PetscInt *[],PetscInt *[],PetscBool*);
202: PETSC_INTERN PetscErrorCode MatConjugate_SeqSELL(Mat A);
203: PETSC_INTERN PetscErrorCode MatScale_SeqSELL(Mat,PetscScalar);
204: PETSC_INTERN PetscErrorCode MatDiagonalScale_SeqSELL(Mat,Vec,Vec);
205: #endif