Actual source code: matpreallocator.c

  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:   PetscBool   nooffproc;
  9:   PetscBool   used;
 10: } Mat_Preallocator;

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

 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;

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

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

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

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

 80: PetscErrorCode MatAssemblyBegin_Preallocator(Mat A, MatAssemblyType type)
 81: {
 82:   PetscInt       nstash, reallocs;

 84:   MatStashScatterBegin_Private(A, &A->stash, A->rmap->range);
 85:   MatStashGetInfo_Private(&A->stash, &nstash, &reallocs);
 86:   PetscInfo(A, "Stash has %" PetscInt_FMT " entries, uses %" PetscInt_FMT " mallocs.\n", nstash, reallocs);
 87:   return 0;
 88: }

 90: PetscErrorCode MatAssemblyEnd_Preallocator(Mat A, MatAssemblyType type)
 91: {
 92:   PetscScalar      *val;
 93:   PetscInt         *row, *col;
 94:   PetscInt         i, j, rstart, ncols, flg;
 95:   PetscMPIInt      n;
 96:   Mat_Preallocator *p = (Mat_Preallocator *) A->data;

 98:   p->nooffproc = PETSC_TRUE;
 99:   while (1) {
100:     MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg);
101:     if (flg) p->nooffproc = PETSC_FALSE;
102:     if (!flg) break;

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

121: PetscErrorCode MatView_Preallocator(Mat A, PetscViewer viewer)
122: {
123:   return 0;
124: }

126: PetscErrorCode MatSetOption_Preallocator(Mat A, MatOption op, PetscBool flg)
127: {
128:   return 0;
129: }

131: PetscErrorCode MatPreallocatorPreallocate_Preallocator(Mat mat, PetscBool fill, Mat A)
132: {
133:   Mat_Preallocator *p = (Mat_Preallocator *) mat->data;
134:   PetscInt          bs;

137:   p->used = PETSC_TRUE;
138:   if (!fill) PetscHSetIJDestroy(&p->ht);
139:   MatGetBlockSize(mat, &bs);
140:   MatXAIJSetPreallocation(A, bs, p->dnz, p->onz, p->dnzu, p->onzu);
141:   MatSetUp(A);
142:   MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
143:   if (fill) {
144:     MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, p->nooffproc);
145:     PetscHashIter  hi;
146:     PetscHashIJKey key;
147:     PetscScalar    *zeros;
148:     PetscInt       n,maxrow=1,*cols,rStart,rEnd,*rowstarts;

150:     MatGetOwnershipRange(A, &rStart, &rEnd);
151:     // Ownership range is in terms of scalar entries, but we deal with blocks
152:     rStart /= bs;
153:     rEnd /= bs;
154:     PetscHSetIJGetSize(p->ht,&n);
155:     PetscMalloc2(n,&cols,rEnd-rStart+1,&rowstarts);
156:     rowstarts[0] = 0;
157:     for (PetscInt i=0; i<rEnd-rStart; i++) {
158:       rowstarts[i+1] = rowstarts[i] + p->dnz[i] + p->onz[i];
159:       maxrow = PetscMax(maxrow, p->dnz[i] + p->onz[i]);
160:     }

163:     PetscHashIterBegin(p->ht,hi);
164:     for (PetscInt i=0; !PetscHashIterAtEnd(p->ht,hi); i++) {
165:       PetscHashIterGetKey(p->ht,hi,key);
166:       PetscInt lrow = key.i - rStart;
167:       cols[rowstarts[lrow]] = key.j;
168:       rowstarts[lrow]++;
169:       PetscHashIterNext(p->ht,hi);
170:     }
171:     PetscHSetIJDestroy(&p->ht);

173:     PetscCalloc1(maxrow*bs*bs,&zeros);
174:     for (PetscInt i=0; i<rEnd-rStart; i++) {
175:       PetscInt grow = rStart + i;
176:       PetscInt end = rowstarts[i], start = end - p->dnz[i] - p->onz[i];
177:       PetscSortInt(end-start,&cols[start]);
178:       MatSetValuesBlocked(A, 1, &grow, end-start, &cols[start], zeros, INSERT_VALUES);
179:     }
180:     PetscFree(zeros);
181:     PetscFree2(cols,rowstarts);

183:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
184:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
185:     MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_FALSE);
186:   }
187:   return 0;
188: }

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

193:   Input Parameters:
194: + mat  - the preallocator
195: . fill - fill the matrix with zeros
196: - A    - the matrix to be preallocated

198:   Notes:
199:   This Mat implementation provides a helper utility to define the correct
200:   preallocation data for a given nonzero structure. Use this object like a
201:   regular matrix, e.g. loop over the nonzero structure of the matrix and
202:   call MatSetValues() or MatSetValuesBlocked() to indicate the nonzero locations.
203:   The matrix entries provided to MatSetValues() will be ignored, it only uses
204:   the row / col indices provided to determine the information required to be
205:   passed to MatXAIJSetPreallocation(). Once you have looped over the nonzero
206:   structure, you must call MatAssemblyBegin(), MatAssemblyEnd() on mat.

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

213:   This function may only be called once for a given MatPreallocator object. If
214:   multiple Mats need to be preallocated, consider using MatDuplicate() after
215:   this function.

217:   Level: advanced

219: .seealso: MATPREALLOCATOR
220: @*/
221: PetscErrorCode MatPreallocatorPreallocate(Mat mat, PetscBool fill, Mat A)
222: {
226:   PetscUseMethod(mat,"MatPreallocatorPreallocate_C",(Mat,PetscBool,Mat),(mat,fill,A));
227:   return 0;
228: }

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

233:    Operations Provided:
234: .vb
235:   MatSetValues()
236: .ve

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

241:   Level: advanced

243: .seealso: Mat, MatPreallocatorPreallocate()

245: M*/

247: PETSC_EXTERN PetscErrorCode MatCreate_Preallocator(Mat A)
248: {
249:   Mat_Preallocator *p;

251:   PetscNewLog(A, &p);
252:   A->data = (void *) p;

254:   p->ht   = NULL;
255:   p->dnz  = NULL;
256:   p->onz  = NULL;
257:   p->dnzu = NULL;
258:   p->onzu = NULL;
259:   p->used = PETSC_FALSE;

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

264:   A->ops->destroy       = MatDestroy_Preallocator;
265:   A->ops->setup         = MatSetUp_Preallocator;
266:   A->ops->setvalues     = MatSetValues_Preallocator;
267:   A->ops->assemblybegin = MatAssemblyBegin_Preallocator;
268:   A->ops->assemblyend   = MatAssemblyEnd_Preallocator;
269:   A->ops->view          = MatView_Preallocator;
270:   A->ops->setoption     = MatSetOption_Preallocator;
271:   A->ops->setblocksizes = MatSetBlockSizes_Default; /* once set, user is not allowed to change the block sizes */

273:   /* special MATPREALLOCATOR functions */
274:   PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", MatPreallocatorPreallocate_Preallocator);
275:   PetscObjectChangeTypeName((PetscObject) A, MATPREALLOCATOR);
276:   return 0;
277: }