1 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
2
MatTransposeAXPY_Private(Mat Y,PetscScalar a,Mat X,MatStructure str,Mat T)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 *);
8
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 }
42
MatAXPY_BasicWithTypeCompare(Mat Y,PetscScalar a,Mat X,MatStructure str)43 static PetscErrorCode MatAXPY_BasicWithTypeCompare(Mat Y, PetscScalar a, Mat X, MatStructure str)
44 {
45 PetscBool flg;
46
47 PetscFunctionBegin;
48 PetscCall(MatIsShell(Y, &flg));
49 if (flg) { /* MatShell has special support for AXPY */
50 PetscErrorCode (*f)(Mat, PetscScalar, Mat, MatStructure);
51
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;
70
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 }
80
81 /*@
82 MatAXPY - Computes Y = a*X + Y.
83
84 Logically Collective
85
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)
91
92 Level: intermediate
93
94 .seealso: [](ch_matrices), `Mat`, `MatAYPX()`
95 @*/
MatAXPY(Mat Y,PetscScalar a,Mat X,MatStructure str)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;
101
102 PetscFunctionBegin;
103 PetscValidHeaderSpecific(Y, MAT_CLASSID, 1);
104 PetscValidLogicalCollectiveScalar(Y, a, 2);
105 PetscValidHeaderSpecific(X, MAT_CLASSID, 3);
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 }
140
MatAXPY_Basic_Preallocate(Mat Y,Mat X,Mat * B)141 PetscErrorCode MatAXPY_Basic_Preallocate(Mat Y, Mat X, Mat *B)
142 {
143 PetscErrorCode (*preall)(Mat, Mat, Mat *) = NULL;
144
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;
154
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;
168
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));
181
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 }
190
MatAXPY_Basic(Mat Y,PetscScalar a,Mat X,MatStructure str)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;
196
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;
207
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;
241
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 }
248
MatAXPY_BasicWithPreallocation(Mat B,Mat Y,PetscScalar a,Mat X,MatStructure str)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;
256
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));
267
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));
280
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 }
301
302 /*@
303 MatShift - Computes `Y = Y + a I`, where `a` is a `PetscScalar`
304
305 Neighbor-wise Collective
306
307 Input Parameters:
308 + Y - the matrix
309 - a - the `PetscScalar`
310
311 Level: intermediate
312
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)
315
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.
319
320 To form Y = Y + diag(V) use `MatDiagonalSet()`
321
322 .seealso: [](ch_matrices), `Mat`, `MatDiagonalSet()`, `MatScale()`, `MatDiagonalScale()`
323 @*/
MatShift(Mat Y,PetscScalar a)324 PetscErrorCode MatShift(Mat Y, PetscScalar a)
325 {
326 PetscFunctionBegin;
327 PetscValidHeaderSpecific(Y, MAT_CLASSID, 1);
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);
332
333 if (Y->ops->shift) PetscUseTypeMethod(Y, shift, a);
334 else PetscCall(MatShift_Basic(Y, a));
335
336 PetscCall(PetscObjectStateIncrease((PetscObject)Y));
337 PetscFunctionReturn(PETSC_SUCCESS);
338 }
339
MatDiagonalSet_Default(Mat Y,Vec D,InsertMode is)340 PetscErrorCode MatDiagonalSet_Default(Mat Y, Vec D, InsertMode is)
341 {
342 PetscInt i, start, end;
343 const PetscScalar *v;
344
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 }
354
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`.
359
360 Neighbor-wise Collective
361
362 Input Parameters:
363 + Y - the input matrix
364 . D - the diagonal matrix, represented as a vector
365 - is - `INSERT_VALUES` or `ADD_VALUES`
366
367 Level: intermediate
368
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).
373
374 .seealso: [](ch_matrices), `Mat`, `MatShift()`, `MatScale()`, `MatDiagonalScale()`
375 @*/
MatDiagonalSet(Mat Y,Vec D,InsertMode is)376 PetscErrorCode MatDiagonalSet(Mat Y, Vec D, InsertMode is)
377 {
378 PetscInt matlocal, veclocal;
379
380 PetscFunctionBegin;
381 PetscValidHeaderSpecific(Y, MAT_CLASSID, 1);
382 PetscValidHeaderSpecific(D, VEC_CLASSID, 2);
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 }
392
393 /*@
394 MatAYPX - Computes Y = a*Y + X.
395
396 Logically Collective
397
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)
403
404 Level: intermediate
405
406 .seealso: [](ch_matrices), `Mat`, `MatAXPY()`
407 @*/
MatAYPX(Mat Y,PetscScalar a,Mat X,MatStructure str)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 }
415
416 /*@
417 MatComputeOperator - Computes the explicit matrix
418
419 Collective
420
421 Input Parameters:
422 + inmat - the matrix
423 - mattype - the matrix type for the explicit operator
424
425 Output Parameter:
426 . mat - the explicit operator
427
428 Level: advanced
429
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`.
434
435 .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperatorTranspose()`
436 @*/
MatComputeOperator(Mat inmat,MatType mattype,Mat * mat)437 PetscErrorCode MatComputeOperator(Mat inmat, MatType mattype, Mat *mat)
438 {
439 PetscFunctionBegin;
440 PetscValidHeaderSpecific(inmat, MAT_CLASSID, 1);
441 PetscAssertPointer(mat, 3);
442 PetscCall(MatConvert_Shell(inmat, mattype ? mattype : MATDENSE, MAT_INITIAL_MATRIX, mat));
443 PetscFunctionReturn(PETSC_SUCCESS);
444 }
445
446 /*@
447 MatComputeOperatorTranspose - Computes the explicit matrix representation of
448 a give matrix that can apply `MatMultTranspose()`
449
450 Collective
451
452 Input Parameters:
453 + inmat - the matrix
454 - mattype - the matrix type for the explicit operator
455
456 Output Parameter:
457 . mat - the explicit operator transposed
458
459 Level: advanced
460
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`.
465
466 .seealso: [](ch_matrices), `Mat`, `MatConvert()`, `MatMult()`, `MatComputeOperator()`
467 @*/
MatComputeOperatorTranspose(Mat inmat,MatType mattype,Mat * mat)468 PetscErrorCode MatComputeOperatorTranspose(Mat inmat, MatType mattype, Mat *mat)
469 {
470 Mat A;
471
472 PetscFunctionBegin;
473 PetscValidHeaderSpecific(inmat, MAT_CLASSID, 1);
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 }
480
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
483
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
489
490 Level: intermediate
491
492 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatZeroEntries()`, `MatEliminateZeros()`, `VecFilter()`
493 @*/
MatFilter(Mat A,PetscReal tol,PetscBool compress,PetscBool keep)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;
500
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;
515
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;
524
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;
541
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;
563
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 }
573