Actual source code: sell.h


  2: #ifndef __SELL_H

  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: static inline PetscErrorCode MatSeqXSELLFreeSELL(Mat AA, MatScalar **val, PetscInt **colidx)
 56: {
 57:   Mat_SeqSELL *A = (Mat_SeqSELL *)AA->data;
 58:   if (A->singlemalloc) {
 59:     PetscFree2(*val, *colidx);
 60:   } else {
 61:     if (A->free_val) PetscFree(*val);
 62:     if (A->free_colidx) PetscFree(*colidx);
 63:   }
 64:   return 0;
 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: \
 75:     /* malloc new storage space */ \
 76:     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:     PetscArraycpy(new_val, VAL, SIDX[SID + 1]); \
 80:     PetscArraycpy(new_colidx, COLIDX, SIDX[SID + 1]); \
 81:     PetscArraycpy(new_val + SIDX[SID + 1] + 8, VAL + SIDX[SID + 1], SIDX[AM >> 3] - SIDX[SID + 1]); \
 82:     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:     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) { \
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; \
129:           /* malloc new storage space */ \
130:           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:           PetscArraycpy(new_val, a->val, a->sliidx[row / 8 + 1]); \
133:           PetscArraycpy(new_colidx, a->colidx, a->sliidx[row / 8 + 1]); \
134:           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:           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:           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:           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:           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