Actual source code: matpreallocator.c
petsc-3.11.4 2019-09-28
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: /* 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: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
147: if (fill) {
148: PetscHashIter hi;
149: PetscHashIJKey key;
150: PetscScalar *zeros;
152: PetscCalloc1(bs*bs,&zeros);
154: PetscHashIterBegin(p->ht,hi);
155: while (!PetscHashIterAtEnd(p->ht,hi)) {
156: PetscHashIterGetKey(p->ht,hi,key);
157: PetscHashIterNext(p->ht,hi);
158: MatSetValuesBlocked(A,1,&key.i,1,&key.j,zeros,INSERT_VALUES);
159: }
160: PetscFree(zeros);
162: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
163: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
164: }
165: return(0);
166: }
168: /*@
169: MatPreallocatorPreallocate - Preallocates the A matrix, using information from mat, optionally filling A with zeros
171: Input Parameter:
172: + mat - the preallocator
173: . fill - fill the matrix with zeros
174: - A - the matrix to be preallocated
176: Notes:
177: This Mat implementation provides a helper utility to define the correct
178: preallocation data for a given nonzero structure. Use this object like a
179: regular matrix, e.g. loop over the nonzero structure of the matrix and
180: call MatSetValues() or MatSetValuesBlocked() to indicate the nonzero locations.
181: The matrix entires provided to MatSetValues() will be ignored, it only uses
182: the row / col indices provided to determine the information required to be
183: passed to MatXAIJSetPreallocation(). Once you have looped over the nonzero
184: structure, you must call MatAssemblyBegin(), MatAssemblyEnd() on mat.
186: After you have assembled the preallocator matrix (mat), call MatPreallocatorPreallocate()
187: to define the preallocation information on the matrix (A). Setting the parameter
188: fill = PETSC_TRUE will insert zeros into the matrix A. Internally MatPreallocatorPreallocate()
189: will call MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
191: Level: advanced
193: .seealso: MATPREALLOCATOR
194: @*/
195: PetscErrorCode MatPreallocatorPreallocate(Mat mat, PetscBool fill, Mat A)
196: {
202: PetscUseMethod(mat, "MatPreallocatorPreallocate_C", (Mat,PetscBool,Mat),(mat,fill,A));
203: return(0);
204: }
206: /*MC
207: MATPREALLOCATOR - MATPREALLOCATOR = "preallocator" - A matrix type to be used for computing a matrix preallocation.
209: Operations Provided:
210: . MatSetValues()
212: Options Database Keys:
213: . -mat_type preallocator - sets the matrix type to "preallocator" during a call to MatSetFromOptions()
215: Level: advanced
217: .seealso: Mat
219: M*/
221: PETSC_EXTERN PetscErrorCode MatCreate_Preallocator(Mat A)
222: {
223: Mat_Preallocator *p;
224: PetscErrorCode ierr;
227: PetscNewLog(A, &p);
228: A->data = (void *) p;
230: p->ht = NULL;
231: p->dnz = NULL;
232: p->onz = NULL;
233: p->dnzu = NULL;
234: p->onzu = NULL;
236: /* matrix ops */
237: PetscMemzero(A->ops, sizeof(struct _MatOps));
239: A->ops->destroy = MatDestroy_Preallocator;
240: A->ops->setup = MatSetUp_Preallocator;
241: A->ops->setvalues = MatSetValues_Preallocator;
242: A->ops->assemblybegin = MatAssemblyBegin_Preallocator;
243: A->ops->assemblyend = MatAssemblyEnd_Preallocator;
244: A->ops->view = MatView_Preallocator;
245: A->ops->setoption = MatSetOption_Preallocator;
246: A->ops->setblocksizes = MatSetBlockSizes_Default; /* once set, user is not allowed to change the block sizes */
248: /* special MATPREALLOCATOR functions */
249: PetscObjectComposeFunction((PetscObject) A, "MatPreallocatorPreallocate_C", MatPreallocatorPreallocate_Preallocator);
250: PetscObjectChangeTypeName((PetscObject) A, MATPREALLOCATOR);
251: return(0);
252: }