Actual source code: matpreallocator.c

petsc-3.14.6 2021-03-30
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, NULL);
 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:   /* Do not bother bstash since MatPreallocator does not implement MatSetValuesBlocked */
 38:   MatStashCreate_Private(PetscObjectComm((PetscObject) A), 1, &A->stash);
 39:   /* arrays are for blocked rows/cols */
 40:   mbs  = m/bs;
 41:   PetscCalloc4(mbs, &p->dnz, mbs, &p->onz, mbs, &p->dnzu, mbs, &p->onzu);
 42:   return(0);
 43: }

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

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

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

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

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

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

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

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

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

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

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

144:   MatGetBlockSize(mat, &bs);
145:   MatXAIJSetPreallocation(A, bs, p->dnz, p->onz, p->dnzu, p->onzu);
146:   MatSetUp(A);
147:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
148:   if (fill) {
149:     PetscHashIter  hi;
150:     PetscHashIJKey key;
151:     PetscScalar    *zeros;

153:     PetscCalloc1(bs*bs,&zeros);

155:     PetscHashIterBegin(p->ht,hi);
156:     while (!PetscHashIterAtEnd(p->ht,hi)) {
157:       PetscHashIterGetKey(p->ht,hi,key);
158:       PetscHashIterNext(p->ht,hi);
159:       MatSetValuesBlocked(A,1,&key.i,1,&key.j,zeros,INSERT_VALUES);
160:     }
161:     PetscFree(zeros);

163:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
164:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
165:   }
166:   return(0);
167: }

169: /*@
170:   MatPreallocatorPreallocate - Preallocates the A matrix, using information from mat, optionally filling A with zeros

172:   Input Parameter:
173: + mat  - the preallocator
174: . fill - fill the matrix with zeros
175: - A    - the matrix to be preallocated

177:   Notes:
178:   This Mat implementation provides a helper utility to define the correct
179:   preallocation data for a given nonzero structure. Use this object like a
180:   regular matrix, e.g. loop over the nonzero structure of the matrix and
181:   call MatSetValues() or MatSetValuesBlocked() to indicate the nonzero locations.
182:   The matrix entires provided to MatSetValues() will be ignored, it only uses
183:   the row / col indices provided to determine the information required to be
184:   passed to MatXAIJSetPreallocation(). Once you have looped over the nonzero
185:   structure, you must call MatAssemblyBegin(), MatAssemblyEnd() on mat.

187:   After you have assembled the preallocator matrix (mat), call MatPreallocatorPreallocate()
188:   to define the preallocation information on the matrix (A). Setting the parameter
189:   fill = PETSC_TRUE will insert zeros into the matrix A. Internally MatPreallocatorPreallocate()
190:   will call MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);

192:   Level: advanced

194: .seealso: MATPREALLOCATOR
195: @*/
196: PetscErrorCode MatPreallocatorPreallocate(Mat mat, PetscBool fill, Mat A)
197: {

203:   PetscUseMethod(mat, "MatPreallocatorPreallocate_C", (Mat,PetscBool,Mat),(mat,fill,A));
204:   return(0);
205: }

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

210:    Operations Provided:
211: .  MatSetValues()

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

216:   Level: advanced

218: .seealso: Mat

220: M*/

222: PETSC_EXTERN PetscErrorCode MatCreate_Preallocator(Mat A)
223: {
224:   Mat_Preallocator *p;
225:   PetscErrorCode    ierr;

228:   PetscNewLog(A, &p);
229:   A->data = (void *) p;

231:   p->ht   = NULL;
232:   p->dnz  = NULL;
233:   p->onz  = NULL;
234:   p->dnzu = NULL;
235:   p->onzu = NULL;

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

240:   A->ops->destroy       = MatDestroy_Preallocator;
241:   A->ops->setup         = MatSetUp_Preallocator;
242:   A->ops->setvalues     = MatSetValues_Preallocator;
243:   A->ops->assemblybegin = MatAssemblyBegin_Preallocator;
244:   A->ops->assemblyend   = MatAssemblyEnd_Preallocator;
245:   A->ops->view          = MatView_Preallocator;
246:   A->ops->setoption     = MatSetOption_Preallocator;
247:   A->ops->setblocksizes = MatSetBlockSizes_Default; /* once set, user is not allowed to change the block sizes */

249:   /* special MATPREALLOCATOR functions */
250:   PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", MatPreallocatorPreallocate_Preallocator);
251:   PetscObjectChangeTypeName((PetscObject) A, MATPREALLOCATOR);
252:   return(0);
253: }