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: } Mat_Preallocator;
11: PetscErrorCode MatDestroy_Preallocator(Mat A)
12: {
13: Mat_Preallocator *p = (Mat_Preallocator *) A->data;
14: PetscErrorCode ierr;
17: MatStashDestroy_Private(&A->stash);
18: PetscHSetIJDestroy(&p->ht);
19: PetscFree4(p->dnz, p->onz, p->dnzu, p->onzu);
20: PetscFree(A->data);
21: PetscObjectChangeTypeName((PetscObject) A, NULL);
22: PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", NULL);
23: return(0);
24: }
26: PetscErrorCode MatSetUp_Preallocator(Mat A)
27: {
28: Mat_Preallocator *p = (Mat_Preallocator *) A->data;
29: PetscInt m, bs, mbs;
30: PetscErrorCode ierr;
33: PetscLayoutSetUp(A->rmap);
34: PetscLayoutSetUp(A->cmap);
35: MatGetLocalSize(A, &m, NULL);
36: PetscHSetIJCreate(&p->ht);
37: MatGetBlockSize(A, &bs);
38: /* Do not bother bstash since MatPreallocator does not implement MatSetValuesBlocked */
39: MatStashCreate_Private(PetscObjectComm((PetscObject) A), 1, &A->stash);
40: /* arrays are for blocked rows/cols */
41: mbs = m/bs;
42: PetscCalloc4(mbs, &p->dnz, mbs, &p->onz, mbs, &p->dnzu, mbs, &p->onzu);
43: return(0);
44: }
46: PetscErrorCode MatSetValues_Preallocator(Mat A, PetscInt m, const PetscInt *rows, PetscInt n, const PetscInt *cols, const PetscScalar *values, InsertMode addv)
47: {
48: Mat_Preallocator *p = (Mat_Preallocator *) A->data;
49: PetscInt rStart, rEnd, r, cStart, cEnd, c, bs;
50: PetscErrorCode ierr;
53: MatGetBlockSize(A, &bs);
54: MatGetOwnershipRange(A, &rStart, &rEnd);
55: MatGetOwnershipRangeColumn(A, &cStart, &cEnd);
56: for (r = 0; r < m; ++r) {
57: PetscHashIJKey key;
58: PetscBool missing;
60: key.i = rows[r];
61: if (key.i < 0) continue;
62: if ((key.i < rStart) || (key.i >= rEnd)) {
63: MatStashValuesRow_Private(&A->stash, key.i, n, cols, values, PETSC_FALSE);
64: } else { /* Hash table is for blocked rows/cols */
65: key.i = rows[r]/bs;
66: for (c = 0; c < n; ++c) {
67: key.j = cols[c]/bs;
68: if (key.j < 0) continue;
69: PetscHSetIJQueryAdd(p->ht, key, &missing);
70: if (missing) {
71: if ((key.j >= cStart/bs) && (key.j < cEnd/bs)) {
72: ++p->dnz[key.i-rStart/bs];
73: if (key.j >= key.i) ++p->dnzu[key.i-rStart/bs];
74: } else {
75: ++p->onz[key.i-rStart/bs];
76: if (key.j >= key.i) ++p->onzu[key.i-rStart/bs];
77: }
78: }
79: }
80: }
81: }
82: return(0);
83: }
85: PetscErrorCode MatAssemblyBegin_Preallocator(Mat A, MatAssemblyType type)
86: {
87: PetscInt nstash, reallocs;
91: MatStashScatterBegin_Private(A, &A->stash, A->rmap->range);
92: MatStashGetInfo_Private(&A->stash, &nstash, &reallocs);
93: PetscInfo2(A, "Stash has %D entries, uses %D mallocs.\n", nstash, reallocs);
94: return(0);
95: }
97: PetscErrorCode MatAssemblyEnd_Preallocator(Mat A, MatAssemblyType type)
98: {
99: PetscScalar *val;
100: PetscInt *row, *col;
101: PetscInt i, j, rstart, ncols, flg;
102: PetscMPIInt n;
103: Mat_Preallocator *p = (Mat_Preallocator *) A->data;
104: PetscErrorCode ierr;
107: p->nooffproc = PETSC_TRUE;
108: while (1) {
109: MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg);
110: if (flg) p->nooffproc = PETSC_FALSE;
111: if (!flg) break;
113: for (i = 0; i < n;) {
114: /* Now identify the consecutive vals belonging to the same row */
115: for (j = i, rstart = row[j]; j < n; j++) {
116: if (row[j] != rstart) break;
117: }
118: if (j < n) ncols = j-i;
119: else ncols = n-i;
120: /* Now assemble all these values with a single function call */
121: MatSetValues_Preallocator(A, 1, row+i, ncols, col+i, val+i, INSERT_VALUES);
122: i = j;
123: }
124: }
125: MatStashScatterEnd_Private(&A->stash);
126: MPIU_Allreduce(MPI_IN_PLACE,&p->nooffproc,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)A));
127: return(0);
128: }
130: PetscErrorCode MatView_Preallocator(Mat A, PetscViewer viewer)
131: {
133: return(0);
134: }
136: PetscErrorCode MatSetOption_Preallocator(Mat A, MatOption op, PetscBool flg)
137: {
139: return(0);
140: }
142: PetscErrorCode MatPreallocatorPreallocate_Preallocator(Mat mat, PetscBool fill, Mat A)
143: {
144: Mat_Preallocator *p = (Mat_Preallocator *) mat->data;
145: PetscInt bs;
146: PetscErrorCode ierr;
149: MatGetBlockSize(mat, &bs);
150: MatXAIJSetPreallocation(A, bs, p->dnz, p->onz, p->dnzu, p->onzu);
151: MatSetUp(A);
152: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
153: MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, p->nooffproc);
154: if (fill) {
155: PetscHashIter hi;
156: PetscHashIJKey key;
157: PetscScalar *zeros;
159: PetscCalloc1(bs*bs,&zeros);
161: PetscHashIterBegin(p->ht,hi);
162: while (!PetscHashIterAtEnd(p->ht,hi)) {
163: PetscHashIterGetKey(p->ht,hi,key);
164: PetscHashIterNext(p->ht,hi);
165: MatSetValuesBlocked(A,1,&key.i,1,&key.j,zeros,INSERT_VALUES);
166: }
167: PetscFree(zeros);
169: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
170: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
171: }
172: return(0);
173: }
175: /*@
176: MatPreallocatorPreallocate - Preallocates the A matrix, using information from mat, optionally filling A with zeros
178: Input Parameter:
179: + mat - the preallocator
180: . fill - fill the matrix with zeros
181: - A - the matrix to be preallocated
183: Notes:
184: This Mat implementation provides a helper utility to define the correct
185: preallocation data for a given nonzero structure. Use this object like a
186: regular matrix, e.g. loop over the nonzero structure of the matrix and
187: call MatSetValues() or MatSetValuesBlocked() to indicate the nonzero locations.
188: The matrix entires provided to MatSetValues() will be ignored, it only uses
189: the row / col indices provided to determine the information required to be
190: passed to MatXAIJSetPreallocation(). Once you have looped over the nonzero
191: structure, you must call MatAssemblyBegin(), MatAssemblyEnd() on mat.
193: After you have assembled the preallocator matrix (mat), call MatPreallocatorPreallocate()
194: to define the preallocation information on the matrix (A). Setting the parameter
195: fill = PETSC_TRUE will insert zeros into the matrix A. Internally MatPreallocatorPreallocate()
196: will call MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
198: Level: advanced
200: .seealso: MATPREALLOCATOR
201: @*/
202: PetscErrorCode MatPreallocatorPreallocate(Mat mat, PetscBool fill, Mat A)
203: {
209: PetscUseMethod(mat, "MatPreallocatorPreallocate_C", (Mat,PetscBool,Mat),(mat,fill,A));
210: return(0);
211: }
213: /*MC
214: MATPREALLOCATOR - MATPREALLOCATOR = "preallocator" - A matrix type to be used for computing a matrix preallocation.
216: Operations Provided:
217: . MatSetValues()
219: Options Database Keys:
220: . -mat_type preallocator - sets the matrix type to "preallocator" during a call to MatSetFromOptions()
222: Level: advanced
224: .seealso: Mat
226: M*/
228: PETSC_EXTERN PetscErrorCode MatCreate_Preallocator(Mat A)
229: {
230: Mat_Preallocator *p;
231: PetscErrorCode ierr;
234: PetscNewLog(A, &p);
235: A->data = (void *) p;
237: p->ht = NULL;
238: p->dnz = NULL;
239: p->onz = NULL;
240: p->dnzu = NULL;
241: p->onzu = NULL;
243: /* matrix ops */
244: PetscMemzero(A->ops, sizeof(struct _MatOps));
246: A->ops->destroy = MatDestroy_Preallocator;
247: A->ops->setup = MatSetUp_Preallocator;
248: A->ops->setvalues = MatSetValues_Preallocator;
249: A->ops->assemblybegin = MatAssemblyBegin_Preallocator;
250: A->ops->assemblyend = MatAssemblyEnd_Preallocator;
251: A->ops->view = MatView_Preallocator;
252: A->ops->setoption = MatSetOption_Preallocator;
253: A->ops->setblocksizes = MatSetBlockSizes_Default; /* once set, user is not allowed to change the block sizes */
255: /* special MATPREALLOCATOR functions */
256: PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", MatPreallocatorPreallocate_Preallocator);
257: PetscObjectChangeTypeName((PetscObject) A, MATPREALLOCATOR);
258: return(0);
259: }