Actual source code: axpy.c
1: #include <petsc/private/matimpl.h>
3: static PetscErrorCode MatTransposeAXPY_Private(Mat Y, PetscScalar a, Mat X, MatStructure str, Mat T)
4: {
5: Mat A, F;
6: PetscScalar vshift, vscale;
7: PetscErrorCode (*f)(Mat, Mat *);
9: PetscFunctionBegin;
10: if (T == X) PetscCall(MatShellGetScalingShifts(T, &vshift, &vscale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
11: else {
12: vshift = 0.0;
13: vscale = 1.0;
14: }
15: PetscCall(PetscObjectQueryFunction((PetscObject)T, "MatTransposeGetMat_C", &f));
16: if (f) {
17: PetscCall(MatTransposeGetMat(T, &A));
18: if (T == X) {
19: PetscCall(PetscInfo(NULL, "Explicitly transposing X of type MATTRANSPOSEVIRTUAL to perform MatAXPY()\n"));
20: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &F));
21: A = Y;
22: } else {
23: PetscCall(PetscInfo(NULL, "Transposing X because Y of type MATTRANSPOSEVIRTUAL to perform MatAXPY()\n"));
24: PetscCall(MatTranspose(X, MAT_INITIAL_MATRIX, &F));
25: }
26: } else {
27: PetscCall(MatHermitianTransposeGetMat(T, &A));
28: if (T == X) {
29: PetscCall(PetscInfo(NULL, "Explicitly Hermitian transposing X of type MATHERMITIANTRANSPOSEVIRTUAL to perform MatAXPY()\n"));
30: PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &F));
31: A = Y;
32: } else {
33: PetscCall(PetscInfo(NULL, "Hermitian transposing X because Y of type MATHERMITIANTRANSPOSEVIRTUAL to perform MatAXPY()\n"));
34: PetscCall(MatHermitianTranspose(X, MAT_INITIAL_MATRIX, &F));
35: }
36: }
37: PetscCall(MatAXPY(A, a * vscale, F, str));
38: PetscCall(MatShift(A, a * vshift));
39: PetscCall(MatDestroy(&F));
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode MatAXPY_BasicWithTypeCompare(Mat Y, PetscScalar a, Mat X, MatStructure str)
44: {
45: PetscBool flg;
47: PetscFunctionBegin;
48: PetscCall(MatIsShell(Y, &flg));
49: if (flg) { /* MatShell has special support for AXPY */
50: PetscErrorCode (*f)(Mat, PetscScalar, Mat, MatStructure);
52: PetscCall(MatGetOperation(Y, MATOP_AXPY, (void (**)(void))&f));
53: if (f) {
54: PetscCall((*f)(Y, a, X, str));
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
57: } else {
58: /* no need to preallocate if Y is dense */
59: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &flg, MATSEQDENSE, MATMPIDENSE, ""));
60: if (flg) {
61: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATNEST, &flg));
62: if (flg) {
63: PetscCall(MatAXPY_Dense_Nest(Y, a, X));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
66: }
67: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &flg, MATSCALAPACK, MATELEMENTAL, ""));
68: if (flg) { /* Avoid MatAXPY_Basic() due to missing MatGetRow() */
69: Mat C;
71: PetscCall(MatConvert(X, ((PetscObject)Y)->type_name, MAT_INITIAL_MATRIX, &C));
72: PetscCall(MatAXPY(Y, a, C, str));
73: PetscCall(MatDestroy(&C));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
76: }
77: PetscCall(MatAXPY_Basic(Y, a, X, str));
78: PetscFunctionReturn(PETSC_SUCCESS);
79: }
81: /*@
82: MatAXPY - Computes Y = a*X + Y.
84: Logically Collective
86: Input Parameters:
87: + a - the scalar multiplier
88: . X - the first matrix
89: . Y - the second matrix
90: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN`, `UNKNOWN_NONZERO_PATTERN`, or `SUBSET_NONZERO_PATTERN` (nonzeros of `X` is a subset of `Y`'s)
92: Level: intermediate
94: .seealso: [](ch_matrices), `Mat`, `MatAYPX()`
95: @*/
96: PetscErrorCode MatAXPY(Mat Y, PetscScalar a, Mat X, MatStructure str)
97: {
98: PetscInt M1, M2, N1, N2;
99: PetscInt m1, m2, n1, n2;
100: PetscBool sametype, transpose;
102: PetscFunctionBegin;
106: PetscCheckSameComm(Y, 1, X, 3);
107: PetscCall(MatGetSize(X, &M1, &N1));
108: PetscCall(MatGetSize(Y, &M2, &N2));
109: PetscCall(MatGetLocalSize(X, &m1, &n1));
110: PetscCall(MatGetLocalSize(Y, &m2, &n2));
111: PetscCheck(M1 == M2 && N1 == N2, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_SIZ, "Non conforming matrix add: global sizes X %" PetscInt_FMT " x %" PetscInt_FMT ", Y %" PetscInt_FMT " x %" PetscInt_FMT, M1, N1, M2, N2);
112: PetscCheck(m1 == m2 && n1 == n2, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Non conforming matrix add: local sizes X %" PetscInt_FMT " x %" PetscInt_FMT ", Y %" PetscInt_FMT " x %" PetscInt_FMT, m1, n1, m2, n2);
113: PetscCheck(Y->assembled, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix (Y)");
114: PetscCheck(X->assembled, PetscObjectComm((PetscObject)X), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix (X)");
115: if (a == (PetscScalar)0.0) PetscFunctionReturn(PETSC_SUCCESS);
116: if (Y == X) {
117: PetscCall(MatScale(Y, 1.0 + a));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
120: PetscCall(PetscObjectObjectTypeCompare((PetscObject)X, (PetscObject)Y, &sametype));
121: PetscCall(PetscLogEventBegin(MAT_AXPY, Y, 0, 0, 0));
122: if (Y->ops->axpy && (sametype || X->ops->axpy == Y->ops->axpy)) {
123: PetscUseTypeMethod(Y, axpy, a, X, str);
124: } else {
125: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &transpose, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
126: if (transpose) {
127: PetscCall(MatTransposeAXPY_Private(Y, a, X, str, X));
128: } else {
129: PetscCall(PetscObjectTypeCompareAny((PetscObject)Y, &transpose, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
130: if (transpose) {
131: PetscCall(MatTransposeAXPY_Private(Y, a, X, str, Y));
132: } else {
133: PetscCall(MatAXPY_BasicWithTypeCompare(Y, a, X, str));
134: }
135: }
136: }
137: PetscCall(PetscLogEventEnd(MAT_AXPY, Y, 0, 0, 0));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: PetscErrorCode MatAXPY_Basic_Preallocate(Mat Y, Mat X, Mat *B)
142: {
143: PetscErrorCode (*preall)(Mat, Mat, Mat *) = NULL;
145: PetscFunctionBegin;
146: /* look for any available faster alternative to the general preallocator */
147: PetscCall(PetscObjectQueryFunction((PetscObject)Y, "MatAXPYGetPreallocation_C", &preall));
148: if (preall) {
149: PetscCall((*preall)(Y, X, B));
150: } else { /* Use MatPrellocator, assumes same row-col distribution */
151: Mat preallocator;
152: PetscInt r, rstart, rend;
153: PetscInt m, n, M, N;
155: PetscCall(MatGetRowUpperTriangular(Y));
156: PetscCall(MatGetRowUpperTriangular(X));
157: PetscCall(MatGetSize(Y, &M, &N));
158: PetscCall(MatGetLocalSize(Y, &m, &n));
159: PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &preallocator));
160: PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
161: PetscCall(MatSetLayouts(preallocator, Y->rmap, Y->cmap));
162: PetscCall(MatSetUp(preallocator));
163: PetscCall(MatGetOwnershipRange(preallocator, &rstart, &rend));
164: for (r = rstart; r < rend; ++r) {
165: PetscInt ncols;
166: const PetscInt *row;
167: const PetscScalar *vals;
169: PetscCall(MatGetRow(Y, r, &ncols, &row, &vals));
170: PetscCall(MatSetValues(preallocator, 1, &r, ncols, row, vals, INSERT_VALUES));
171: PetscCall(MatRestoreRow(Y, r, &ncols, &row, &vals));
172: PetscCall(MatGetRow(X, r, &ncols, &row, &vals));
173: PetscCall(MatSetValues(preallocator, 1, &r, ncols, row, vals, INSERT_VALUES));
174: PetscCall(MatRestoreRow(X, r, &ncols, &row, &vals));
175: }
176: PetscCall(MatSetOption(preallocator, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
177: PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
178: PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
179: PetscCall(MatRestoreRowUpperTriangular(Y));
180: PetscCall(MatRestoreRowUpperTriangular(X));
182: PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), B));
183: PetscCall(MatSetType(*B, ((PetscObject)Y)->type_name));
184: PetscCall(MatSetLayouts(*B, Y->rmap, Y->cmap));
185: PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_FALSE, *B));
186: PetscCall(MatDestroy(&preallocator));
187: }
188: PetscFunctionReturn(PETSC_SUCCESS);
189: }
191: PetscErrorCode MatAXPY_Basic(Mat Y, PetscScalar a, Mat X, MatStructure str)
192: {
193: PetscFunctionBegin;
194: if (str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN) {
195: PetscBool isdense;
197: /* no need to preallocate if Y is dense */
198: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
199: if (isdense) str = SUBSET_NONZERO_PATTERN;
200: }
201: if (str != DIFFERENT_NONZERO_PATTERN && str != UNKNOWN_NONZERO_PATTERN) {
202: PetscInt i, start, end, j, ncols, m, n;
203: const PetscInt *row;
204: PetscScalar *val;
205: const PetscScalar *vals;
207: PetscCall(MatGetSize(X, &m, &n));
208: PetscCall(MatGetOwnershipRange(X, &start, &end));
209: PetscCall(MatGetRowUpperTriangular(X));
210: if (a == 1.0) {
211: for (i = start; i < end; i++) {
212: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
213: PetscCall(MatSetValues(Y, 1, &i, ncols, row, vals, ADD_VALUES));
214: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
215: }
216: } else {
217: PetscInt vs = 100;
218: /* realloc if needed, as this function may be used in parallel */
219: PetscCall(PetscMalloc1(vs, &val));
220: for (i = start; i < end; i++) {
221: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
222: if (vs < ncols) {
223: vs = PetscMin(2 * ncols, n);
224: PetscCall(PetscRealloc(vs * sizeof(*val), &val));
225: }
226: for (j = 0; j < ncols; j++) val[j] = a * vals[j];
227: PetscCall(MatSetValues(Y, 1, &i, ncols, row, val, ADD_VALUES));
228: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
229: }
230: PetscCall(PetscFree(val));
231: }
232: PetscCall(MatRestoreRowUpperTriangular(X));
233: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
234: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
235: } else {
236: Mat B;
238: PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B));
239: PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
240: PetscCall(MatHeaderMerge(Y, &B));
241: }
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B, Mat Y, PetscScalar a, Mat X, MatStructure str)
246: {
247: PetscInt i, start, end, j, ncols, m, n;
248: const PetscInt *row;
249: PetscScalar *val;
250: const PetscScalar *vals;
252: PetscFunctionBegin;
253: PetscCall(MatGetSize(X, &m, &n));
254: PetscCall(MatGetOwnershipRange(X, &start, &end));
255: PetscCall(MatGetRowUpperTriangular(Y));
256: PetscCall(MatGetRowUpperTriangular(X));
257: if (a == 1.0) {
258: for (i = start; i < end; i++) {
259: PetscCall(MatGetRow(Y, i, &ncols, &row, &vals));
260: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
261: PetscCall(MatRestoreRow(Y, i, &ncols, &row, &vals));
263: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
264: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
265: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
266: }
267: } else {
268: PetscInt vs = 100;
269: /* realloc if needed, as this function may be used in parallel */
270: PetscCall(PetscMalloc1(vs, &val));
271: for (i = start; i < end; i++) {
272: PetscCall(MatGetRow(Y, i, &ncols, &row, &vals));
273: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
274: PetscCall(MatRestoreRow(Y, i, &ncols, &row, &vals));
276: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
277: if (vs < ncols) {
278: vs = PetscMin(2 * ncols, n);
279: PetscCall(PetscRealloc(vs * sizeof(*val), &val));
280: }
281: for (j = 0; j < ncols; j++) val[j] = a * vals[j];
282: PetscCall(MatSetValues(B, 1, &i, ncols, row, val, ADD_VALUES));
283: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
284: }
285: PetscCall(PetscFree(val));
286: }
287: PetscCall(MatRestoreRowUpperTriangular(Y));
288: PetscCall(MatRestoreRowUpperTriangular(X));
289: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
290: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /*@
295: MatShift - Computes `Y = Y + a I`, where `a` is a `PetscScalar`
297: Neighbor-wise Collective
299: Input Parameters:
300: + Y - the matrix
301: - a - the `PetscScalar`
303: Level: intermediate
305: Notes:
306: If `Y` is a rectangular matrix, the shift is done on the main diagonal of the matrix (https://en.wikipedia.org/wiki/Main_diagonal)
308: If the matrix `Y` is missing some diagonal entries this routine can be very slow. To make it fast one should initially
309: fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
310: entry). No operation is performed when a is zero.
312: To form Y = Y + diag(V) use `MatDiagonalSet()`
314: .seealso: [](ch_matrices), `Mat`, `MatDiagonalSet()`, `MatScale()`, `MatDiagonalScale()`
315: @*/
316: PetscErrorCode MatShift(Mat Y, PetscScalar a)
317: {
318: PetscFunctionBegin;
320: PetscCheck(Y->assembled, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
321: PetscCheck(!Y->factortype, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
322: MatCheckPreallocated(Y, 1);
323: if (a == 0.0) PetscFunctionReturn(PETSC_SUCCESS);
325: if (Y->ops->shift) PetscUseTypeMethod(Y, shift, a);
326: else PetscCall(MatShift_Basic(Y, a));
328: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: PetscErrorCode MatDiagonalSet_Default(Mat Y, Vec D, InsertMode is)
333: {
334: PetscInt i, start, end;
335: const PetscScalar *v;
337: PetscFunctionBegin;
338: PetscCall(MatGetOwnershipRange(Y, &start, &end));
339: PetscCall(VecGetArrayRead(D, &v));
340: for (i = start; i < end; i++) PetscCall(MatSetValues(Y, 1, &i, 1, &i, v + i - start, is));
341: PetscCall(VecRestoreArrayRead(D, &v));
342: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
343: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: /*@
348: MatDiagonalSet - Computes `Y` = `Y` + `D`, where `D` is a diagonal matrix
349: that is represented as a vector. Or Y[i,i] = D[i] if `InsertMode` is
350: `INSERT_VALUES`.
352: Neighbor-wise Collective
354: Input Parameters:
355: + Y - the input matrix
356: . D - the diagonal matrix, represented as a vector
357: - is - `INSERT_VALUES` or `ADD_VALUES`
359: Level: intermediate
361: Note:
362: If the matrix `Y` is missing some diagonal entries this routine can be very slow. To make it fast one should initially
363: fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
364: entry).
366: .seealso: [](ch_matrices), `Mat`, `MatShift()`, `MatScale()`, `MatDiagonalScale()`
367: @*/
368: PetscErrorCode MatDiagonalSet(Mat Y, Vec D, InsertMode is)
369: {
370: PetscInt matlocal, veclocal;
372: PetscFunctionBegin;
375: MatCheckPreallocated(Y, 1);
376: PetscCall(MatGetLocalSize(Y, &matlocal, NULL));
377: PetscCall(VecGetLocalSize(D, &veclocal));
378: PetscCheck(matlocal == veclocal, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Number local rows of matrix %" PetscInt_FMT " does not match that of vector for diagonal %" PetscInt_FMT, matlocal, veclocal);
379: if (Y->ops->diagonalset) PetscUseTypeMethod(Y, diagonalset, D, is);
380: else PetscCall(MatDiagonalSet_Default(Y, D, is));
381: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: /*@
386: MatAYPX - Computes Y = a*Y + X.
388: Logically Collective
390: Input Parameters:
391: + a - the `PetscScalar` multiplier
392: . Y - the first matrix
393: . X - the second matrix
394: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN`, `UNKNOWN_NONZERO_PATTERN`, or `SUBSET_NONZERO_PATTERN` (nonzeros of `X` is a subset of `Y`'s)
396: Level: intermediate
398: .seealso: [](ch_matrices), `Mat`, `MatAXPY()`
399: @*/
400: PetscErrorCode MatAYPX(Mat Y, PetscScalar a, Mat X, MatStructure str)
401: {
402: PetscFunctionBegin;
403: PetscCall(MatScale(Y, a));
404: PetscCall(MatAXPY(Y, 1.0, X, str));
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@
409: MatComputeOperator - Computes the explicit matrix
411: Collective
413: Input Parameters:
414: + inmat - the matrix
415: - mattype - the matrix type for the explicit operator
417: Output Parameter:
418: . mat - the explicit operator
420: Level: advanced
422: Note:
423: This computation is done by applying the operator to columns of the identity matrix.
424: This routine is costly in general, and is recommended for use only with relatively small systems.
425: Currently, this routine uses a dense matrix format if `mattype` == `NULL`.
427: .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperatorTranspose()`
428: @*/
429: PetscErrorCode MatComputeOperator(Mat inmat, MatType mattype, Mat *mat)
430: {
431: PetscFunctionBegin;
433: PetscAssertPointer(mat, 3);
434: PetscCall(MatConvert_Shell(inmat, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@
439: MatComputeOperatorTranspose - Computes the explicit matrix representation of
440: a give matrix that can apply `MatMultTranspose()`
442: Collective
444: Input Parameters:
445: + inmat - the matrix
446: - mattype - the matrix type for the explicit operator
448: Output Parameter:
449: . mat - the explicit operator transposed
451: Level: advanced
453: Note:
454: This computation is done by applying the transpose of the operator to columns of the identity matrix.
455: This routine is costly in general, and is recommended for use only with relatively small systems.
456: Currently, this routine uses a dense matrix format if `mattype` == `NULL`.
458: .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperator()`
459: @*/
460: PetscErrorCode MatComputeOperatorTranspose(Mat inmat, MatType mattype, Mat *mat)
461: {
462: Mat A;
464: PetscFunctionBegin;
466: PetscAssertPointer(mat, 3);
467: PetscCall(MatCreateTranspose(inmat, &A));
468: PetscCall(MatConvert_Shell(A, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
469: PetscCall(MatDestroy(&A));
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: /*@
474: MatFilter - Set all values in the matrix with an absolute value less than or equal to the tolerance to zero, and optionally compress the underlying storage
476: Input Parameters:
477: + A - The matrix
478: . tol - The zero tolerance
479: . compress - Whether the storage from the input matrix `A` should be compressed once values less than or equal to `tol` are set to zero
480: - keep - If `compress` is true and for a given row of `A`, the diagonal coefficient is less than or equal to `tol`, indicates whether it should be left in the structure or eliminated as well
482: Level: intermediate
484: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatZeroEntries()`, `MatEliminateZeros()`, `VecFilter()`
485: @*/
486: PetscErrorCode MatFilter(Mat A, PetscReal tol, PetscBool compress, PetscBool keep)
487: {
488: Mat a;
489: PetscScalar *newVals;
490: PetscInt *newCols, rStart, rEnd, maxRows, r, colMax = 0, nnz0 = 0, nnz1 = 0;
491: PetscBool flg;
493: PetscFunctionBegin;
494: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)A, &flg, MATSEQDENSE, MATMPIDENSE, ""));
495: if (flg) {
496: PetscCall(MatDenseGetLocalMatrix(A, &a));
497: PetscCall(MatDenseGetLDA(a, &r));
498: PetscCall(MatGetSize(a, &rStart, &rEnd));
499: PetscCall(MatDenseGetArray(a, &newVals));
500: for (; colMax < rEnd; ++colMax) {
501: for (maxRows = 0; maxRows < rStart; ++maxRows) newVals[maxRows + colMax * r] = PetscAbsScalar(newVals[maxRows + colMax * r]) <= tol ? 0.0 : newVals[maxRows + colMax * r];
502: }
503: PetscCall(MatDenseRestoreArray(a, &newVals));
504: } else {
505: const PetscInt *ranges;
506: PetscMPIInt rank, size;
508: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
509: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
510: PetscCall(MatGetOwnershipRanges(A, &ranges));
511: rStart = ranges[rank];
512: rEnd = ranges[rank + 1];
513: PetscCall(MatGetRowUpperTriangular(A));
514: for (r = rStart; r < rEnd; ++r) {
515: PetscInt ncols;
517: PetscCall(MatGetRow(A, r, &ncols, NULL, NULL));
518: colMax = PetscMax(colMax, ncols);
519: PetscCall(MatRestoreRow(A, r, &ncols, NULL, NULL));
520: }
521: maxRows = 0;
522: for (r = 0; r < size; ++r) maxRows = PetscMax(maxRows, ranges[r + 1] - ranges[r]);
523: PetscCall(PetscCalloc2(colMax, &newCols, colMax, &newVals));
524: PetscCall(MatGetOption(A, MAT_NO_OFF_PROC_ENTRIES, &flg)); /* cache user-defined value */
525: PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
526: /* short-circuit code in MatAssemblyBegin() and MatAssemblyEnd() */
527: /* that are potentially called many times depending on the distribution of A */
528: for (r = rStart; r < rStart + maxRows; ++r) {
529: if (r < rEnd) {
530: const PetscScalar *vals;
531: const PetscInt *cols;
532: PetscInt ncols, newcols = 0, c;
534: PetscCall(MatGetRow(A, r, &ncols, &cols, &vals));
535: nnz0 += ncols - 1;
536: for (c = 0; c < ncols; ++c) {
537: if (PetscUnlikely(PetscAbsScalar(vals[c]) <= tol)) newCols[newcols++] = cols[c];
538: }
539: nnz1 += ncols - newcols - 1;
540: PetscCall(MatRestoreRow(A, r, &ncols, &cols, &vals));
541: PetscCall(MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES));
542: }
543: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
544: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
545: }
546: PetscCall(MatRestoreRowUpperTriangular(A));
547: PetscCall(PetscFree2(newCols, newVals));
548: PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, flg)); /* reset option to its user-defined value */
549: if (nnz0 > 0) PetscCall(PetscInfo(NULL, "Filtering left %g%% edges in graph\n", 100 * (double)nnz1 / (double)nnz0));
550: else PetscCall(PetscInfo(NULL, "Warning: %" PetscInt_FMT " edges to filter with %" PetscInt_FMT " rows\n", nnz0, maxRows));
551: }
552: if (compress && A->ops->eliminatezeros) {
553: Mat B;
554: PetscBool flg;
556: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
557: if (!flg) {
558: PetscCall(MatEliminateZeros(A, keep));
559: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
560: PetscCall(MatHeaderReplace(A, &B));
561: }
562: }
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }