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, (PetscErrorCodeFn **)&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;
206: PetscBool option;
208: PetscCall(MatGetSize(X, &m, &n));
209: PetscCall(MatGetOwnershipRange(X, &start, &end));
210: PetscCall(MatGetRowUpperTriangular(X));
211: if (a == 1.0) {
212: for (i = start; i < end; i++) {
213: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
214: PetscCall(MatSetValues(Y, 1, &i, ncols, row, vals, ADD_VALUES));
215: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
216: }
217: } else {
218: PetscInt vs = 100;
219: /* realloc if needed, as this function may be used in parallel */
220: PetscCall(PetscMalloc1(vs, &val));
221: for (i = start; i < end; i++) {
222: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
223: if (vs < ncols) {
224: vs = PetscMin(2 * ncols, n);
225: PetscCall(PetscRealloc(vs * sizeof(*val), &val));
226: }
227: for (j = 0; j < ncols; j++) val[j] = a * vals[j];
228: PetscCall(MatSetValues(Y, 1, &i, ncols, row, val, ADD_VALUES));
229: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
230: }
231: PetscCall(PetscFree(val));
232: }
233: PetscCall(MatRestoreRowUpperTriangular(X));
234: PetscCall(MatGetOption(Y, MAT_NO_OFF_PROC_ENTRIES, &option));
235: PetscCall(MatSetOption(Y, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
236: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
237: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
238: PetscCall(MatSetOption(Y, MAT_NO_OFF_PROC_ENTRIES, option));
239: } else {
240: Mat B;
242: PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B));
243: PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
244: PetscCall(MatHeaderMerge(Y, &B));
245: }
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B, Mat Y, PetscScalar a, Mat X, MatStructure str)
250: {
251: PetscInt i, start, end, j, ncols, m, n;
252: const PetscInt *row;
253: PetscScalar *val;
254: const PetscScalar *vals;
255: PetscBool option;
257: PetscFunctionBegin;
258: PetscCall(MatGetSize(X, &m, &n));
259: PetscCall(MatGetOwnershipRange(X, &start, &end));
260: PetscCall(MatGetRowUpperTriangular(Y));
261: PetscCall(MatGetRowUpperTriangular(X));
262: if (a == 1.0) {
263: for (i = start; i < end; i++) {
264: PetscCall(MatGetRow(Y, i, &ncols, &row, &vals));
265: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
266: PetscCall(MatRestoreRow(Y, i, &ncols, &row, &vals));
268: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
269: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
270: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
271: }
272: } else {
273: PetscInt vs = 100;
274: /* realloc if needed, as this function may be used in parallel */
275: PetscCall(PetscMalloc1(vs, &val));
276: for (i = start; i < end; i++) {
277: PetscCall(MatGetRow(Y, i, &ncols, &row, &vals));
278: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
279: PetscCall(MatRestoreRow(Y, i, &ncols, &row, &vals));
281: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
282: if (vs < ncols) {
283: vs = PetscMin(2 * ncols, n);
284: PetscCall(PetscRealloc(vs * sizeof(*val), &val));
285: }
286: for (j = 0; j < ncols; j++) val[j] = a * vals[j];
287: PetscCall(MatSetValues(B, 1, &i, ncols, row, val, ADD_VALUES));
288: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
289: }
290: PetscCall(PetscFree(val));
291: }
292: PetscCall(MatRestoreRowUpperTriangular(Y));
293: PetscCall(MatRestoreRowUpperTriangular(X));
294: PetscCall(MatGetOption(B, MAT_NO_OFF_PROC_ENTRIES, &option));
295: PetscCall(MatSetOption(B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
296: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
297: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
298: PetscCall(MatSetOption(B, MAT_NO_OFF_PROC_ENTRIES, option));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: /*@
303: MatShift - Computes `Y = Y + a I`, where `a` is a `PetscScalar`
305: Neighbor-wise Collective
307: Input Parameters:
308: + Y - the matrix
309: - a - the `PetscScalar`
311: Level: intermediate
313: Notes:
314: If `Y` is a rectangular matrix, the shift is done on the main diagonal of the matrix (https://en.wikipedia.org/wiki/Main_diagonal)
316: If the matrix `Y` is missing some diagonal entries this routine can be very slow. To make it fast one should initially
317: fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
318: entry). No operation is performed when a is zero.
320: To form Y = Y + diag(V) use `MatDiagonalSet()`
322: .seealso: [](ch_matrices), `Mat`, `MatDiagonalSet()`, `MatScale()`, `MatDiagonalScale()`
323: @*/
324: PetscErrorCode MatShift(Mat Y, PetscScalar a)
325: {
326: PetscFunctionBegin;
328: PetscCheck(Y->assembled, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
329: PetscCheck(!Y->factortype, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
330: MatCheckPreallocated(Y, 1);
331: if (a == 0.0) PetscFunctionReturn(PETSC_SUCCESS);
333: if (Y->ops->shift) PetscUseTypeMethod(Y, shift, a);
334: else PetscCall(MatShift_Basic(Y, a));
336: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: PetscErrorCode MatDiagonalSet_Default(Mat Y, Vec D, InsertMode is)
341: {
342: PetscInt i, start, end;
343: const PetscScalar *v;
345: PetscFunctionBegin;
346: PetscCall(MatGetOwnershipRange(Y, &start, &end));
347: PetscCall(VecGetArrayRead(D, &v));
348: for (i = start; i < end; i++) PetscCall(MatSetValues(Y, 1, &i, 1, &i, v + i - start, is));
349: PetscCall(VecRestoreArrayRead(D, &v));
350: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
351: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@
356: MatDiagonalSet - Computes `Y` = `Y` + `D`, where `D` is a diagonal matrix
357: that is represented as a vector. Or Y[i,i] = D[i] if `InsertMode` is
358: `INSERT_VALUES`.
360: Neighbor-wise Collective
362: Input Parameters:
363: + Y - the input matrix
364: . D - the diagonal matrix, represented as a vector
365: - is - `INSERT_VALUES` or `ADD_VALUES`
367: Level: intermediate
369: Note:
370: If the matrix `Y` is missing some diagonal entries this routine can be very slow. To make it fast one should initially
371: fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
372: entry).
374: .seealso: [](ch_matrices), `Mat`, `MatShift()`, `MatScale()`, `MatDiagonalScale()`
375: @*/
376: PetscErrorCode MatDiagonalSet(Mat Y, Vec D, InsertMode is)
377: {
378: PetscInt matlocal, veclocal;
380: PetscFunctionBegin;
383: MatCheckPreallocated(Y, 1);
384: PetscCall(MatGetLocalSize(Y, &matlocal, NULL));
385: PetscCall(VecGetLocalSize(D, &veclocal));
386: 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);
387: if (Y->ops->diagonalset) PetscUseTypeMethod(Y, diagonalset, D, is);
388: else PetscCall(MatDiagonalSet_Default(Y, D, is));
389: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: /*@
394: MatAYPX - Computes Y = a*Y + X.
396: Logically Collective
398: Input Parameters:
399: + a - the `PetscScalar` multiplier
400: . Y - the first matrix
401: . X - the second matrix
402: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN`, `UNKNOWN_NONZERO_PATTERN`, or `SUBSET_NONZERO_PATTERN` (nonzeros of `X` is a subset of `Y`'s)
404: Level: intermediate
406: .seealso: [](ch_matrices), `Mat`, `MatAXPY()`
407: @*/
408: PetscErrorCode MatAYPX(Mat Y, PetscScalar a, Mat X, MatStructure str)
409: {
410: PetscFunctionBegin;
411: PetscCall(MatScale(Y, a));
412: PetscCall(MatAXPY(Y, 1.0, X, str));
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: /*@
417: MatComputeOperator - Computes the explicit matrix
419: Collective
421: Input Parameters:
422: + inmat - the matrix
423: - mattype - the matrix type for the explicit operator
425: Output Parameter:
426: . mat - the explicit operator
428: Level: advanced
430: Note:
431: This computation is done by applying the operator to columns of the identity matrix.
432: This routine is costly in general, and is recommended for use only with relatively small systems.
433: Currently, this routine uses a dense matrix format if `mattype` == `NULL`.
435: .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperatorTranspose()`
436: @*/
437: PetscErrorCode MatComputeOperator(Mat inmat, MatType mattype, Mat *mat)
438: {
439: PetscFunctionBegin;
441: PetscAssertPointer(mat, 3);
442: PetscCall(MatConvert_Shell(inmat, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /*@
447: MatComputeOperatorTranspose - Computes the explicit matrix representation of
448: a give matrix that can apply `MatMultTranspose()`
450: Collective
452: Input Parameters:
453: + inmat - the matrix
454: - mattype - the matrix type for the explicit operator
456: Output Parameter:
457: . mat - the explicit operator transposed
459: Level: advanced
461: Note:
462: This computation is done by applying the transpose of the operator to columns of the identity matrix.
463: This routine is costly in general, and is recommended for use only with relatively small systems.
464: Currently, this routine uses a dense matrix format if `mattype` == `NULL`.
466: .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperator()`
467: @*/
468: PetscErrorCode MatComputeOperatorTranspose(Mat inmat, MatType mattype, Mat *mat)
469: {
470: Mat A;
472: PetscFunctionBegin;
474: PetscAssertPointer(mat, 3);
475: PetscCall(MatCreateTranspose(inmat, &A));
476: PetscCall(MatConvert_Shell(A, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
477: PetscCall(MatDestroy(&A));
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@
482: 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
484: Input Parameters:
485: + A - The matrix
486: . tol - The zero tolerance
487: . compress - Whether the storage from the input matrix `A` should be compressed once values less than or equal to `tol` are set to zero
488: - 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
490: Level: intermediate
492: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatZeroEntries()`, `MatEliminateZeros()`, `VecFilter()`
493: @*/
494: PetscErrorCode MatFilter(Mat A, PetscReal tol, PetscBool compress, PetscBool keep)
495: {
496: Mat a;
497: PetscScalar *newVals;
498: PetscInt *newCols, rStart, rEnd, maxRows, r, colMax = 0, nnz0 = 0, nnz1 = 0;
499: PetscBool flg;
501: PetscFunctionBegin;
502: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)A, &flg, MATSEQDENSE, MATMPIDENSE, ""));
503: if (flg) {
504: PetscCall(MatDenseGetLocalMatrix(A, &a));
505: PetscCall(MatDenseGetLDA(a, &r));
506: PetscCall(MatGetSize(a, &rStart, &rEnd));
507: PetscCall(MatDenseGetArray(a, &newVals));
508: for (; colMax < rEnd; ++colMax) {
509: for (maxRows = 0; maxRows < rStart; ++maxRows) newVals[maxRows + colMax * r] = PetscAbsScalar(newVals[maxRows + colMax * r]) <= tol ? 0.0 : newVals[maxRows + colMax * r];
510: }
511: PetscCall(MatDenseRestoreArray(a, &newVals));
512: } else {
513: const PetscInt *ranges;
514: PetscMPIInt rank, size;
516: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
517: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
518: PetscCall(MatGetOwnershipRanges(A, &ranges));
519: rStart = ranges[rank];
520: rEnd = ranges[rank + 1];
521: PetscCall(MatGetRowUpperTriangular(A));
522: for (r = rStart; r < rEnd; ++r) {
523: PetscInt ncols;
525: PetscCall(MatGetRow(A, r, &ncols, NULL, NULL));
526: colMax = PetscMax(colMax, ncols);
527: PetscCall(MatRestoreRow(A, r, &ncols, NULL, NULL));
528: }
529: maxRows = 0;
530: for (r = 0; r < size; ++r) maxRows = PetscMax(maxRows, ranges[r + 1] - ranges[r]);
531: PetscCall(PetscCalloc2(colMax, &newCols, colMax, &newVals));
532: PetscCall(MatGetOption(A, MAT_NO_OFF_PROC_ENTRIES, &flg)); /* cache user-defined value */
533: PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
534: /* short-circuit code in MatAssemblyBegin() and MatAssemblyEnd() */
535: /* that are potentially called many times depending on the distribution of A */
536: for (r = rStart; r < rStart + maxRows; ++r) {
537: if (r < rEnd) {
538: const PetscScalar *vals;
539: const PetscInt *cols;
540: PetscInt ncols, newcols = 0, c;
542: PetscCall(MatGetRow(A, r, &ncols, &cols, &vals));
543: nnz0 += ncols - 1;
544: for (c = 0; c < ncols; ++c) {
545: if (PetscUnlikely(PetscAbsScalar(vals[c]) <= tol)) newCols[newcols++] = cols[c];
546: }
547: nnz1 += ncols - newcols - 1;
548: PetscCall(MatRestoreRow(A, r, &ncols, &cols, &vals));
549: PetscCall(MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES));
550: }
551: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
552: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
553: }
554: PetscCall(MatRestoreRowUpperTriangular(A));
555: PetscCall(PetscFree2(newCols, newVals));
556: PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, flg)); /* reset option to its user-defined value */
557: if (nnz0 > 0) PetscCall(PetscInfo(NULL, "Filtering left %g%% edges in graph\n", 100 * (double)nnz1 / (double)nnz0));
558: else PetscCall(PetscInfo(NULL, "Warning: %" PetscInt_FMT " edges to filter with %" PetscInt_FMT " rows\n", nnz0, maxRows));
559: }
560: if (compress && A->ops->eliminatezeros) {
561: Mat B;
562: PetscBool flg;
564: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
565: if (!flg) {
566: PetscCall(MatEliminateZeros(A, keep));
567: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
568: PetscCall(MatHeaderReplace(A, &B));
569: }
570: }
571: PetscFunctionReturn(PETSC_SUCCESS);
572: }