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: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATCONSTANTDIAGONAL, &flg));
77: if (flg) {
78: PetscScalar b;
80: PetscCall(MatConstantDiagonalGetConstant(X, &b));
81: PetscCall(MatShift(Y, a * b));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
84: PetscCall(PetscObjectTypeCompare((PetscObject)X, MATDIAGONAL, &flg));
85: if (flg) {
86: Vec d_X;
88: PetscCall(MatDiagonalGetDiagonal(X, &d_X));
89: if (a == 1.0) PetscCall(MatDiagonalSet(Y, d_X, ADD_VALUES));
90: else {
91: Vec d_Y;
93: PetscCall(PetscObjectQuery((PetscObject)Y, "__MatAXPY_BasicWithTypeCompare_Diagonal", (PetscObject *)&d_Y));
94: if (!d_Y) {
95: Vec _d_Y;
96: PetscCall(VecDuplicate(d_X, &_d_Y));
97: PetscCall(PetscObjectCompose((PetscObject)Y, "__MatAXPY_BasicWithTypeCompare_Diagonal", (PetscObject)_d_Y));
98: d_Y = _d_Y;
99: PetscCall(VecDestroy(&_d_Y));
100: }
101: PetscCall(VecAXPBY(d_Y, a, 0.0, d_X));
102: PetscCall(MatDiagonalSet(Y, d_Y, ADD_VALUES));
103: }
104: PetscCall(MatDiagonalRestoreDiagonal(X, &d_X));
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
107: }
108: PetscCall(MatAXPY_Basic(Y, a, X, str));
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: /*@
113: MatAXPY - Computes Y = a*X + Y.
115: Logically Collective
117: Input Parameters:
118: + a - the scalar multiplier
119: . X - the first matrix
120: . Y - the second matrix
121: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN`, `UNKNOWN_NONZERO_PATTERN`, or `SUBSET_NONZERO_PATTERN` (nonzeros of `X` is a subset of `Y`'s)
123: Level: intermediate
125: Note:
126: If `Y` `MAT_STRUCTURE_ONLY` option is set to true, only the union of the nonzero patterns is computed, without accessing numerical values.
127: As with numerical matrices, `a` equal to zero leaves `Y` unchanged.
129: .seealso: [](ch_matrices), `Mat`, `MatAYPX()`
130: @*/
131: PetscErrorCode MatAXPY(Mat Y, PetscScalar a, Mat X, MatStructure str)
132: {
133: PetscInt M1, M2, N1, N2;
134: PetscInt m1, m2, n1, n2;
135: PetscBool sametype, transpose;
137: PetscFunctionBegin;
141: PetscCheckSameComm(Y, 1, X, 3);
142: PetscCall(MatGetSize(X, &M1, &N1));
143: PetscCall(MatGetSize(Y, &M2, &N2));
144: PetscCall(MatGetLocalSize(X, &m1, &n1));
145: PetscCall(MatGetLocalSize(Y, &m2, &n2));
146: 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);
147: 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);
148: PetscCheck(Y->assembled, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix (Y)");
149: PetscCheck(X->assembled, PetscObjectComm((PetscObject)X), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix (X)");
150: if (a == 0.0) PetscFunctionReturn(PETSC_SUCCESS);
151: PetscCheck(!X->structure_only || Y->structure_only, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_INCOMP, "Cannot add a structure-only matrix to a matrix with numerical values");
152: if (Y->structure_only && (Y == X || str == SAME_NONZERO_PATTERN || str == SUBSET_NONZERO_PATTERN)) PetscFunctionReturn(PETSC_SUCCESS);
153: if (Y == X) {
154: PetscCall(MatScale(Y, 1.0 + a));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
157: PetscCall(PetscObjectObjectTypeCompare((PetscObject)X, (PetscObject)Y, &sametype));
158: PetscCall(PetscLogEventBegin(MAT_AXPY, Y, 0, 0, 0));
159: if (Y->structure_only) PetscCall(MatAXPY_Basic(Y, a, X, str));
160: else if (Y->ops->axpy && (sametype || X->ops->axpy == Y->ops->axpy)) PetscUseTypeMethod(Y, axpy, a, X, str);
161: else {
162: PetscCall(PetscObjectTypeCompareAny((PetscObject)X, &transpose, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
163: if (transpose) PetscCall(MatTransposeAXPY_Private(Y, a, X, str, X));
164: else {
165: PetscCall(PetscObjectTypeCompareAny((PetscObject)Y, &transpose, MATTRANSPOSEVIRTUAL, MATHERMITIANTRANSPOSEVIRTUAL, ""));
166: if (transpose) PetscCall(MatTransposeAXPY_Private(Y, a, X, str, Y));
167: else PetscCall(MatAXPY_BasicWithTypeCompare(Y, a, X, str));
168: }
169: }
170: PetscCall(PetscLogEventEnd(MAT_AXPY, Y, 0, 0, 0));
171: PetscFunctionReturn(PETSC_SUCCESS);
172: }
174: PetscErrorCode MatAXPY_Basic_Preallocate(Mat Y, Mat X, Mat *B)
175: {
176: PetscErrorCode (*preall)(Mat, Mat, Mat *) = NULL;
178: PetscFunctionBegin;
179: /* look for any available faster alternative to the general preallocator */
180: PetscCall(PetscObjectQueryFunction((PetscObject)Y, "MatAXPYGetPreallocation_C", &preall));
181: if (preall) PetscCall((*preall)(Y, X, B));
182: else {
183: /* use MatPrellocator, assumes same row-col distribution */
184: Mat preallocator;
185: PetscInt rstart, rend;
186: PetscInt m, n, M, N;
188: PetscCall(MatGetRowUpperTriangular(Y));
189: PetscCall(MatGetRowUpperTriangular(X));
190: PetscCall(MatGetSize(Y, &M, &N));
191: PetscCall(MatGetLocalSize(Y, &m, &n));
192: PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), &preallocator));
193: PetscCall(MatSetType(preallocator, MATPREALLOCATOR));
194: PetscCall(MatSetLayouts(preallocator, Y->rmap, Y->cmap));
195: PetscCall(MatSetUp(preallocator));
196: PetscCall(MatGetOwnershipRange(preallocator, &rstart, &rend));
197: for (PetscInt r = rstart; r < rend; ++r) {
198: PetscInt ncols;
199: const PetscInt *row;
201: PetscCall(MatGetRow(Y, r, &ncols, &row, NULL));
202: PetscCall(MatSetValues(preallocator, 1, &r, ncols, row, NULL, INSERT_VALUES));
203: PetscCall(MatRestoreRow(Y, r, &ncols, &row, NULL));
204: PetscCall(MatGetRow(X, r, &ncols, &row, NULL));
205: PetscCall(MatSetValues(preallocator, 1, &r, ncols, row, NULL, INSERT_VALUES));
206: PetscCall(MatRestoreRow(X, r, &ncols, &row, NULL));
207: }
208: PetscCall(MatSetOption(preallocator, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
209: PetscCall(MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY));
210: PetscCall(MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY));
211: PetscCall(MatRestoreRowUpperTriangular(Y));
212: PetscCall(MatRestoreRowUpperTriangular(X));
214: PetscCall(MatCreate(PetscObjectComm((PetscObject)Y), B));
215: PetscCall(MatSetType(*B, ((PetscObject)Y)->type_name));
216: PetscCall(MatSetLayouts(*B, Y->rmap, Y->cmap));
217: PetscCall(MatSetOption(*B, MAT_STRUCTURE_ONLY, Y->structure_only));
218: PetscCall(MatPreallocatorPreallocate(preallocator, PETSC_FALSE, *B));
219: PetscCall(MatDestroy(&preallocator));
220: }
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: PetscErrorCode MatAXPY_Basic(Mat Y, PetscScalar a, Mat X, MatStructure str)
225: {
226: PetscFunctionBegin;
227: if (str == DIFFERENT_NONZERO_PATTERN || str == UNKNOWN_NONZERO_PATTERN) {
228: PetscBool isdense;
230: /* no need to preallocate if Y is dense */
231: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)Y, &isdense, MATSEQDENSE, MATMPIDENSE, ""));
232: if (isdense) str = SUBSET_NONZERO_PATTERN;
233: }
234: if (str != DIFFERENT_NONZERO_PATTERN && str != UNKNOWN_NONZERO_PATTERN) {
235: PetscInt start, end, ncols, m, n;
236: const PetscInt *row;
237: PetscScalar *val;
238: const PetscScalar *vals;
239: PetscBool option;
241: PetscCall(MatGetSize(X, &m, &n));
242: PetscCall(MatGetOwnershipRange(X, &start, &end));
243: PetscCall(MatGetRowUpperTriangular(X));
244: if (a == 1.0) {
245: for (PetscInt i = start; i < end; i++) {
246: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
247: PetscCall(MatSetValues(Y, 1, &i, ncols, row, vals, ADD_VALUES));
248: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
249: }
250: } else {
251: PetscInt vs = 100;
252: /* realloc if needed, as this function may be used in parallel */
253: PetscCall(PetscMalloc1(vs, &val));
254: for (PetscInt i = start; i < end; i++) {
255: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
256: if (vs < ncols) {
257: vs = PetscMin(2 * ncols, n);
258: PetscCall(PetscRealloc(vs * sizeof(*val), &val));
259: }
260: for (PetscInt j = 0; j < ncols; j++) val[j] = a * vals[j];
261: PetscCall(MatSetValues(Y, 1, &i, ncols, row, val, ADD_VALUES));
262: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
263: }
264: PetscCall(PetscFree(val));
265: }
266: PetscCall(MatRestoreRowUpperTriangular(X));
267: PetscCall(MatGetOption(Y, MAT_NO_OFF_PROC_ENTRIES, &option));
268: PetscCall(MatSetOption(Y, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
269: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
270: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
271: PetscCall(MatSetOption(Y, MAT_NO_OFF_PROC_ENTRIES, option));
272: } else {
273: Mat B;
275: PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B));
276: PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
277: PetscCall(MatHeaderMerge(Y, &B));
278: }
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: PetscErrorCode MatAXPY_BasicWithPreallocation(Mat B, Mat Y, PetscScalar a, Mat X, MatStructure str)
283: {
284: PetscInt start, end, ncols, m, n;
285: const PetscInt *row;
286: PetscScalar *val;
287: const PetscScalar *vals;
288: PetscBool option;
290: PetscFunctionBegin;
291: PetscCall(MatGetSize(X, &m, &n));
292: PetscCall(MatGetOwnershipRange(X, &start, &end));
293: PetscCall(MatGetRowUpperTriangular(Y));
294: PetscCall(MatGetRowUpperTriangular(X));
295: if (Y->structure_only) {
296: PetscCall(MatSetOption(B, MAT_STRUCTURE_ONLY, PETSC_TRUE));
297: for (PetscInt i = start; i < end; i++) {
298: PetscCall(MatGetRow(Y, i, &ncols, &row, NULL));
299: PetscCall(MatSetValues(B, 1, &i, ncols, row, NULL, INSERT_VALUES));
300: PetscCall(MatRestoreRow(Y, i, &ncols, &row, NULL));
302: PetscCall(MatGetRow(X, i, &ncols, &row, NULL));
303: PetscCall(MatSetValues(B, 1, &i, ncols, row, NULL, INSERT_VALUES));
304: PetscCall(MatRestoreRow(X, i, &ncols, &row, NULL));
305: }
306: } else if (a == 1.0) {
307: for (PetscInt i = start; i < end; i++) {
308: PetscCall(MatGetRow(Y, i, &ncols, &row, &vals));
309: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
310: PetscCall(MatRestoreRow(Y, i, &ncols, &row, &vals));
312: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
313: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
314: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
315: }
316: } else {
317: PetscInt vs = 100;
319: /* realloc if needed, as this function may be used in parallel */
320: PetscCall(PetscMalloc1(vs, &val));
321: for (PetscInt i = start; i < end; i++) {
322: PetscCall(MatGetRow(Y, i, &ncols, &row, &vals));
323: PetscCall(MatSetValues(B, 1, &i, ncols, row, vals, ADD_VALUES));
324: PetscCall(MatRestoreRow(Y, i, &ncols, &row, &vals));
326: PetscCall(MatGetRow(X, i, &ncols, &row, &vals));
327: if (vs < ncols) {
328: vs = PetscMin(2 * ncols, n);
329: PetscCall(PetscRealloc(vs * sizeof(*val), &val));
330: }
331: for (PetscInt j = 0; j < ncols; j++) val[j] = a * vals[j];
332: PetscCall(MatSetValues(B, 1, &i, ncols, row, val, ADD_VALUES));
333: PetscCall(MatRestoreRow(X, i, &ncols, &row, &vals));
334: }
335: PetscCall(PetscFree(val));
336: }
337: PetscCall(MatRestoreRowUpperTriangular(Y));
338: PetscCall(MatRestoreRowUpperTriangular(X));
339: PetscCall(MatGetOption(B, MAT_NO_OFF_PROC_ENTRIES, &option));
340: PetscCall(MatSetOption(B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
341: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
342: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
343: PetscCall(MatSetOption(B, MAT_NO_OFF_PROC_ENTRIES, option));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: /*@
348: MatShift - Computes `Y = Y + a I`, where `a` is a `PetscScalar`
350: Neighbor-wise Collective
352: Input Parameters:
353: + Y - the matrix
354: - a - the `PetscScalar`
356: Level: intermediate
358: Notes:
359: If `Y` is a rectangular matrix, the shift is done on the main diagonal of the matrix (https://en.wikipedia.org/wiki/Main_diagonal)
361: If the matrix `Y` is missing some diagonal entries this routine can be very slow. To make it fast one should initially
362: fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
363: entry). No operation is performed when a is zero.
365: To form Y = Y + diag(V) use `MatDiagonalSet()`
367: .seealso: [](ch_matrices), `Mat`, `MatDiagonalSet()`, `MatScale()`, `MatDiagonalScale()`
368: @*/
369: PetscErrorCode MatShift(Mat Y, PetscScalar a)
370: {
371: PetscFunctionBegin;
373: PetscCheck(Y->assembled, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for unassembled matrix");
374: PetscCheck(!Y->factortype, PetscObjectComm((PetscObject)Y), PETSC_ERR_ARG_WRONGSTATE, "Not for factored matrix");
375: MatCheckPreallocated(Y, 1);
376: if (a == 0.0) PetscFunctionReturn(PETSC_SUCCESS);
378: if (Y->ops->shift) PetscUseTypeMethod(Y, shift, a);
379: else PetscCall(MatShift_Basic(Y, a));
381: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: PetscErrorCode MatDiagonalSet_Default(Mat Y, Vec D, InsertMode is)
386: {
387: PetscInt i, start, end;
388: const PetscScalar *v;
390: PetscFunctionBegin;
391: PetscCall(MatGetOwnershipRange(Y, &start, &end));
392: PetscCall(VecGetArrayRead(D, &v));
393: for (i = start; i < end; i++) PetscCall(MatSetValues(Y, 1, &i, 1, &i, v + i - start, is));
394: PetscCall(VecRestoreArrayRead(D, &v));
395: PetscCall(MatAssemblyBegin(Y, MAT_FINAL_ASSEMBLY));
396: PetscCall(MatAssemblyEnd(Y, MAT_FINAL_ASSEMBLY));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: /*@
401: MatDiagonalSet - Computes `Y` = `Y` + `D`, where `D` is a diagonal matrix
402: that is represented as a vector. Or Y[i,i] = D[i] if `InsertMode` is
403: `INSERT_VALUES`.
405: Neighbor-wise Collective
407: Input Parameters:
408: + Y - the input matrix
409: . D - the diagonal matrix, represented as a vector
410: - is - `INSERT_VALUES` or `ADD_VALUES`
412: Level: intermediate
414: Note:
415: If the matrix `Y` is missing some diagonal entries this routine can be very slow. To make it fast one should initially
416: fill the matrix so that all diagonal entries have a value (with a value of zero for those locations that would not have an
417: entry).
419: .seealso: [](ch_matrices), `Mat`, `MatShift()`, `MatScale()`, `MatDiagonalScale()`
420: @*/
421: PetscErrorCode MatDiagonalSet(Mat Y, Vec D, InsertMode is)
422: {
423: PetscInt matlocal, veclocal;
425: PetscFunctionBegin;
428: MatCheckPreallocated(Y, 1);
429: PetscCall(MatGetLocalSize(Y, &matlocal, NULL));
430: PetscCall(VecGetLocalSize(D, &veclocal));
431: 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);
432: if (Y->ops->diagonalset) PetscUseTypeMethod(Y, diagonalset, D, is);
433: else PetscCall(MatDiagonalSet_Default(Y, D, is));
434: PetscCall(PetscObjectStateIncrease((PetscObject)Y));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@
439: MatAYPX - Computes Y = a*Y + X.
441: Logically Collective
443: Input Parameters:
444: + a - the `PetscScalar` multiplier
445: . Y - the first matrix
446: . X - the second matrix
447: - str - either `SAME_NONZERO_PATTERN`, `DIFFERENT_NONZERO_PATTERN`, `UNKNOWN_NONZERO_PATTERN`, or `SUBSET_NONZERO_PATTERN` (nonzeros of `X` is a subset of `Y`'s)
449: Level: intermediate
451: Note:
452: `X` and `Y` may be the same matrix, in which case this computes `Y` = (`a` + 1) * `Y`. Passing an `X` that
453: merely shares storage with `Y` (for example a `MATTRANSPOSEVIRTUAL` wrapping `Y`) is not supported and
454: gives an undefined result.
456: .seealso: [](ch_matrices), `Mat`, `MatAXPY()`
457: @*/
458: PetscErrorCode MatAYPX(Mat Y, PetscScalar a, Mat X, MatStructure str)
459: {
460: PetscFunctionBegin;
461: if (Y == X) {
462: PetscCall(MatScale(Y, a + 1.0));
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
465: PetscCall(MatScale(Y, a));
466: PetscCall(MatAXPY(Y, 1.0, X, str));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*@
471: MatComputeOperator - Computes the explicit matrix
473: Collective
475: Input Parameters:
476: + inmat - the matrix
477: - mattype - the matrix type for the explicit operator
479: Output Parameter:
480: . mat - the explicit operator
482: Level: advanced
484: Note:
485: This computation is done by applying the operator to columns of the identity matrix.
486: This routine is costly in general, and is recommended for use only with relatively small systems.
487: Currently, this routine uses a dense matrix format if `mattype` == `NULL`.
489: .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperatorTranspose()`
490: @*/
491: PetscErrorCode MatComputeOperator(Mat inmat, MatType mattype, Mat *mat)
492: {
493: PetscFunctionBegin;
495: PetscAssertPointer(mat, 3);
496: PetscCall(MatConvert_Shell(inmat, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@
501: MatComputeOperatorTranspose - Computes the explicit matrix representation of
502: a give matrix that can apply `MatMultTranspose()`
504: Collective
506: Input Parameters:
507: + inmat - the matrix
508: - mattype - the matrix type for the explicit operator
510: Output Parameter:
511: . mat - the explicit operator transposed
513: Level: advanced
515: Note:
516: This computation is done by applying the transpose of the operator to columns of the identity matrix.
517: This routine is costly in general, and is recommended for use only with relatively small systems.
518: Currently, this routine uses a dense matrix format if `mattype` == `NULL`.
520: .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperator()`
521: @*/
522: PetscErrorCode MatComputeOperatorTranspose(Mat inmat, MatType mattype, Mat *mat)
523: {
524: Mat A;
526: PetscFunctionBegin;
528: PetscAssertPointer(mat, 3);
529: PetscCall(MatCreateTranspose(inmat, &A));
530: PetscCall(MatConvert_Shell(A, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
531: PetscCall(MatDestroy(&A));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /*@
536: 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
538: Input Parameters:
539: + A - The matrix
540: . tol - The zero tolerance
541: . compress - Whether the storage from the input matrix `A` should be compressed once values less than or equal to `tol` are set to zero
542: - 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
544: Level: intermediate
546: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatZeroEntries()`, `MatEliminateZeros()`, `VecFilter()`
547: @*/
548: PetscErrorCode MatFilter(Mat A, PetscReal tol, PetscBool compress, PetscBool keep)
549: {
550: Mat a;
551: PetscScalar *newVals;
552: PetscInt *newCols, rStart, rEnd, maxRows, r, colMax = 0, nnz0 = 0, nnz1 = 0;
553: PetscBool flg;
555: PetscFunctionBegin;
556: PetscCall(PetscObjectBaseTypeCompareAny((PetscObject)A, &flg, MATSEQDENSE, MATMPIDENSE, ""));
557: if (flg) {
558: PetscCall(MatDenseGetLocalMatrix(A, &a));
559: PetscCall(MatDenseGetLDA(a, &r));
560: PetscCall(MatGetSize(a, &rStart, &rEnd));
561: PetscCall(MatDenseGetArray(a, &newVals));
562: for (; colMax < rEnd; ++colMax) {
563: for (maxRows = 0; maxRows < rStart; ++maxRows) newVals[maxRows + colMax * r] = PetscAbsScalar(newVals[maxRows + colMax * r]) <= tol ? 0.0 : newVals[maxRows + colMax * r];
564: }
565: PetscCall(MatDenseRestoreArray(a, &newVals));
566: } else {
567: const PetscInt *ranges;
568: PetscMPIInt rank, size;
570: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)A), &rank));
571: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
572: PetscCall(MatGetOwnershipRanges(A, &ranges));
573: rStart = ranges[rank];
574: rEnd = ranges[rank + 1];
575: PetscCall(MatGetRowUpperTriangular(A));
576: for (r = rStart; r < rEnd; ++r) {
577: PetscInt ncols;
579: PetscCall(MatGetRow(A, r, &ncols, NULL, NULL));
580: colMax = PetscMax(colMax, ncols);
581: PetscCall(MatRestoreRow(A, r, &ncols, NULL, NULL));
582: }
583: maxRows = 0;
584: for (r = 0; r < size; ++r) maxRows = PetscMax(maxRows, ranges[r + 1] - ranges[r]);
585: PetscCall(PetscCalloc2(colMax, &newCols, colMax, &newVals));
586: PetscCall(MatGetOption(A, MAT_NO_OFF_PROC_ENTRIES, &flg)); /* cache user-defined value */
587: PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
588: /* short-circuit code in MatAssemblyBegin() and MatAssemblyEnd() */
589: /* that are potentially called many times depending on the distribution of A */
590: for (r = rStart; r < rStart + maxRows; ++r) {
591: if (r < rEnd) {
592: const PetscScalar *vals;
593: const PetscInt *cols;
594: PetscInt ncols, newcols = 0, c;
596: PetscCall(MatGetRow(A, r, &ncols, &cols, &vals));
597: nnz0 += ncols - 1;
598: for (c = 0; c < ncols; ++c) {
599: if (PetscUnlikely(PetscAbsScalar(vals[c]) <= tol)) newCols[newcols++] = cols[c];
600: }
601: nnz1 += ncols - newcols - 1;
602: PetscCall(MatRestoreRow(A, r, &ncols, &cols, &vals));
603: PetscCall(MatSetValues(A, 1, &r, newcols, newCols, newVals, INSERT_VALUES));
604: }
605: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
606: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
607: }
608: PetscCall(MatRestoreRowUpperTriangular(A));
609: PetscCall(PetscFree2(newCols, newVals));
610: PetscCall(MatSetOption(A, MAT_NO_OFF_PROC_ENTRIES, flg)); /* reset option to its user-defined value */
611: if (nnz0 > 0) PetscCall(PetscInfo(NULL, "Filtering left %g%% edges in graph\n", 100 * (double)nnz1 / (double)nnz0));
612: else PetscCall(PetscInfo(NULL, "Warning: %" PetscInt_FMT " edges to filter with %" PetscInt_FMT " rows\n", nnz0, maxRows));
613: }
614: if (compress && A->ops->eliminatezeros) {
615: Mat B;
616: PetscBool flg;
618: PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &flg, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
619: if (!flg) {
620: PetscCall(MatEliminateZeros(A, keep));
621: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
622: PetscCall(MatHeaderReplace(A, &B));
623: }
624: }
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }