Actual source code: matpreallocator.c

petsc-3.11.2 2019-05-18
Report Typos and Errors
  1:  #include <petsc/private/matimpl.h>
  2:  #include <petsc/private/hashsetij.h>

  4: typedef struct {
  5:   PetscHSetIJ ht;
  6:   PetscInt   *dnz, *onz;
  7:   PetscInt   *dnzu, *onzu;
  8: } Mat_Preallocator;

 10: PetscErrorCode MatDestroy_Preallocator(Mat A)
 11: {
 12:   Mat_Preallocator *p = (Mat_Preallocator *) A->data;
 13:   PetscErrorCode    ierr;

 16:   MatStashDestroy_Private(&A->stash);
 17:   PetscHSetIJDestroy(&p->ht);
 18:   PetscFree4(p->dnz, p->onz, p->dnzu, p->onzu);
 19:   PetscFree(A->data);
 20:   PetscObjectChangeTypeName((PetscObject) A, 0);
 21:   PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", NULL);
 22:   return(0);
 23: }

 25: PetscErrorCode MatSetUp_Preallocator(Mat A)
 26: {
 27:   Mat_Preallocator *p = (Mat_Preallocator *) A->data;
 28:   PetscInt          m, bs, mbs;
 29:   PetscErrorCode    ierr;

 32:   PetscLayoutSetUp(A->rmap);
 33:   PetscLayoutSetUp(A->cmap);
 34:   MatGetLocalSize(A, &m, NULL);
 35:   PetscHSetIJCreate(&p->ht);
 36:   MatGetBlockSize(A, &bs);
 37:   MatStashCreate_Private(PetscObjectComm((PetscObject) A), bs, &A->stash);
 38:   /* arrays are for blocked rows/cols */
 39:   mbs  = m/bs;
 40:   PetscCalloc4(mbs, &p->dnz, mbs, &p->onz, mbs, &p->dnzu, mbs, &p->onzu);
 41:   return(0);
 42: }

 44: PetscErrorCode MatSetValues_Preallocator(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
 45: {
 46:   Mat_Preallocator *p = (Mat_Preallocator *) A->data;
 47:   PetscInt          rStart, rEnd, r, cStart, cEnd, c, bs;
 48:   PetscErrorCode    ierr;

 51:   MatGetBlockSize(A, &bs);
 52:   MatGetOwnershipRange(A, &rStart, &rEnd);
 53:   MatGetOwnershipRangeColumn(A, &cStart, &cEnd);
 54:   for (r = 0; r < m; ++r) {
 55:     PetscHashIJKey key;
 56:     PetscBool      missing;

 58:     key.i = rows[r];
 59:     if (key.i < 0) continue;
 60:     if ((key.i < rStart) || (key.i >= rEnd)) {
 61:       MatStashValuesRow_Private(&A->stash, key.i, n, cols, values, PETSC_FALSE);
 62:     } else { /* Hash table is for blocked rows/cols */
 63:       key.i = rows[r]/bs;
 64:       for (c = 0; c < n; ++c) {
 65:         key.j = cols[c]/bs;
 66:         if (key.j < 0) continue;
 67:         PetscHSetIJQueryAdd(p->ht, key, &missing);
 68:         if (missing) {
 69:           if ((key.j >= cStart/bs) && (key.j < cEnd/bs)) {
 70:             ++p->dnz[key.i-rStart/bs];
 71:             if (key.j >= key.i) ++p->dnzu[key.i-rStart/bs];
 72:           } else {
 73:             ++p->onz[key.i-rStart/bs];
 74:             if (key.j >= key.i) ++p->onzu[key.i-rStart/bs];
 75:           }
 76:         }
 77:       }
 78:     }
 79:   }
 80:   return(0);
 81: }

 83: PetscErrorCode MatAssemblyBegin_Preallocator(Mat A, MatAssemblyType type)
 84: {
 85:   PetscInt       nstash, reallocs;

 89:   MatStashScatterBegin_Private(A, &A->stash, A->rmap->range);
 90:   MatStashGetInfo_Private(&A->stash, &nstash, &reallocs);
 91:   PetscInfo2(A, "Stash has %D entries, uses %D mallocs.\n", nstash, reallocs);
 92:   return(0);
 93: }

 95: PetscErrorCode MatAssemblyEnd_Preallocator(Mat A, MatAssemblyType type)
 96: {
 97:   PetscScalar   *val;
 98:   PetscInt      *row, *col;
 99:   PetscInt       i, j, rstart, ncols, flg;
100:   PetscMPIInt    n;

104:   while (1) {
105:     MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg);
106:     if (!flg) break;

108:     for (i = 0; i < n; ) {
109:       /* Now identify the consecutive vals belonging to the same row */
110:       for (j = i, rstart = row[j]; j < n; j++) {
111:         if (row[j] != rstart) break;
112:       }
113:       if (j < n) ncols = j-i;
114:       else       ncols = n-i;
115:       /* Now assemble all these values with a single function call */
116:       MatSetValues_Preallocator(A, 1, row+i, ncols, col+i, val+i, INSERT_VALUES);
117:       i = j;
118:     }
119:   }
120:   MatStashScatterEnd_Private(&A->stash);
121:   return(0);
122: }

124: PetscErrorCode MatView_Preallocator(Mat A, PetscViewer viewer)
125: {
127:   return(0);
128: }

130: PetscErrorCode MatSetOption_Preallocator(Mat A, MatOption op, PetscBool flg)
131: {
133:   return(0);
134: }

136: PetscErrorCode MatPreallocatorPreallocate_Preallocator(Mat mat, PetscBool fill, Mat A)
137: {
138:   Mat_Preallocator *p = (Mat_Preallocator *) mat->data;
139:   PetscInt          bs;
140:   PetscErrorCode    ierr;

143:   MatGetBlockSize(mat, &bs);
144:   MatXAIJSetPreallocation(A, bs, p->dnz, p->onz, p->dnzu, p->onzu);
145:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
146:   return(0);
147: }

149: /*@
150:   MatPreallocatorPreallocate - Preallocates the input matrix, optionally filling it with zeros

152:   Input Parameter:
153: + mat  - the preallocator
154: - fill - fill the matrix with zeros

156:   Output Parameter:
157: . A    - the matrix

159:   Level: advanced

161: .seealso: MATPREALLOCATOR
162: @*/
163: PetscErrorCode MatPreallocatorPreallocate(Mat mat, PetscBool fill, Mat A)
164: {

170:   PetscUseMethod(mat, "MatPreallocatorPreallocate_C", (Mat,PetscBool,Mat),(mat,fill,A));
171:   return(0);
172: }

174: /*MC
175:    MATPREALLOCATOR - MATPREALLOCATOR = "preallocator" - A matrix type to be used for computing a matrix preallocation.

177:    Operations Provided:
178: .  MatSetValues()

180:    Options Database Keys:
181: . -mat_type preallocator - sets the matrix type to "preallocator" during a call to MatSetFromOptions()

183:   Level: advanced

185: .seealso: Mat

187: M*/

189: PETSC_EXTERN PetscErrorCode MatCreate_Preallocator(Mat A)
190: {
191:   Mat_Preallocator *p;
192:   PetscErrorCode    ierr;

195:   PetscNewLog(A, &p);
196:   A->data = (void *) p;

198:   p->ht   = NULL;
199:   p->dnz  = NULL;
200:   p->onz  = NULL;
201:   p->dnzu = NULL;
202:   p->onzu = NULL;

204:   /* matrix ops */
205:   PetscMemzero(A->ops, sizeof(struct _MatOps));

207:   A->ops->destroy       = MatDestroy_Preallocator;
208:   A->ops->setup         = MatSetUp_Preallocator;
209:   A->ops->setvalues     = MatSetValues_Preallocator;
210:   A->ops->assemblybegin = MatAssemblyBegin_Preallocator;
211:   A->ops->assemblyend   = MatAssemblyEnd_Preallocator;
212:   A->ops->view          = MatView_Preallocator;
213:   A->ops->setoption     = MatSetOption_Preallocator;
214:   A->ops->setblocksizes = MatSetBlockSizes_Default; /* once set, user is not allowed to change the block sizes */

216:   /* special MATPREALLOCATOR functions */
217:   PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", MatPreallocatorPreallocate_Preallocator);
218:   PetscObjectChangeTypeName((PetscObject) A, MATPREALLOCATOR);
219:   return(0);
220: }