Actual source code: matpreallocator.c
petsc-3.14.6 2021-03-30
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: }