1 /*
2 This provides a simple shell for Fortran (and C programmers) to
3 create a very simple matrix class for use with KSP without coding
4 much of anything.
5 */
6
7 #include <../src/mat/impls/shell/shell.h> /*I "petscmat.h" I*/
8
9 /*
10 Store and scale values on zeroed rows
11 xx = [x_1, 0], 0 on zeroed columns
12 */
MatShellPreZeroRight(Mat A,Vec x,Vec * xx)13 static PetscErrorCode MatShellPreZeroRight(Mat A, Vec x, Vec *xx)
14 {
15 Mat_Shell *shell = (Mat_Shell *)A->data;
16
17 PetscFunctionBegin;
18 *xx = x;
19 if (shell->zrows) {
20 PetscCall(VecSet(shell->zvals_w, 0.0));
21 PetscCall(VecScatterBegin(shell->zvals_sct_c, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
22 PetscCall(VecScatterEnd(shell->zvals_sct_c, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
23 PetscCall(VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals));
24 }
25 if (shell->zcols) {
26 if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
27 PetscCall(VecCopy(x, shell->right_work));
28 PetscCall(VecISSet(shell->right_work, shell->zcols, 0.0));
29 *xx = shell->right_work;
30 }
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33
34 /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
MatShellPostZeroLeft(Mat A,Vec x)35 static PetscErrorCode MatShellPostZeroLeft(Mat A, Vec x)
36 {
37 Mat_Shell *shell = (Mat_Shell *)A->data;
38
39 PetscFunctionBegin;
40 if (shell->zrows) {
41 PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE));
42 PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, x, INSERT_VALUES, SCATTER_REVERSE));
43 }
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
47 /*
48 Store and scale values on zeroed rows
49 xx = [x_1, 0], 0 on zeroed rows
50 */
MatShellPreZeroLeft(Mat A,Vec x,Vec * xx)51 static PetscErrorCode MatShellPreZeroLeft(Mat A, Vec x, Vec *xx)
52 {
53 Mat_Shell *shell = (Mat_Shell *)A->data;
54
55 PetscFunctionBegin;
56 *xx = NULL;
57 if (!shell->zrows) {
58 *xx = x;
59 } else {
60 if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
61 PetscCall(VecCopy(x, shell->left_work));
62 PetscCall(VecSet(shell->zvals_w, 0.0));
63 PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE));
64 PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals_w, shell->left_work, INSERT_VALUES, SCATTER_REVERSE));
65 PetscCall(VecScatterBegin(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
66 PetscCall(VecScatterEnd(shell->zvals_sct_r, x, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
67 PetscCall(VecPointwiseMult(shell->zvals_w, shell->zvals_w, shell->zvals));
68 *xx = shell->left_work;
69 }
70 PetscFunctionReturn(PETSC_SUCCESS);
71 }
72
73 /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */
MatShellPostZeroRight(Mat A,Vec x)74 static PetscErrorCode MatShellPostZeroRight(Mat A, Vec x)
75 {
76 Mat_Shell *shell = (Mat_Shell *)A->data;
77
78 PetscFunctionBegin;
79 if (shell->zcols) PetscCall(VecISSet(x, shell->zcols, 0.0));
80 if (shell->zrows) {
81 PetscCall(VecScatterBegin(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE));
82 PetscCall(VecScatterEnd(shell->zvals_sct_c, shell->zvals_w, x, ADD_VALUES, SCATTER_REVERSE));
83 }
84 PetscFunctionReturn(PETSC_SUCCESS);
85 }
86
87 /*
88 xx = diag(left)*x
89 */
MatShellPreScaleLeft(Mat A,Vec x,Vec * xx,PetscBool conjugate)90 static PetscErrorCode MatShellPreScaleLeft(Mat A, Vec x, Vec *xx, PetscBool conjugate)
91 {
92 Mat_Shell *shell = (Mat_Shell *)A->data;
93
94 PetscFunctionBegin;
95 *xx = NULL;
96 if (!shell->left) {
97 *xx = x;
98 } else {
99 if (!shell->left_work) PetscCall(VecDuplicate(shell->left, &shell->left_work));
100 if (conjugate) { /* get arrays because there is no VecPointwiseMultConj() */
101 PetscInt i, m;
102 const PetscScalar *d, *xarray;
103 PetscScalar *w;
104 PetscCall(VecGetLocalSize(x, &m));
105 PetscCall(VecGetArrayRead(shell->left, &d));
106 PetscCall(VecGetArrayRead(x, &xarray));
107 PetscCall(VecGetArrayWrite(shell->left_work, &w));
108 for (i = 0; i < m; i++) w[i] = PetscConj(d[i]) * xarray[i];
109 PetscCall(VecRestoreArrayRead(shell->dshift, &d));
110 PetscCall(VecRestoreArrayRead(x, &xarray));
111 PetscCall(VecRestoreArrayWrite(shell->left_work, &w));
112 } else PetscCall(VecPointwiseMult(shell->left_work, x, shell->left));
113 *xx = shell->left_work;
114 }
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117
118 /*
119 xx = diag(right)*x
120 */
MatShellPreScaleRight(Mat A,Vec x,Vec * xx)121 static PetscErrorCode MatShellPreScaleRight(Mat A, Vec x, Vec *xx)
122 {
123 Mat_Shell *shell = (Mat_Shell *)A->data;
124
125 PetscFunctionBegin;
126 *xx = NULL;
127 if (!shell->right) {
128 *xx = x;
129 } else {
130 if (!shell->right_work) PetscCall(VecDuplicate(shell->right, &shell->right_work));
131 PetscCall(VecPointwiseMult(shell->right_work, x, shell->right));
132 *xx = shell->right_work;
133 }
134 PetscFunctionReturn(PETSC_SUCCESS);
135 }
136
137 /*
138 x = diag(left)*x
139 */
MatShellPostScaleLeft(Mat A,Vec x)140 static PetscErrorCode MatShellPostScaleLeft(Mat A, Vec x)
141 {
142 Mat_Shell *shell = (Mat_Shell *)A->data;
143
144 PetscFunctionBegin;
145 if (shell->left) PetscCall(VecPointwiseMult(x, x, shell->left));
146 PetscFunctionReturn(PETSC_SUCCESS);
147 }
148
149 /*
150 x = diag(right)*x
151 */
MatShellPostScaleRight(Mat A,Vec x,PetscBool conjugate)152 static PetscErrorCode MatShellPostScaleRight(Mat A, Vec x, PetscBool conjugate)
153 {
154 Mat_Shell *shell = (Mat_Shell *)A->data;
155
156 PetscFunctionBegin;
157 if (shell->right) {
158 if (conjugate) { /* get arrays because there is no VecPointwiseMultConj() */
159 PetscInt i, m;
160 const PetscScalar *d;
161 PetscScalar *xarray;
162 PetscCall(VecGetLocalSize(x, &m));
163 PetscCall(VecGetArrayRead(shell->right, &d));
164 PetscCall(VecGetArray(x, &xarray));
165 for (i = 0; i < m; i++) xarray[i] = PetscConj(d[i]) * xarray[i];
166 PetscCall(VecRestoreArrayRead(shell->dshift, &d));
167 PetscCall(VecRestoreArray(x, &xarray));
168 } else PetscCall(VecPointwiseMult(x, x, shell->right));
169 }
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
173 /*
174 Y = vscale*Y + diag(dshift)*X + vshift*X
175
176 On input Y already contains A*x
177
178 If conjugate=PETSC_TRUE then vscale, dshift, and vshift are conjugated
179 */
MatShellShiftAndScale(Mat A,Vec X,Vec Y,PetscBool conjugate)180 static PetscErrorCode MatShellShiftAndScale(Mat A, Vec X, Vec Y, PetscBool conjugate)
181 {
182 Mat_Shell *shell = (Mat_Shell *)A->data;
183 PetscScalar vscale = conjugate ? PetscConj(shell->vscale) : shell->vscale;
184 PetscScalar vshift = conjugate ? PetscConj(shell->vshift) : shell->vshift;
185
186 PetscFunctionBegin;
187 if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
188 PetscInt i, m;
189 const PetscScalar *x, *d;
190 PetscScalar *y;
191 PetscCall(VecGetLocalSize(X, &m));
192 PetscCall(VecGetArrayRead(shell->dshift, &d));
193 PetscCall(VecGetArrayRead(X, &x));
194 PetscCall(VecGetArray(Y, &y));
195 if (conjugate)
196 for (i = 0; i < m; i++) y[i] = vscale * y[i] + PetscConj(d[i]) * x[i];
197 else
198 for (i = 0; i < m; i++) y[i] = vscale * y[i] + d[i] * x[i];
199 PetscCall(VecRestoreArrayRead(shell->dshift, &d));
200 PetscCall(VecRestoreArrayRead(X, &x));
201 PetscCall(VecRestoreArray(Y, &y));
202 } else {
203 PetscCall(VecScale(Y, vscale));
204 }
205 if (vshift != 0.0) PetscCall(VecAXPY(Y, vshift, X)); /* if test is for non-square matrices */
206 PetscFunctionReturn(PETSC_SUCCESS);
207 }
208
MatShellGetContext_Shell(Mat mat,PetscCtxRt ctx)209 static PetscErrorCode MatShellGetContext_Shell(Mat mat, PetscCtxRt ctx)
210 {
211 Mat_Shell *shell = (Mat_Shell *)mat->data;
212
213 PetscFunctionBegin;
214 if (shell->ctxcontainer) PetscCall(PetscContainerGetPointer(shell->ctxcontainer, ctx));
215 else *(void **)ctx = NULL;
216 PetscFunctionReturn(PETSC_SUCCESS);
217 }
218
219 /*@
220 MatShellGetContext - Returns the user-provided context associated with a `MATSHELL` shell matrix.
221
222 Not Collective
223
224 Input Parameter:
225 . mat - the matrix, should have been created with `MatCreateShell()`
226
227 Output Parameter:
228 . ctx - the user provided context
229
230 Level: advanced
231
232 Fortran Notes:
233 This only works when the context is a Fortran derived type or a `PetscObject`. Declare `ctx` with
234 .vb
235 type(tUsertype), pointer :: ctx
236 .ve
237
238 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetOperation()`, `MatShellSetContext()`
239 @*/
MatShellGetContext(Mat mat,PetscCtxRt ctx)240 PetscErrorCode MatShellGetContext(Mat mat, PetscCtxRt ctx)
241 {
242 PetscFunctionBegin;
243 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
244 PetscAssertPointer(ctx, 2);
245 PetscUseMethod(mat, "MatShellGetContext_C", (Mat, void *), (mat, ctx));
246 PetscFunctionReturn(PETSC_SUCCESS);
247 }
248
MatZeroRowsColumns_Local_Shell(Mat mat,PetscInt nr,PetscInt rows[],PetscInt nc,PetscInt cols[],PetscScalar diag,PetscBool rc)249 static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat, PetscInt nr, PetscInt rows[], PetscInt nc, PetscInt cols[], PetscScalar diag, PetscBool rc)
250 {
251 Mat_Shell *shell = (Mat_Shell *)mat->data;
252 Vec x = NULL, b = NULL;
253 IS is1, is2;
254 const PetscInt *ridxs;
255 PetscInt *idxs, *gidxs;
256 PetscInt cum, rst, cst, i;
257
258 PetscFunctionBegin;
259 if (!shell->zvals) PetscCall(MatCreateVecs(mat, NULL, &shell->zvals));
260 if (!shell->zvals_w) PetscCall(VecDuplicate(shell->zvals, &shell->zvals_w));
261 PetscCall(MatGetOwnershipRange(mat, &rst, NULL));
262 PetscCall(MatGetOwnershipRangeColumn(mat, &cst, NULL));
263
264 /* Expand/create index set of zeroed rows */
265 PetscCall(PetscMalloc1(nr, &idxs));
266 for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
267 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), nr, idxs, PETSC_OWN_POINTER, &is1));
268 PetscCall(ISSort(is1));
269 PetscCall(VecISSet(shell->zvals, is1, diag));
270 if (shell->zrows) {
271 PetscCall(ISSum(shell->zrows, is1, &is2));
272 PetscCall(ISDestroy(&shell->zrows));
273 PetscCall(ISDestroy(&is1));
274 shell->zrows = is2;
275 } else shell->zrows = is1;
276
277 /* Create scatters for diagonal values communications */
278 PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
279 PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
280
281 /* row scatter: from/to left vector */
282 PetscCall(MatCreateVecs(mat, &x, &b));
283 PetscCall(VecScatterCreate(b, shell->zrows, shell->zvals_w, shell->zrows, &shell->zvals_sct_r));
284
285 /* col scatter: from right vector to left vector */
286 PetscCall(ISGetIndices(shell->zrows, &ridxs));
287 PetscCall(ISGetLocalSize(shell->zrows, &nr));
288 PetscCall(PetscMalloc1(nr, &gidxs));
289 for (i = 0, cum = 0; i < nr; i++) {
290 if (ridxs[i] >= mat->cmap->N) continue;
291 gidxs[cum] = ridxs[i];
292 cum++;
293 }
294 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), cum, gidxs, PETSC_OWN_POINTER, &is1));
295 PetscCall(VecScatterCreate(x, is1, shell->zvals_w, is1, &shell->zvals_sct_c));
296 PetscCall(ISDestroy(&is1));
297 PetscCall(VecDestroy(&x));
298 PetscCall(VecDestroy(&b));
299
300 /* Expand/create index set of zeroed columns */
301 if (rc) {
302 PetscCall(PetscMalloc1(nc, &idxs));
303 for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
304 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)mat), nc, idxs, PETSC_OWN_POINTER, &is1));
305 PetscCall(ISSort(is1));
306 if (shell->zcols) {
307 PetscCall(ISSum(shell->zcols, is1, &is2));
308 PetscCall(ISDestroy(&shell->zcols));
309 PetscCall(ISDestroy(&is1));
310 shell->zcols = is2;
311 } else shell->zcols = is1;
312 }
313 PetscFunctionReturn(PETSC_SUCCESS);
314 }
315
MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)316 static PetscErrorCode MatZeroRows_Shell(Mat mat, PetscInt n, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
317 {
318 Mat_Shell *shell = (Mat_Shell *)mat->data;
319 PetscInt nr, *lrows;
320
321 PetscFunctionBegin;
322 if (x && b) {
323 Vec xt;
324 PetscScalar *vals;
325 PetscInt *gcols, i, st, nl, nc;
326
327 PetscCall(PetscMalloc1(n, &gcols));
328 for (i = 0, nc = 0; i < n; i++)
329 if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];
330
331 PetscCall(MatCreateVecs(mat, &xt, NULL));
332 PetscCall(VecCopy(x, xt));
333 PetscCall(PetscCalloc1(nc, &vals));
334 PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
335 PetscCall(PetscFree(vals));
336 PetscCall(VecAssemblyBegin(xt));
337 PetscCall(VecAssemblyEnd(xt));
338 PetscCall(VecAYPX(xt, -1.0, x)); /* xt = [0, x2] */
339
340 PetscCall(VecGetOwnershipRange(xt, &st, NULL));
341 PetscCall(VecGetLocalSize(xt, &nl));
342 PetscCall(VecGetArray(xt, &vals));
343 for (i = 0; i < nl; i++) {
344 PetscInt g = i + st;
345 if (g > mat->rmap->N) continue;
346 if (PetscAbsScalar(vals[i]) == 0.0) continue;
347 PetscCall(VecSetValue(b, g, diag * vals[i], INSERT_VALUES));
348 }
349 PetscCall(VecRestoreArray(xt, &vals));
350 PetscCall(VecAssemblyBegin(b));
351 PetscCall(VecAssemblyEnd(b)); /* b = [b1, x2 * diag] */
352 PetscCall(VecDestroy(&xt));
353 PetscCall(PetscFree(gcols));
354 }
355 PetscCall(PetscLayoutMapLocal(mat->rmap, n, rows, &nr, &lrows, NULL));
356 PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, 0, NULL, diag, PETSC_FALSE));
357 if (shell->axpy) PetscCall(MatZeroRows(shell->axpy, n, rows, 0.0, NULL, NULL));
358 PetscCall(PetscFree(lrows));
359 PetscFunctionReturn(PETSC_SUCCESS);
360 }
361
MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b)362 static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat, PetscInt n, const PetscInt rowscols[], PetscScalar diag, Vec x, Vec b)
363 {
364 Mat_Shell *shell = (Mat_Shell *)mat->data;
365 PetscInt *lrows, *lcols;
366 PetscInt nr, nc;
367 PetscBool congruent;
368
369 PetscFunctionBegin;
370 if (x && b) {
371 Vec xt, bt;
372 PetscScalar *vals;
373 PetscInt *grows, *gcols, i, st, nl;
374
375 PetscCall(PetscMalloc2(n, &grows, n, &gcols));
376 for (i = 0, nr = 0; i < n; i++)
377 if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
378 for (i = 0, nc = 0; i < n; i++)
379 if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
380 PetscCall(PetscCalloc1(n, &vals));
381
382 PetscCall(MatCreateVecs(mat, &xt, &bt));
383 PetscCall(VecCopy(x, xt));
384 PetscCall(VecSetValues(xt, nc, gcols, vals, INSERT_VALUES)); /* xt = [x1, 0] */
385 PetscCall(VecAssemblyBegin(xt));
386 PetscCall(VecAssemblyEnd(xt));
387 PetscCall(VecAXPY(xt, -1.0, x)); /* xt = [0, -x2] */
388 PetscCall(MatMult(mat, xt, bt)); /* bt = [-A12*x2,-A22*x2] */
389 PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* bt = [-A12*x2,0] */
390 PetscCall(VecAssemblyBegin(bt));
391 PetscCall(VecAssemblyEnd(bt));
392 PetscCall(VecAXPY(b, 1.0, bt)); /* b = [b1 - A12*x2, b2] */
393 PetscCall(VecSetValues(bt, nr, grows, vals, INSERT_VALUES)); /* b = [b1 - A12*x2, 0] */
394 PetscCall(VecAssemblyBegin(bt));
395 PetscCall(VecAssemblyEnd(bt));
396 PetscCall(PetscFree(vals));
397
398 PetscCall(VecGetOwnershipRange(xt, &st, NULL));
399 PetscCall(VecGetLocalSize(xt, &nl));
400 PetscCall(VecGetArray(xt, &vals));
401 for (i = 0; i < nl; i++) {
402 PetscInt g = i + st;
403 if (g > mat->rmap->N) continue;
404 if (PetscAbsScalar(vals[i]) == 0.0) continue;
405 PetscCall(VecSetValue(b, g, -diag * vals[i], INSERT_VALUES));
406 }
407 PetscCall(VecRestoreArray(xt, &vals));
408 PetscCall(VecAssemblyBegin(b));
409 PetscCall(VecAssemblyEnd(b)); /* b = [b1 - A12*x2, x2 * diag] */
410 PetscCall(VecDestroy(&xt));
411 PetscCall(VecDestroy(&bt));
412 PetscCall(PetscFree2(grows, gcols));
413 }
414 PetscCall(PetscLayoutMapLocal(mat->rmap, n, rowscols, &nr, &lrows, NULL));
415 PetscCall(MatHasCongruentLayouts(mat, &congruent));
416 if (congruent) {
417 nc = nr;
418 lcols = lrows;
419 } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
420 PetscInt i, nt, *t;
421
422 PetscCall(PetscMalloc1(n, &t));
423 for (i = 0, nt = 0; i < n; i++)
424 if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
425 PetscCall(PetscLayoutMapLocal(mat->cmap, nt, t, &nc, &lcols, NULL));
426 PetscCall(PetscFree(t));
427 }
428 PetscCall(MatZeroRowsColumns_Local_Shell(mat, nr, lrows, nc, lcols, diag, PETSC_TRUE));
429 if (!congruent) PetscCall(PetscFree(lcols));
430 PetscCall(PetscFree(lrows));
431 if (shell->axpy) PetscCall(MatZeroRowsColumns(shell->axpy, n, rowscols, 0.0, NULL, NULL));
432 PetscFunctionReturn(PETSC_SUCCESS);
433 }
434
MatDestroy_Shell(Mat mat)435 static PetscErrorCode MatDestroy_Shell(Mat mat)
436 {
437 Mat_Shell *shell = (Mat_Shell *)mat->data;
438 MatShellMatFunctionList matmat;
439
440 PetscFunctionBegin;
441 if (shell->ops->destroy) PetscCall((*shell->ops->destroy)(mat));
442 PetscCall(PetscMemzero(shell->ops, sizeof(struct _MatShellOps)));
443 PetscCall(VecDestroy(&shell->left));
444 PetscCall(VecDestroy(&shell->right));
445 PetscCall(VecDestroy(&shell->dshift));
446 PetscCall(VecDestroy(&shell->left_work));
447 PetscCall(VecDestroy(&shell->right_work));
448 PetscCall(VecDestroy(&shell->left_add_work));
449 PetscCall(VecDestroy(&shell->right_add_work));
450 PetscCall(VecDestroy(&shell->axpy_left));
451 PetscCall(VecDestroy(&shell->axpy_right));
452 PetscCall(MatDestroy(&shell->axpy));
453 PetscCall(VecDestroy(&shell->zvals_w));
454 PetscCall(VecDestroy(&shell->zvals));
455 PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
456 PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
457 PetscCall(ISDestroy(&shell->zrows));
458 PetscCall(ISDestroy(&shell->zcols));
459
460 matmat = shell->matmat;
461 while (matmat) {
462 MatShellMatFunctionList next = matmat->next;
463
464 PetscCall(PetscObjectComposeFunction((PetscObject)mat, matmat->composedname, NULL));
465 PetscCall(PetscFree(matmat->composedname));
466 PetscCall(PetscFree(matmat->resultname));
467 PetscCall(PetscFree(matmat));
468 matmat = next;
469 }
470 PetscCall(MatShellSetContext(mat, NULL));
471 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetContext_C", NULL));
472 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContext_C", NULL));
473 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetContextDestroy_C", NULL));
474 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetVecType_C", NULL));
475 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetManageScalingShifts_C", NULL));
476 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetScalingShifts_C", NULL));
477 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetOperation_C", NULL));
478 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellGetOperation_C", NULL));
479 PetscCall(PetscObjectComposeFunction((PetscObject)mat, "MatShellSetMatProductOperation_C", NULL));
480 PetscCall(PetscFree(mat->data));
481 PetscFunctionReturn(PETSC_SUCCESS);
482 }
483
484 typedef struct {
485 PetscErrorCode (*numeric)(Mat, Mat, Mat, void *);
486 PetscCtxDestroyFn *destroy;
487 void *ctx;
488 Mat B;
489 Mat Bt;
490 Mat axpy;
491 } MatProductCtx_MatMatShell;
492
MatProductCtxDestroy_MatMatShell(PetscCtxRt data)493 static PetscErrorCode MatProductCtxDestroy_MatMatShell(PetscCtxRt data)
494 {
495 MatProductCtx_MatMatShell *mmdata = *(MatProductCtx_MatMatShell **)data;
496
497 PetscFunctionBegin;
498 if (mmdata->destroy) PetscCall((*mmdata->destroy)(&mmdata->ctx));
499 PetscCall(MatDestroy(&mmdata->B));
500 PetscCall(MatDestroy(&mmdata->Bt));
501 PetscCall(MatDestroy(&mmdata->axpy));
502 PetscCall(PetscFree(mmdata));
503 PetscFunctionReturn(PETSC_SUCCESS);
504 }
505
MatProductNumeric_Shell_X(Mat D)506 static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
507 {
508 Mat_Product *product;
509 Mat A, B;
510 MatProductCtx_MatMatShell *mdata;
511 Mat_Shell *shell;
512 PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE;
513 PetscErrorCode (*stashsym)(Mat), (*stashnum)(Mat);
514
515 PetscFunctionBegin;
516 MatCheckProduct(D, 1);
517 product = D->product;
518 PetscCheck(product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data empty");
519 A = product->A;
520 B = product->B;
521 mdata = (MatProductCtx_MatMatShell *)product->data;
522 PetscCheck(mdata->numeric, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Missing numeric operation");
523 shell = (Mat_Shell *)A->data;
524 stashsym = D->ops->productsymbolic;
525 stashnum = D->ops->productnumeric;
526 if (shell->managescalingshifts) {
527 PetscCheck(!shell->zcols && !shell->zrows, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProduct not supported with zeroed rows/columns");
528 if (shell->right || shell->left) {
529 useBmdata = PETSC_TRUE;
530 if (!mdata->B) {
531 PetscCall(MatDuplicate(B, MAT_SHARE_NONZERO_PATTERN, &mdata->B));
532 } else {
533 newB = PETSC_FALSE;
534 }
535 PetscCall(MatCopy(B, mdata->B, SAME_NONZERO_PATTERN));
536 }
537 switch (product->type) {
538 case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
539 if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
540 break;
541 case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
542 if (shell->left) PetscCall(MatDiagonalScale(mdata->B, shell->left, NULL));
543 break;
544 case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
545 if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
546 break;
547 case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
548 if (shell->right && shell->left) {
549 PetscBool flg;
550
551 PetscCall(VecEqual(shell->right, shell->left, &flg));
552 PetscCheck(flg, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling", MatProductTypes[product->type], ((PetscObject)A)->type_name,
553 ((PetscObject)B)->type_name);
554 }
555 if (shell->right) PetscCall(MatDiagonalScale(mdata->B, NULL, shell->right));
556 break;
557 case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
558 if (shell->right && shell->left) {
559 PetscBool flg;
560
561 PetscCall(VecEqual(shell->right, shell->left, &flg));
562 PetscCheck(flg, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling", MatProductTypes[product->type], ((PetscObject)A)->type_name,
563 ((PetscObject)B)->type_name);
564 }
565 if (shell->right) PetscCall(MatDiagonalScale(mdata->B, shell->right, NULL));
566 break;
567 default:
568 SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
569 }
570 }
571 /* allow the user to call MatMat operations on D */
572 D->product = NULL;
573 D->ops->productsymbolic = NULL;
574 D->ops->productnumeric = NULL;
575
576 PetscCall((*mdata->numeric)(A, useBmdata ? mdata->B : B, D, mdata->ctx));
577
578 /* clear any leftover user data and restore D pointers */
579 PetscCall(MatProductClear(D));
580 D->ops->productsymbolic = stashsym;
581 D->ops->productnumeric = stashnum;
582 D->product = product;
583
584 if (shell->managescalingshifts) {
585 PetscCall(MatScale(D, shell->vscale));
586 switch (product->type) {
587 case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
588 case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
589 if (shell->left) {
590 PetscCall(MatDiagonalScale(D, shell->left, NULL));
591 if (shell->dshift || shell->vshift != 0.0) {
592 if (!shell->left_work) PetscCall(MatCreateVecs(A, NULL, &shell->left_work));
593 if (shell->dshift) {
594 PetscCall(VecCopy(shell->dshift, shell->left_work));
595 PetscCall(VecShift(shell->left_work, shell->vshift));
596 PetscCall(VecPointwiseMult(shell->left_work, shell->left_work, shell->left));
597 } else {
598 PetscCall(VecSet(shell->left_work, shell->vshift));
599 }
600 if (product->type == MATPRODUCT_ABt) {
601 MatReuse reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
602 MatStructure str = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;
603
604 PetscCall(MatTranspose(mdata->B, reuse, &mdata->Bt));
605 PetscCall(MatDiagonalScale(mdata->Bt, shell->left_work, NULL));
606 PetscCall(MatAXPY(D, 1.0, mdata->Bt, str));
607 } else {
608 MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
609
610 PetscCall(MatDiagonalScale(mdata->B, shell->left_work, NULL));
611 PetscCall(MatAXPY(D, 1.0, mdata->B, str));
612 }
613 }
614 }
615 break;
616 case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
617 if (shell->right) {
618 PetscCall(MatDiagonalScale(D, shell->right, NULL));
619 if (shell->dshift || shell->vshift != 0.0) {
620 MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
621
622 if (!shell->right_work) PetscCall(MatCreateVecs(A, &shell->right_work, NULL));
623 if (shell->dshift) {
624 PetscCall(VecCopy(shell->dshift, shell->right_work));
625 PetscCall(VecShift(shell->right_work, shell->vshift));
626 PetscCall(VecPointwiseMult(shell->right_work, shell->right_work, shell->right));
627 } else {
628 PetscCall(VecSet(shell->right_work, shell->vshift));
629 }
630 PetscCall(MatDiagonalScale(mdata->B, shell->right_work, NULL));
631 PetscCall(MatAXPY(D, 1.0, mdata->B, str));
632 }
633 }
634 break;
635 case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
636 case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
637 PetscCheck(!shell->dshift && shell->vshift == 0.0, PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices with diagonal shift", MatProductTypes[product->type], ((PetscObject)A)->type_name,
638 ((PetscObject)B)->type_name);
639 break;
640 default:
641 SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
642 }
643 if (shell->axpy && shell->axpy_vscale != 0.0) {
644 Mat X;
645 PetscObjectState axpy_state;
646 MatStructure str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */
647
648 PetscCall(MatShellGetContext(shell->axpy, &X));
649 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
650 PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
651 if (!mdata->axpy) {
652 str = DIFFERENT_NONZERO_PATTERN;
653 PetscCall(MatProductCreate(shell->axpy, B, NULL, &mdata->axpy));
654 PetscCall(MatProductSetType(mdata->axpy, product->type));
655 PetscCall(MatProductSetFromOptions(mdata->axpy));
656 PetscCall(MatProductSymbolic(mdata->axpy));
657 } else { /* May be that shell->axpy has changed */
658 PetscBool flg;
659
660 PetscCall(MatProductReplaceMats(shell->axpy, B, NULL, mdata->axpy));
661 PetscCall(MatHasOperation(mdata->axpy, MATOP_PRODUCTSYMBOLIC, &flg));
662 if (!flg) {
663 str = DIFFERENT_NONZERO_PATTERN;
664 PetscCall(MatProductSetFromOptions(mdata->axpy));
665 PetscCall(MatProductSymbolic(mdata->axpy));
666 }
667 }
668 PetscCall(MatProductNumeric(mdata->axpy));
669 PetscCall(MatAXPY(D, shell->axpy_vscale, mdata->axpy, str));
670 }
671 }
672 PetscFunctionReturn(PETSC_SUCCESS);
673 }
674
MatProductSymbolic_Shell_X(Mat D)675 static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
676 {
677 Mat_Product *product;
678 Mat A, B;
679 MatShellMatFunctionList matmat;
680 Mat_Shell *shell;
681 PetscBool flg = PETSC_FALSE;
682 char composedname[256];
683 MatProductCtx_MatMatShell *mdata;
684
685 PetscFunctionBegin;
686 MatCheckProduct(D, 1);
687 product = D->product;
688 PetscCheck(!product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "Product data not empty");
689 A = product->A;
690 B = product->B;
691 shell = (Mat_Shell *)A->data;
692 matmat = shell->matmat;
693 PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
694 while (matmat) {
695 PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
696 flg = (PetscBool)(flg && (matmat->ptype == product->type));
697 if (flg) break;
698 matmat = matmat->next;
699 }
700 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Composedname \"%s\" for product type %s not found", composedname, MatProductTypes[product->type]);
701 switch (product->type) {
702 case MATPRODUCT_AB:
703 PetscCall(MatSetSizes(D, A->rmap->n, B->cmap->n, A->rmap->N, B->cmap->N));
704 break;
705 case MATPRODUCT_AtB:
706 PetscCall(MatSetSizes(D, A->cmap->n, B->cmap->n, A->cmap->N, B->cmap->N));
707 break;
708 case MATPRODUCT_ABt:
709 PetscCall(MatSetSizes(D, A->rmap->n, B->rmap->n, A->rmap->N, B->rmap->N));
710 break;
711 case MATPRODUCT_RARt:
712 PetscCall(MatSetSizes(D, B->rmap->n, B->rmap->n, B->rmap->N, B->rmap->N));
713 break;
714 case MATPRODUCT_PtAP:
715 PetscCall(MatSetSizes(D, B->cmap->n, B->cmap->n, B->cmap->N, B->cmap->N));
716 break;
717 default:
718 SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "MatProductSymbolic type %s not supported for %s and %s matrices", MatProductTypes[product->type], ((PetscObject)A)->type_name, ((PetscObject)B)->type_name);
719 }
720 /* respect users who passed in a matrix for which resultname is the base type */
721 if (matmat->resultname) {
722 PetscCall(PetscObjectBaseTypeCompare((PetscObject)D, matmat->resultname, &flg));
723 if (!flg) PetscCall(MatSetType(D, matmat->resultname));
724 }
725 /* If matrix type was not set or different, we need to reset this pointers */
726 D->ops->productsymbolic = MatProductSymbolic_Shell_X;
727 D->ops->productnumeric = MatProductNumeric_Shell_X;
728 /* attach product data */
729 PetscCall(PetscNew(&mdata));
730 mdata->numeric = matmat->numeric;
731 mdata->destroy = matmat->destroy;
732 if (matmat->symbolic) {
733 PetscCall((*matmat->symbolic)(A, B, D, &mdata->ctx));
734 } else { /* call general setup if symbolic operation not provided */
735 PetscCall(MatSetUp(D));
736 }
737 PetscCheck(D->product, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product disappeared after user symbolic phase");
738 PetscCheck(!D->product->data, PetscObjectComm((PetscObject)D), PETSC_ERR_COR, "Product data not empty after user symbolic phase");
739 D->product->data = mdata;
740 D->product->destroy = MatProductCtxDestroy_MatMatShell;
741 /* Be sure to reset these pointers if the user did something unexpected */
742 D->ops->productsymbolic = MatProductSymbolic_Shell_X;
743 D->ops->productnumeric = MatProductNumeric_Shell_X;
744 PetscFunctionReturn(PETSC_SUCCESS);
745 }
746
MatProductSetFromOptions_Shell_X(Mat D)747 static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
748 {
749 Mat_Product *product;
750 Mat A, B;
751 MatShellMatFunctionList matmat;
752 Mat_Shell *shell;
753 PetscBool flg;
754 char composedname[256];
755
756 PetscFunctionBegin;
757 MatCheckProduct(D, 1);
758 product = D->product;
759 A = product->A;
760 B = product->B;
761 PetscCall(MatIsShell(A, &flg));
762 if (!flg) PetscFunctionReturn(PETSC_SUCCESS);
763 shell = (Mat_Shell *)A->data;
764 matmat = shell->matmat;
765 PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, ((PetscObject)B)->type_name));
766 while (matmat) {
767 PetscCall(PetscStrcmp(composedname, matmat->composedname, &flg));
768 flg = (PetscBool)(flg && (matmat->ptype == product->type));
769 if (flg) break;
770 matmat = matmat->next;
771 }
772 if (flg) {
773 D->ops->productsymbolic = MatProductSymbolic_Shell_X;
774 } else PetscCall(PetscInfo(D, " symbolic product %s not registered for product type %s\n", composedname, MatProductTypes[product->type]));
775 PetscFunctionReturn(PETSC_SUCCESS);
776 }
777
MatShellSetMatProductOperation_Private(Mat A,MatProductType ptype,PetscErrorCode (* symbolic)(Mat,Mat,Mat,void **),PetscErrorCode (* numeric)(Mat,Mat,Mat,void *),PetscCtxDestroyFn * destroy,char * composedname,const char * resultname)778 static PetscErrorCode MatShellSetMatProductOperation_Private(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscCtxDestroyFn *destroy, char *composedname, const char *resultname)
779 {
780 PetscBool flg;
781 Mat_Shell *shell;
782 MatShellMatFunctionList matmat;
783
784 PetscFunctionBegin;
785 PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing numeric routine");
786 PetscCheck(composedname, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing composed name");
787
788 /* add product callback */
789 shell = (Mat_Shell *)A->data;
790 matmat = shell->matmat;
791 if (!matmat) {
792 PetscCall(PetscNew(&shell->matmat));
793 matmat = shell->matmat;
794 } else {
795 MatShellMatFunctionList entry = matmat;
796 while (entry) {
797 PetscCall(PetscStrcmp(composedname, entry->composedname, &flg));
798 flg = (PetscBool)(flg && (entry->ptype == ptype));
799 matmat = entry;
800 if (flg) goto set;
801 entry = entry->next;
802 }
803 PetscCall(PetscNew(&matmat->next));
804 matmat = matmat->next;
805 }
806
807 set:
808 matmat->symbolic = symbolic;
809 matmat->numeric = numeric;
810 matmat->destroy = destroy;
811 matmat->ptype = ptype;
812 PetscCall(PetscFree(matmat->composedname));
813 PetscCall(PetscFree(matmat->resultname));
814 PetscCall(PetscStrallocpy(composedname, &matmat->composedname));
815 PetscCall(PetscStrallocpy(resultname, &matmat->resultname));
816 PetscCall(PetscInfo(A, "Composing %s for product type %s with result %s\n", matmat->composedname, MatProductTypes[matmat->ptype], matmat->resultname ? matmat->resultname : "not specified"));
817 PetscCall(PetscObjectComposeFunction((PetscObject)A, matmat->composedname, MatProductSetFromOptions_Shell_X));
818 PetscFunctionReturn(PETSC_SUCCESS);
819 }
820
821 /*@C
822 MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a `MATSHELL` shell matrix.
823
824 Logically Collective; No Fortran Support
825
826 Input Parameters:
827 + A - the `MATSHELL` shell matrix
828 . ptype - the product type
829 . symbolic - the function for the symbolic phase (can be `NULL`)
830 . numeric - the function for the numerical phase
831 . destroy - the function for the destruction of the needed data generated during the symbolic phase (can be `NULL`)
832 . Btype - the matrix type for the matrix to be multiplied against
833 - Ctype - the matrix type for the result (can be `NULL`)
834
835 Level: advanced
836
837 Example Usage:
838 .vb
839 extern PetscErrorCode usersymbolic(Mat, Mat, Mat, void**);
840 extern PetscErrorCode usernumeric(Mat, Mat, Mat, void*);
841 PetscCtxDestroyFn *ctxdestroy;
842
843 MatCreateShell(comm, m, n, M, N, ctx, &A);
844 MatShellSetMatProductOperation(
845 A, MATPRODUCT_AB, usersymbolic, usernumeric, ctxdestroy,MATSEQAIJ, MATDENSE
846 );
847 // create B of type SEQAIJ etc..
848 MatProductCreate(A, B, PETSC_NULLPTR, &C);
849 MatProductSetType(C, MATPRODUCT_AB);
850 MatProductSetFromOptions(C);
851 MatProductSymbolic(C); // actually runs the user defined symbolic operation
852 MatProductNumeric(C); // actually runs the user defined numeric operation
853 // use C = A * B
854 .ve
855
856 Notes:
857 `MATPRODUCT_ABC` is not supported yet.
858
859 If the symbolic phase is not specified, `MatSetUp()` is called on the result matrix that must have its type set if `Ctype` is `NULL`.
860
861 Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
862 PETSc will take care of calling the user-defined callbacks.
863 It is allowed to specify the same callbacks for different `Btype` matrix types.
864 The couple (`Btype`,`ptype`) uniquely identifies the operation, the last specified callbacks takes precedence.
865
866 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatProductType`, `MatType`, `MatSetUp()`
867 @*/
MatShellSetMatProductOperation(Mat A,MatProductType ptype,PetscErrorCode (* symbolic)(Mat,Mat,Mat,void **),PetscErrorCode (* numeric)(Mat,Mat,Mat,void *),PetscCtxDestroyFn * destroy,MatType Btype,MatType Ctype)868 PetscErrorCode MatShellSetMatProductOperation(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscCtxDestroyFn *destroy, MatType Btype, MatType Ctype)
869 {
870 PetscFunctionBegin;
871 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
872 PetscValidLogicalCollectiveEnum(A, ptype, 2);
873 PetscCheck(ptype != MATPRODUCT_ABC, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for product type %s", MatProductTypes[ptype]);
874 PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Missing numeric routine, argument 4");
875 PetscAssertPointer(Btype, 6);
876 if (Ctype) PetscAssertPointer(Ctype, 7);
877 PetscTryMethod(A, "MatShellSetMatProductOperation_C", (Mat, MatProductType, PetscErrorCode (*)(Mat, Mat, Mat, void **), PetscErrorCode (*)(Mat, Mat, Mat, void *), PetscCtxDestroyFn *, MatType, MatType), (A, ptype, symbolic, numeric, destroy, Btype, Ctype));
878 PetscFunctionReturn(PETSC_SUCCESS);
879 }
880
MatShellSetMatProductOperation_Shell(Mat A,MatProductType ptype,PetscErrorCode (* symbolic)(Mat,Mat,Mat,void **),PetscErrorCode (* numeric)(Mat,Mat,Mat,void *),PetscCtxDestroyFn * destroy,MatType Btype,MatType Ctype)881 static PetscErrorCode MatShellSetMatProductOperation_Shell(Mat A, MatProductType ptype, PetscErrorCode (*symbolic)(Mat, Mat, Mat, void **), PetscErrorCode (*numeric)(Mat, Mat, Mat, void *), PetscCtxDestroyFn *destroy, MatType Btype, MatType Ctype)
882 {
883 PetscBool flg;
884 char composedname[256];
885 MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
886 PetscMPIInt size;
887
888 PetscFunctionBegin;
889 PetscValidType(A, 1);
890 while (Bnames) { /* user passed in the root name */
891 PetscCall(PetscStrcmp(Btype, Bnames->rname, &flg));
892 if (flg) break;
893 Bnames = Bnames->next;
894 }
895 while (Cnames) { /* user passed in the root name */
896 PetscCall(PetscStrcmp(Ctype, Cnames->rname, &flg));
897 if (flg) break;
898 Cnames = Cnames->next;
899 }
900 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
901 Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
902 Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
903 PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, Btype));
904 PetscCall(MatShellSetMatProductOperation_Private(A, ptype, symbolic, numeric, destroy, composedname, Ctype));
905 PetscFunctionReturn(PETSC_SUCCESS);
906 }
907
MatCopy_Shell(Mat A,Mat B,MatStructure str)908 static PetscErrorCode MatCopy_Shell(Mat A, Mat B, MatStructure str)
909 {
910 Mat_Shell *shellA = (Mat_Shell *)A->data, *shellB = (Mat_Shell *)B->data;
911 PetscBool matflg;
912 MatShellMatFunctionList matmatA;
913
914 PetscFunctionBegin;
915 PetscCall(MatIsShell(B, &matflg));
916 PetscCheck(matflg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix %s not derived from MATSHELL", ((PetscObject)B)->type_name);
917
918 B->ops[0] = A->ops[0];
919 shellB->ops[0] = shellA->ops[0];
920
921 if (shellA->ops->copy) PetscCall((*shellA->ops->copy)(A, B, str));
922 shellB->vscale = shellA->vscale;
923 shellB->vshift = shellA->vshift;
924 if (shellA->dshift) {
925 if (!shellB->dshift) PetscCall(VecDuplicate(shellA->dshift, &shellB->dshift));
926 PetscCall(VecCopy(shellA->dshift, shellB->dshift));
927 } else {
928 PetscCall(VecDestroy(&shellB->dshift));
929 }
930 if (shellA->left) {
931 if (!shellB->left) PetscCall(VecDuplicate(shellA->left, &shellB->left));
932 PetscCall(VecCopy(shellA->left, shellB->left));
933 } else {
934 PetscCall(VecDestroy(&shellB->left));
935 }
936 if (shellA->right) {
937 if (!shellB->right) PetscCall(VecDuplicate(shellA->right, &shellB->right));
938 PetscCall(VecCopy(shellA->right, shellB->right));
939 } else {
940 PetscCall(VecDestroy(&shellB->right));
941 }
942 PetscCall(MatDestroy(&shellB->axpy));
943 shellB->axpy_vscale = 0.0;
944 shellB->axpy_state = 0;
945 if (shellA->axpy) {
946 PetscCall(PetscObjectReference((PetscObject)shellA->axpy));
947 shellB->axpy = shellA->axpy;
948 shellB->axpy_vscale = shellA->axpy_vscale;
949 shellB->axpy_state = shellA->axpy_state;
950 }
951 if (shellA->zrows) {
952 PetscCall(ISDuplicate(shellA->zrows, &shellB->zrows));
953 if (shellA->zcols) PetscCall(ISDuplicate(shellA->zcols, &shellB->zcols));
954 PetscCall(VecDuplicate(shellA->zvals, &shellB->zvals));
955 PetscCall(VecCopy(shellA->zvals, shellB->zvals));
956 PetscCall(VecDuplicate(shellA->zvals_w, &shellB->zvals_w));
957 PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_r));
958 PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_c));
959 shellB->zvals_sct_r = shellA->zvals_sct_r;
960 shellB->zvals_sct_c = shellA->zvals_sct_c;
961 }
962
963 matmatA = shellA->matmat;
964 if (matmatA) {
965 while (matmatA->next) {
966 PetscCall(MatShellSetMatProductOperation_Private(B, matmatA->ptype, matmatA->symbolic, matmatA->numeric, matmatA->destroy, matmatA->composedname, matmatA->resultname));
967 matmatA = matmatA->next;
968 }
969 }
970 PetscFunctionReturn(PETSC_SUCCESS);
971 }
972
MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat * M)973 static PetscErrorCode MatDuplicate_Shell(Mat mat, MatDuplicateOption op, Mat *M)
974 {
975 PetscFunctionBegin;
976 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)mat), mat->rmap->n, mat->cmap->n, mat->rmap->N, mat->cmap->N, NULL, M));
977 ((Mat_Shell *)(*M)->data)->ctxcontainer = ((Mat_Shell *)mat->data)->ctxcontainer;
978 PetscCall(PetscObjectCompose((PetscObject)*M, "MatShell ctx", (PetscObject)((Mat_Shell *)(*M)->data)->ctxcontainer));
979 PetscCall(PetscObjectChangeTypeName((PetscObject)*M, ((PetscObject)mat)->type_name));
980 if (op == MAT_COPY_VALUES) PetscCall(MatCopy(mat, *M, SAME_NONZERO_PATTERN));
981 PetscCall(PetscObjectCopyFortranFunctionPointers((PetscObject)mat, (PetscObject)*M));
982 PetscFunctionReturn(PETSC_SUCCESS);
983 }
984
MatMult_Shell(Mat A,Vec x,Vec y)985 static PetscErrorCode MatMult_Shell(Mat A, Vec x, Vec y)
986 {
987 Mat_Shell *shell = (Mat_Shell *)A->data;
988 Vec xx;
989 PetscObjectState instate, outstate;
990
991 PetscFunctionBegin;
992 PetscCall(MatShellPreZeroRight(A, x, &xx));
993 PetscCall(MatShellPreScaleRight(A, xx, &xx));
994 PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
995 PetscCall((*shell->ops->mult)(A, xx, y));
996 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
997 if (instate == outstate) {
998 /* increase the state of the output vector since the user did not update its state themself as should have been done */
999 PetscCall(PetscObjectStateIncrease((PetscObject)y));
1000 }
1001 PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_FALSE));
1002 PetscCall(MatShellPostScaleLeft(A, y));
1003 PetscCall(MatShellPostZeroLeft(A, y));
1004
1005 if (shell->axpy) {
1006 Mat X;
1007 PetscObjectState axpy_state;
1008
1009 PetscCall(MatShellGetContext(shell->axpy, &X));
1010 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1011 PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
1012
1013 PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1014 PetscCall(VecCopy(x, shell->axpy_right));
1015 PetscCall(MatMult(shell->axpy, shell->axpy_right, shell->axpy_left));
1016 PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_left));
1017 }
1018 PetscFunctionReturn(PETSC_SUCCESS);
1019 }
1020
MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)1021 static PetscErrorCode MatMultAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1022 {
1023 Mat_Shell *shell = (Mat_Shell *)A->data;
1024
1025 PetscFunctionBegin;
1026 if (y == z) {
1027 if (!shell->right_add_work) PetscCall(VecDuplicate(z, &shell->right_add_work));
1028 PetscCall(MatMult(A, x, shell->right_add_work));
1029 PetscCall(VecAXPY(z, 1.0, shell->right_add_work));
1030 } else {
1031 PetscCall(MatMult(A, x, z));
1032 PetscCall(VecAXPY(z, 1.0, y));
1033 }
1034 PetscFunctionReturn(PETSC_SUCCESS);
1035 }
1036
MatMultTranspose_Shell(Mat A,Vec x,Vec y)1037 static PetscErrorCode MatMultTranspose_Shell(Mat A, Vec x, Vec y)
1038 {
1039 Mat_Shell *shell = (Mat_Shell *)A->data;
1040 Vec xx;
1041 PetscObjectState instate, outstate;
1042
1043 PetscFunctionBegin;
1044 PetscCall(MatShellPreZeroLeft(A, x, &xx));
1045 PetscCall(MatShellPreScaleLeft(A, xx, &xx, PETSC_FALSE));
1046 PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
1047 PetscCall((*shell->ops->multtranspose)(A, xx, y));
1048 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
1049 if (instate == outstate) {
1050 /* increase the state of the output vector since the user did not update its state themself as should have been done */
1051 PetscCall(PetscObjectStateIncrease((PetscObject)y));
1052 }
1053 PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_FALSE));
1054 PetscCall(MatShellPostScaleRight(A, y, PETSC_FALSE));
1055 PetscCall(MatShellPostZeroRight(A, y));
1056
1057 if (shell->axpy) {
1058 Mat X;
1059 PetscObjectState axpy_state;
1060
1061 PetscCall(MatShellGetContext(shell->axpy, &X));
1062 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1063 PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
1064 PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1065 PetscCall(VecCopy(x, shell->axpy_left));
1066 PetscCall(MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
1067 PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_right));
1068 }
1069 PetscFunctionReturn(PETSC_SUCCESS);
1070 }
1071
MatMultHermitianTranspose_Shell(Mat A,Vec x,Vec y)1072 static PetscErrorCode MatMultHermitianTranspose_Shell(Mat A, Vec x, Vec y)
1073 {
1074 Mat_Shell *shell = (Mat_Shell *)A->data;
1075 Vec xx;
1076 PetscObjectState instate, outstate;
1077
1078 PetscFunctionBegin;
1079 PetscCall(MatShellPreZeroLeft(A, x, &xx));
1080 PetscCall(MatShellPreScaleLeft(A, xx, &xx, PETSC_TRUE));
1081 PetscCall(PetscObjectStateGet((PetscObject)y, &instate));
1082 PetscCall((*shell->ops->multhermitiantranspose)(A, xx, y));
1083 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate));
1084 if (instate == outstate) {
1085 /* increase the state of the output vector since the user did not update its state themself as should have been done */
1086 PetscCall(PetscObjectStateIncrease((PetscObject)y));
1087 }
1088 PetscCall(MatShellShiftAndScale(A, xx, y, PETSC_TRUE));
1089 PetscCall(MatShellPostScaleRight(A, y, PETSC_TRUE));
1090 PetscCall(MatShellPostZeroRight(A, y));
1091
1092 if (shell->axpy) {
1093 Mat X;
1094 PetscObjectState axpy_state;
1095
1096 PetscCall(MatShellGetContext(shell->axpy, &X));
1097 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1098 PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
1099 PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left));
1100 PetscCall(VecCopy(x, shell->axpy_left));
1101 PetscCall(MatMultHermitianTranspose(shell->axpy, shell->axpy_left, shell->axpy_right));
1102 PetscCall(VecAXPY(y, PetscConj(shell->axpy_vscale), shell->axpy_right));
1103 }
1104 PetscFunctionReturn(PETSC_SUCCESS);
1105 }
1106
MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)1107 static PetscErrorCode MatMultTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1108 {
1109 Mat_Shell *shell = (Mat_Shell *)A->data;
1110
1111 PetscFunctionBegin;
1112 if (y == z) {
1113 if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work));
1114 PetscCall(MatMultTranspose(A, x, shell->left_add_work));
1115 PetscCall(VecAXPY(z, 1.0, shell->left_add_work));
1116 } else {
1117 PetscCall(MatMultTranspose(A, x, z));
1118 PetscCall(VecAXPY(z, 1.0, y));
1119 }
1120 PetscFunctionReturn(PETSC_SUCCESS);
1121 }
1122
MatMultHermitianTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)1123 static PetscErrorCode MatMultHermitianTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z)
1124 {
1125 Mat_Shell *shell = (Mat_Shell *)A->data;
1126
1127 PetscFunctionBegin;
1128 if (y == z) {
1129 if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work));
1130 PetscCall(MatMultHermitianTranspose(A, x, shell->left_add_work));
1131 PetscCall(VecAXPY(z, 1.0, shell->left_add_work));
1132 } else {
1133 PetscCall(MatMultHermitianTranspose(A, x, z));
1134 PetscCall(VecAXPY(z, 1.0, y));
1135 }
1136 PetscFunctionReturn(PETSC_SUCCESS);
1137 }
1138
1139 /*
1140 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1141 */
MatGetDiagonal_Shell(Mat A,Vec v)1142 static PetscErrorCode MatGetDiagonal_Shell(Mat A, Vec v)
1143 {
1144 Mat_Shell *shell = (Mat_Shell *)A->data;
1145
1146 PetscFunctionBegin;
1147 PetscCheck(shell->ops->getdiagonal, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal using MatShellSetOperation(S, MATOP_GET_DIAGONAL...)");
1148 PetscCall((*shell->ops->getdiagonal)(A, v));
1149 PetscCall(VecScale(v, shell->vscale));
1150 if (shell->dshift) PetscCall(VecAXPY(v, 1.0, shell->dshift));
1151 PetscCall(VecShift(v, shell->vshift));
1152 if (shell->left) PetscCall(VecPointwiseMult(v, v, shell->left));
1153 if (shell->right) PetscCall(VecPointwiseMult(v, v, shell->right));
1154 if (shell->zrows) {
1155 PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
1156 PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE));
1157 }
1158 if (shell->axpy) {
1159 Mat X;
1160 PetscObjectState axpy_state;
1161
1162 PetscCall(MatShellGetContext(shell->axpy, &X));
1163 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state));
1164 PetscCheck(shell->axpy_state == axpy_state, PetscObjectComm((PetscObject)A), PETSC_ERR_ORDER, "Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
1165 PetscCall(MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left));
1166 PetscCall(MatGetDiagonal(shell->axpy, shell->axpy_left));
1167 PetscCall(VecAXPY(v, shell->axpy_vscale, shell->axpy_left));
1168 }
1169 PetscFunctionReturn(PETSC_SUCCESS);
1170 }
1171
MatGetDiagonalBlock_Shell(Mat A,Mat * b)1172 static PetscErrorCode MatGetDiagonalBlock_Shell(Mat A, Mat *b)
1173 {
1174 Mat_Shell *shell = (Mat_Shell *)A->data;
1175 Vec left = NULL, right = NULL;
1176
1177 PetscFunctionBegin;
1178 PetscCheck(!shell->zrows && !shell->zcols, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME
1179 PetscCheck(!shell->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatAXPY() has been called on the input Mat"); // TODO FIXME
1180 PetscCheck(!shell->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME
1181 PetscCheck(shell->ops->getdiagonalblock, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal block using MatShellSetOperation(S, MATOP_GET_DIAGONAL_BLOCK...)");
1182 PetscCall((*shell->ops->getdiagonalblock)(A, b));
1183 PetscCall(MatScale(*b, shell->vscale));
1184 PetscCall(MatShift(*b, shell->vshift));
1185 if (shell->left) {
1186 PetscCall(VecCreateLocalVector(shell->left, &left));
1187 PetscCall(VecGetLocalVectorRead(shell->left, left));
1188 }
1189 if (shell->right) {
1190 PetscCall(VecCreateLocalVector(shell->right, &right));
1191 PetscCall(VecGetLocalVectorRead(shell->right, right));
1192 }
1193 PetscCall(MatDiagonalScale(*b, left, right));
1194 if (shell->left) {
1195 PetscCall(VecRestoreLocalVectorRead(shell->left, left));
1196 PetscCall(VecDestroy(&left));
1197 }
1198 if (shell->right) {
1199 PetscCall(VecRestoreLocalVectorRead(shell->right, right));
1200 PetscCall(VecDestroy(&right));
1201 }
1202 PetscFunctionReturn(PETSC_SUCCESS);
1203 }
1204
MatShift_Shell(Mat Y,PetscScalar a)1205 static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a)
1206 {
1207 Mat_Shell *shell = (Mat_Shell *)Y->data;
1208 PetscBool flg;
1209
1210 PetscFunctionBegin;
1211 PetscCall(MatHasCongruentLayouts(Y, &flg));
1212 PetscCheck(flg, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot shift shell matrix if it is not congruent");
1213 if (shell->left || shell->right) {
1214 if (!shell->dshift) {
1215 PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift));
1216 PetscCall(VecSet(shell->dshift, a));
1217 } else {
1218 if (shell->left) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->left));
1219 if (shell->right) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->right));
1220 PetscCall(VecShift(shell->dshift, a));
1221 }
1222 if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->left));
1223 if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->right));
1224 } else shell->vshift += a;
1225 if (shell->zrows) PetscCall(VecShift(shell->zvals, a));
1226 PetscFunctionReturn(PETSC_SUCCESS);
1227 }
1228
MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)1229 static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s)
1230 {
1231 Mat_Shell *shell = (Mat_Shell *)A->data;
1232
1233 PetscFunctionBegin;
1234 if (!shell->dshift) PetscCall(VecDuplicate(D, &shell->dshift));
1235 if (shell->left || shell->right) {
1236 if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work));
1237 if (shell->left && shell->right) {
1238 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1239 PetscCall(VecPointwiseDivide(shell->right_work, shell->right_work, shell->right));
1240 } else if (shell->left) {
1241 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left));
1242 } else {
1243 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->right));
1244 }
1245 PetscCall(VecAXPY(shell->dshift, s, shell->right_work));
1246 } else {
1247 PetscCall(VecAXPY(shell->dshift, s, D));
1248 }
1249 PetscFunctionReturn(PETSC_SUCCESS);
1250 }
1251
MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)1252 static PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins)
1253 {
1254 Mat_Shell *shell = (Mat_Shell *)A->data;
1255 Vec d;
1256 PetscBool flg;
1257
1258 PetscFunctionBegin;
1259 PetscCall(MatHasCongruentLayouts(A, &flg));
1260 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot diagonal set or shift shell matrix if it is not congruent");
1261 if (ins == INSERT_VALUES) {
1262 PetscCall(VecDuplicate(D, &d));
1263 PetscCall(MatGetDiagonal(A, d));
1264 PetscCall(MatDiagonalSet_Shell_Private(A, d, -1.));
1265 PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1266 PetscCall(VecDestroy(&d));
1267 if (shell->zrows) PetscCall(VecCopy(D, shell->zvals));
1268 } else {
1269 PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.));
1270 if (shell->zrows) PetscCall(VecAXPY(shell->zvals, 1.0, D));
1271 }
1272 PetscFunctionReturn(PETSC_SUCCESS);
1273 }
1274
MatScale_Shell(Mat Y,PetscScalar a)1275 static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a)
1276 {
1277 Mat_Shell *shell = (Mat_Shell *)Y->data;
1278
1279 PetscFunctionBegin;
1280 shell->vscale *= a;
1281 shell->vshift *= a;
1282 if (shell->dshift) PetscCall(VecScale(shell->dshift, a));
1283 shell->axpy_vscale *= a;
1284 if (shell->zrows) PetscCall(VecScale(shell->zvals, a));
1285 PetscFunctionReturn(PETSC_SUCCESS);
1286 }
1287
MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)1288 static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right)
1289 {
1290 Mat_Shell *shell = (Mat_Shell *)Y->data;
1291
1292 PetscFunctionBegin;
1293 if (left) {
1294 if (!shell->left) {
1295 PetscCall(VecDuplicate(left, &shell->left));
1296 PetscCall(VecCopy(left, shell->left));
1297 } else {
1298 PetscCall(VecPointwiseMult(shell->left, shell->left, left));
1299 }
1300 if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, left));
1301 }
1302 if (right) {
1303 if (!shell->right) {
1304 PetscCall(VecDuplicate(right, &shell->right));
1305 PetscCall(VecCopy(right, shell->right));
1306 } else {
1307 PetscCall(VecPointwiseMult(shell->right, shell->right, right));
1308 }
1309 if (shell->zrows) {
1310 if (!shell->left_work) PetscCall(MatCreateVecs(Y, NULL, &shell->left_work));
1311 PetscCall(VecSet(shell->zvals_w, 1.0));
1312 PetscCall(VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1313 PetscCall(VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD));
1314 PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w));
1315 }
1316 }
1317 if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy, left, right));
1318 PetscFunctionReturn(PETSC_SUCCESS);
1319 }
1320
MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)1321 PETSC_INTERN PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t)
1322 {
1323 Mat_Shell *shell = (Mat_Shell *)Y->data;
1324
1325 PetscFunctionBegin;
1326 if (t == MAT_FINAL_ASSEMBLY) {
1327 shell->vshift = 0.0;
1328 shell->vscale = 1.0;
1329 shell->axpy_vscale = 0.0;
1330 shell->axpy_state = 0;
1331 PetscCall(VecDestroy(&shell->dshift));
1332 PetscCall(VecDestroy(&shell->left));
1333 PetscCall(VecDestroy(&shell->right));
1334 PetscCall(MatDestroy(&shell->axpy));
1335 PetscCall(VecDestroy(&shell->axpy_left));
1336 PetscCall(VecDestroy(&shell->axpy_right));
1337 PetscCall(VecScatterDestroy(&shell->zvals_sct_c));
1338 PetscCall(VecScatterDestroy(&shell->zvals_sct_r));
1339 PetscCall(ISDestroy(&shell->zrows));
1340 PetscCall(ISDestroy(&shell->zcols));
1341 }
1342 PetscFunctionReturn(PETSC_SUCCESS);
1343 }
1344
MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)1345 static PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str)
1346 {
1347 Mat_Shell *shell = (Mat_Shell *)Y->data;
1348
1349 PetscFunctionBegin;
1350 if (X == Y) {
1351 PetscCall(MatScale(Y, 1.0 + a));
1352 PetscFunctionReturn(PETSC_SUCCESS);
1353 }
1354 if (!shell->axpy) {
1355 PetscCall(MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy));
1356 shell->axpy_vscale = a;
1357 PetscCall(PetscObjectStateGet((PetscObject)X, &shell->axpy_state));
1358 } else {
1359 PetscCall(MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str));
1360 }
1361 PetscFunctionReturn(PETSC_SUCCESS);
1362 }
1363
1364 static struct _MatOps MatOps_Values = {NULL,
1365 NULL,
1366 NULL,
1367 NULL,
1368 /* 4*/ MatMultAdd_Shell,
1369 NULL,
1370 MatMultTransposeAdd_Shell,
1371 NULL,
1372 NULL,
1373 NULL,
1374 /*10*/ NULL,
1375 NULL,
1376 NULL,
1377 NULL,
1378 NULL,
1379 /*15*/ NULL,
1380 NULL,
1381 NULL,
1382 MatDiagonalScale_Shell,
1383 NULL,
1384 /*20*/ NULL,
1385 MatAssemblyEnd_Shell,
1386 NULL,
1387 NULL,
1388 /*24*/ MatZeroRows_Shell,
1389 NULL,
1390 NULL,
1391 NULL,
1392 NULL,
1393 /*29*/ NULL,
1394 NULL,
1395 NULL,
1396 /*32*/ NULL,
1397 NULL,
1398 /*34*/ MatDuplicate_Shell,
1399 NULL,
1400 NULL,
1401 NULL,
1402 NULL,
1403 /*39*/ MatAXPY_Shell,
1404 NULL,
1405 NULL,
1406 NULL,
1407 MatCopy_Shell,
1408 /*44*/ NULL,
1409 MatScale_Shell,
1410 MatShift_Shell,
1411 MatDiagonalSet_Shell,
1412 MatZeroRowsColumns_Shell,
1413 /*49*/ NULL,
1414 NULL,
1415 NULL,
1416 NULL,
1417 NULL,
1418 /*54*/ NULL,
1419 NULL,
1420 NULL,
1421 NULL,
1422 NULL,
1423 /*59*/ NULL,
1424 MatDestroy_Shell,
1425 NULL,
1426 MatConvertFrom_Shell,
1427 NULL,
1428 /*64*/ NULL,
1429 NULL,
1430 NULL,
1431 NULL,
1432 NULL,
1433 /*69*/ NULL,
1434 MatConvert_Shell,
1435 NULL,
1436 NULL,
1437 NULL,
1438 /*74*/ NULL,
1439 NULL,
1440 NULL,
1441 NULL,
1442 NULL,
1443 /*79*/ NULL,
1444 NULL,
1445 NULL,
1446 NULL,
1447 NULL,
1448 /*84*/ NULL,
1449 NULL,
1450 NULL,
1451 NULL,
1452 NULL,
1453 /*89*/ NULL,
1454 NULL,
1455 NULL,
1456 NULL,
1457 NULL,
1458 /*94*/ NULL,
1459 NULL,
1460 NULL,
1461 NULL,
1462 NULL,
1463 /*99*/ NULL,
1464 NULL,
1465 NULL,
1466 NULL,
1467 NULL,
1468 /*104*/ NULL,
1469 NULL,
1470 NULL,
1471 NULL,
1472 NULL,
1473 /*109*/ NULL,
1474 NULL,
1475 NULL,
1476 MatMultHermitianTransposeAdd_Shell,
1477 NULL,
1478 /*114*/ NULL,
1479 NULL,
1480 NULL,
1481 NULL,
1482 NULL,
1483 /*119*/ NULL,
1484 NULL,
1485 NULL,
1486 NULL,
1487 NULL,
1488 /*124*/ NULL,
1489 NULL,
1490 NULL,
1491 NULL,
1492 NULL,
1493 /*129*/ NULL,
1494 NULL,
1495 NULL,
1496 NULL,
1497 NULL,
1498 /*134*/ NULL,
1499 NULL,
1500 NULL,
1501 NULL,
1502 NULL,
1503 /*139*/ NULL,
1504 NULL,
1505 NULL,
1506 NULL,
1507 NULL};
1508
MatShellSetContext_Shell(Mat mat,PetscCtx ctx)1509 static PetscErrorCode MatShellSetContext_Shell(Mat mat, PetscCtx ctx)
1510 {
1511 Mat_Shell *shell = (Mat_Shell *)mat->data;
1512
1513 PetscFunctionBegin;
1514 if (ctx) {
1515 PetscContainer ctxcontainer;
1516 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer));
1517 PetscCall(PetscContainerSetPointer(ctxcontainer, ctx));
1518 PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer));
1519 shell->ctxcontainer = ctxcontainer;
1520 PetscCall(PetscContainerDestroy(&ctxcontainer));
1521 } else {
1522 PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL));
1523 shell->ctxcontainer = NULL;
1524 }
1525 PetscFunctionReturn(PETSC_SUCCESS);
1526 }
1527
MatShellSetContextDestroy_Shell(Mat mat,PetscCtxDestroyFn * f)1528 static PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscCtxDestroyFn *f)
1529 {
1530 Mat_Shell *shell = (Mat_Shell *)mat->data;
1531
1532 PetscFunctionBegin;
1533 if (shell->ctxcontainer) PetscCall(PetscContainerSetCtxDestroy(shell->ctxcontainer, f));
1534 PetscFunctionReturn(PETSC_SUCCESS);
1535 }
1536
MatShellSetContext_Immutable(Mat mat,PetscCtx ctx)1537 PetscErrorCode MatShellSetContext_Immutable(Mat mat, PetscCtx ctx)
1538 {
1539 PetscFunctionBegin;
1540 SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetContext() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
1541 PetscFunctionReturn(PETSC_SUCCESS);
1542 }
1543
MatShellSetContextDestroy_Immutable(Mat mat,PetscCtxDestroyFn * f)1544 PetscErrorCode MatShellSetContextDestroy_Immutable(Mat mat, PetscCtxDestroyFn *f)
1545 {
1546 PetscFunctionBegin;
1547 SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetContextDestroy() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
1548 PetscFunctionReturn(PETSC_SUCCESS);
1549 }
1550
MatShellSetManageScalingShifts_Immutable(Mat mat)1551 PetscErrorCode MatShellSetManageScalingShifts_Immutable(Mat mat)
1552 {
1553 PetscFunctionBegin;
1554 SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "Cannot call MatShellSetManageScalingShifts() for a %s, it is used internally by the structure", ((PetscObject)mat)->type_name);
1555 PetscFunctionReturn(PETSC_SUCCESS);
1556 }
1557
MatShellSetVecType_Shell(Mat mat,VecType vtype)1558 static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype)
1559 {
1560 PetscFunctionBegin;
1561 PetscCall(PetscFree(mat->defaultvectype));
1562 PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype));
1563 PetscFunctionReturn(PETSC_SUCCESS);
1564 }
1565
MatShellSetManageScalingShifts_Shell(Mat A)1566 static PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1567 {
1568 Mat_Shell *shell = (Mat_Shell *)A->data;
1569
1570 PetscFunctionBegin;
1571 shell->managescalingshifts = PETSC_FALSE;
1572 A->ops->diagonalset = NULL;
1573 A->ops->diagonalscale = NULL;
1574 A->ops->scale = NULL;
1575 A->ops->shift = NULL;
1576 A->ops->axpy = NULL;
1577 PetscFunctionReturn(PETSC_SUCCESS);
1578 }
1579
MatShellGetScalingShifts_Shell(Mat A,PetscScalar * vshift,PetscScalar * vscale,Vec * dshift,Vec * left,Vec * right,Mat * axpy,IS * zrows,IS * zcols)1580 static PetscErrorCode MatShellGetScalingShifts_Shell(Mat A, PetscScalar *vshift, PetscScalar *vscale, Vec *dshift, Vec *left, Vec *right, Mat *axpy, IS *zrows, IS *zcols)
1581 {
1582 Mat_Shell *shell = (Mat_Shell *)A->data;
1583
1584 PetscFunctionBegin;
1585 PetscCheck(shell->managescalingshifts, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot get scaling and shifts if MatShellSetManageScalingShifts() has been called");
1586 if (vshift == MAT_SHELL_NOT_ALLOWED) PetscCheck(shell->vshift == 0.0, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: vshift != 0.0, set via MatShift()");
1587 else if (vshift) *vshift = shell->vshift;
1588 if (vscale == MAT_SHELL_NOT_ALLOWED) PetscCheck(shell->vscale == 1.0, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: vscale != 1.0, set via MatScale()");
1589 else if (vscale) *vscale = shell->vscale;
1590 if (dshift == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: dshift, set via MatDiagonalSet()");
1591 else if (dshift) *dshift = shell->dshift;
1592 if (left == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->left, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: left, set via MatDiagonalScale()");
1593 else if (left) *left = shell->left;
1594 if (right == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: right, set via MatDiagonalScale()");
1595 else if (right) *right = shell->right;
1596 if (axpy == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: axpy, set via MatAXPY()");
1597 else if (axpy) *axpy = shell->axpy;
1598 if (zrows == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->zrows, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: zrows, set via MatZeroRows()");
1599 else if (zrows) *zrows = shell->zrows;
1600 if (zcols == MAT_SHELL_NOT_ALLOWED) PetscCheck(!shell->zcols, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Non-trivial member in parent MatShell: zcols, set via MatZeroRowsColumns()");
1601 else if (zcols) *zcols = shell->zcols;
1602 PetscFunctionReturn(PETSC_SUCCESS);
1603 }
1604
MatShellSetOperation_Shell(Mat mat,MatOperation op,PetscErrorCodeFn * f)1605 static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, PetscErrorCodeFn *f)
1606 {
1607 Mat_Shell *shell = (Mat_Shell *)mat->data;
1608
1609 PetscFunctionBegin;
1610 switch (op) {
1611 case MATOP_DESTROY:
1612 shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1613 break;
1614 case MATOP_VIEW:
1615 if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view;
1616 mat->ops->view = (PetscErrorCode (*)(Mat, PetscViewer))f;
1617 break;
1618 case MATOP_COPY:
1619 shell->ops->copy = (PetscErrorCode (*)(Mat, Mat, MatStructure))f;
1620 break;
1621 case MATOP_DIAGONAL_SET:
1622 case MATOP_DIAGONAL_SCALE:
1623 case MATOP_SHIFT:
1624 case MATOP_SCALE:
1625 case MATOP_AXPY:
1626 case MATOP_ZERO_ROWS:
1627 case MATOP_ZERO_ROWS_COLUMNS:
1628 case MATOP_ZERO_ROWS_LOCAL:
1629 case MATOP_ZERO_ROWS_COLUMNS_LOCAL:
1630 PetscCheck(!shell->managescalingshifts, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1631 (((PetscErrorCodeFn **)mat->ops)[op]) = f;
1632 break;
1633 case MATOP_GET_DIAGONAL:
1634 if (shell->managescalingshifts) {
1635 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat, Vec))f;
1636 mat->ops->getdiagonal = MatGetDiagonal_Shell;
1637 } else {
1638 shell->ops->getdiagonal = NULL;
1639 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat, Vec))f;
1640 }
1641 break;
1642 case MATOP_GET_DIAGONAL_BLOCK:
1643 if (shell->managescalingshifts) {
1644 shell->ops->getdiagonalblock = (PetscErrorCode (*)(Mat, Mat *))f;
1645 mat->ops->getdiagonalblock = MatGetDiagonalBlock_Shell;
1646 } else {
1647 shell->ops->getdiagonalblock = NULL;
1648 mat->ops->getdiagonalblock = (PetscErrorCode (*)(Mat, Mat *))f;
1649 }
1650 break;
1651 case MATOP_MULT:
1652 if (shell->managescalingshifts) {
1653 shell->ops->mult = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1654 mat->ops->mult = MatMult_Shell;
1655 } else {
1656 shell->ops->mult = NULL;
1657 mat->ops->mult = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1658 }
1659 break;
1660 case MATOP_MULT_TRANSPOSE:
1661 if (shell->managescalingshifts) {
1662 shell->ops->multtranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1663 mat->ops->multtranspose = MatMultTranspose_Shell;
1664 } else {
1665 shell->ops->multtranspose = NULL;
1666 mat->ops->multtranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1667 }
1668 break;
1669 case MATOP_MULT_HERMITIAN_TRANSPOSE:
1670 if (shell->managescalingshifts) {
1671 shell->ops->multhermitiantranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1672 mat->ops->multhermitiantranspose = MatMultHermitianTranspose_Shell;
1673 } else {
1674 shell->ops->multhermitiantranspose = NULL;
1675 mat->ops->multhermitiantranspose = (PetscErrorCode (*)(Mat, Vec, Vec))f;
1676 }
1677 break;
1678 default:
1679 (((PetscErrorCodeFn **)mat->ops)[op]) = f;
1680 break;
1681 }
1682 PetscFunctionReturn(PETSC_SUCCESS);
1683 }
1684
MatShellGetOperation_Shell(Mat mat,MatOperation op,PetscErrorCodeFn ** f)1685 static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, PetscErrorCodeFn **f)
1686 {
1687 Mat_Shell *shell = (Mat_Shell *)mat->data;
1688
1689 PetscFunctionBegin;
1690 switch (op) {
1691 case MATOP_DESTROY:
1692 *f = (PetscErrorCodeFn *)shell->ops->destroy;
1693 break;
1694 case MATOP_VIEW:
1695 *f = (PetscErrorCodeFn *)mat->ops->view;
1696 break;
1697 case MATOP_COPY:
1698 *f = (PetscErrorCodeFn *)shell->ops->copy;
1699 break;
1700 case MATOP_DIAGONAL_SET:
1701 case MATOP_DIAGONAL_SCALE:
1702 case MATOP_SHIFT:
1703 case MATOP_SCALE:
1704 case MATOP_AXPY:
1705 case MATOP_ZERO_ROWS:
1706 case MATOP_ZERO_ROWS_COLUMNS:
1707 *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1708 break;
1709 case MATOP_GET_DIAGONAL:
1710 if (shell->ops->getdiagonal) *f = (PetscErrorCodeFn *)shell->ops->getdiagonal;
1711 else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1712 break;
1713 case MATOP_GET_DIAGONAL_BLOCK:
1714 if (shell->ops->getdiagonalblock) *f = (PetscErrorCodeFn *)shell->ops->getdiagonalblock;
1715 else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1716 break;
1717 case MATOP_MULT:
1718 if (shell->ops->mult) *f = (PetscErrorCodeFn *)shell->ops->mult;
1719 else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1720 break;
1721 case MATOP_MULT_TRANSPOSE:
1722 if (shell->ops->multtranspose) *f = (PetscErrorCodeFn *)shell->ops->multtranspose;
1723 else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1724 break;
1725 case MATOP_MULT_HERMITIAN_TRANSPOSE:
1726 if (shell->ops->multhermitiantranspose) *f = (PetscErrorCodeFn *)shell->ops->multhermitiantranspose;
1727 else *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1728 break;
1729 default:
1730 *f = (((PetscErrorCodeFn **)mat->ops)[op]);
1731 }
1732 PetscFunctionReturn(PETSC_SUCCESS);
1733 }
1734
1735 /*MC
1736 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type with its own data structure -- perhaps matrix-free.
1737
1738 Level: advanced
1739
1740 Notes:
1741 See `MatCreateShell()` for details on the usage of `MATSHELL`
1742
1743 `PCSHELL` can be used in conjunction with `MATSHELL` to provide a custom preconditioner appropriate for your `MATSHELL`. Since
1744 many standard preconditioners such as `PCILU` depend on having an explicit representation of the matrix entries they cannot be used
1745 directly with `MATSHELL`.
1746
1747 .seealso: [](ch_matrices), `Mat`, `MatCreateShell()`, `PCSHELL`
1748 M*/
1749
MatCreate_Shell(Mat A)1750 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1751 {
1752 Mat_Shell *b;
1753
1754 PetscFunctionBegin;
1755 PetscCall(PetscNew(&b));
1756 A->data = (void *)b;
1757 A->ops[0] = MatOps_Values;
1758
1759 b->ctxcontainer = NULL;
1760 b->vshift = 0.0;
1761 b->vscale = 1.0;
1762 b->managescalingshifts = PETSC_TRUE;
1763 A->assembled = PETSC_TRUE;
1764 A->preallocated = PETSC_FALSE;
1765
1766 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell));
1767 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell));
1768 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell));
1769 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell));
1770 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell));
1771 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetScalingShifts_C", MatShellGetScalingShifts_Shell));
1772 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell));
1773 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell));
1774 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell));
1775 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL));
1776 PetscFunctionReturn(PETSC_SUCCESS);
1777 }
1778
1779 /*@C
1780 MatCreateShell - Creates a new matrix of `MatType` `MATSHELL` for use with a user-defined
1781 private matrix data storage format.
1782
1783 Collective
1784
1785 Input Parameters:
1786 + comm - MPI communicator
1787 . m - number of local rows (or `PETSC_DECIDE` to have calculated if `M` is given)
1788 . n - number of local columns (or `PETSC_DECIDE` to have calculated if `N` is given)
1789 . M - number of global rows (may be `PETSC_DETERMINE` to have calculated if `m` is given)
1790 . N - number of global columns (may be `PETSC_DETERMINE` to have calculated if `n` is given)
1791 - ctx - pointer to data needed by the shell matrix routines
1792
1793 Output Parameter:
1794 . A - the matrix
1795
1796 Level: advanced
1797
1798 Example Usage:
1799 .vb
1800 extern PetscErrorCode mult(Mat, Vec, Vec);
1801
1802 MatCreateShell(comm, m, n, M, N, ctx, &mat);
1803 MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))mult);
1804 MatShellSetContext(mat,ctx);
1805 // Use matrix for operations that have been set
1806 MatDestroy(mat);
1807 .ve
1808
1809 Notes:
1810 The shell matrix type is intended to provide a simple way for users to write a custom matrix specifically for their application.
1811
1812 `MatCreateShell()` is used in conjunction with `MatShellSetContext()` and `MatShellSetOperation()`.
1813
1814 PETSc requires that matrices and vectors being used for certain
1815 operations are partitioned accordingly. For example, when
1816 creating a shell matrix, `A`, that supports parallel matrix-vector
1817 products using `MatMult`(A,x,y) the user should set the number
1818 of local matrix rows to be the number of local elements of the
1819 corresponding result vector, y. Note that this is information is
1820 required for use of the matrix interface routines, even though
1821 the shell matrix may not actually be physically partitioned.
1822 For example,
1823
1824 .vb
1825 Vec x, y
1826 extern PetscErrorCode mult(Mat,Vec,Vec);
1827 Mat A
1828
1829 VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1830 VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1831 VecGetLocalSize(y,&m);
1832 VecGetLocalSize(x,&n);
1833 MatCreateShell(comm,m,n,M,N,ctx,&A);
1834 MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1835 MatMult(A,x,y);
1836 MatDestroy(&A);
1837 VecDestroy(&y);
1838 VecDestroy(&x);
1839 .ve
1840
1841 `MATSHELL` handles `MatShift()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatScale()`, `MatZeroRows()` and `MatZeroRowsColumns()`
1842 internally, so these operations cannot be overwritten unless `MatShellSetManageScalingShifts()` is called.
1843
1844 Developer Notes:
1845 For rectangular matrices do all the scalings and shifts make sense?
1846
1847 Regarding shifting and scaling. The general form is
1848
1849 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1850
1851 The order you apply the operations is important. For example if you have a dshift then
1852 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
1853 you get s*vscale*A + diag(shift)
1854
1855 A is the user provided function.
1856
1857 `KSP`/`PC` uses changes in the `Mat`'s "state" to decide if preconditioners need to be rebuilt `PCSetUp()` only calls the setup() for
1858 for the `PC` implementation if the `Mat` state has increased from the previous call. Thus to get changes in a `MATSHELL` to trigger
1859 an update in the preconditioner you must call `MatAssemblyBegin()` and `MatAssemblyEnd()` or `PetscObjectStateIncrease`((`PetscObject`)mat);
1860 each time the `MATSHELL` matrix has changed.
1861
1862 Matrix-matrix product operations can be specified using `MatShellSetMatProductOperation()`
1863
1864 Calling `MatAssemblyBegin()`/`MatAssemblyEnd()` on a `MATSHELL` removes any previously supplied shift and scales that were provided
1865 with `MatDiagonalSet()`, `MatShift()`, `MatScale()`, or `MatDiagonalScale()`.
1866
1867 Fortran Notes:
1868 To use this from Fortran with a `ctx` you must write an interface definition for this
1869 function and for `MatShellGetContext()` that tells Fortran the Fortran derived data type you are passing
1870 in as the `ctx` argument.
1871
1872 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
1873 @*/
MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscCtx ctx,Mat * A)1874 PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscCtx ctx, Mat *A)
1875 {
1876 PetscFunctionBegin;
1877 PetscCall(MatCreate(comm, A));
1878 PetscCall(MatSetSizes(*A, m, n, M, N));
1879 PetscCall(MatSetType(*A, MATSHELL));
1880 PetscCall(MatShellSetContext(*A, ctx));
1881 PetscCall(MatSetUp(*A));
1882 PetscFunctionReturn(PETSC_SUCCESS);
1883 }
1884
1885 /*@
1886 MatShellSetContext - sets the context for a `MATSHELL` shell matrix
1887
1888 Logically Collective
1889
1890 Input Parameters:
1891 + mat - the `MATSHELL` shell matrix
1892 - ctx - the context
1893
1894 Level: advanced
1895
1896 Note:
1897 This provides an easy way, along with `MatCreateShell()` and `MatShellSetOperation()` to provide a custom matrix format
1898 specifically for your application.
1899
1900 Fortran Notes:
1901 You must write a Fortran interface definition for this
1902 function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument.
1903
1904 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`
1905 @*/
MatShellSetContext(Mat mat,PetscCtx ctx)1906 PetscErrorCode MatShellSetContext(Mat mat, PetscCtx ctx)
1907 {
1908 PetscFunctionBegin;
1909 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
1910 PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx));
1911 PetscFunctionReturn(PETSC_SUCCESS);
1912 }
1913
1914 /*@C
1915 MatShellSetContextDestroy - sets the destroy function for a `MATSHELL` shell matrix context
1916
1917 Logically Collective
1918
1919 Input Parameters:
1920 + mat - the shell matrix
1921 - f - the context destroy function, see `PetscCtxDestroyFn` for calling sequence
1922
1923 Level: advanced
1924
1925 Note:
1926 If the `MatShell` is never duplicated, the behavior of this function is equivalent
1927 to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()`
1928 ensures proper reference counting for the user provided context data in the case that
1929 the `MATSHELL` is duplicated.
1930
1931 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()`,
1932 `PetscCtxDestroyFn`
1933 @*/
MatShellSetContextDestroy(Mat mat,PetscCtxDestroyFn * f)1934 PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscCtxDestroyFn *f)
1935 {
1936 PetscFunctionBegin;
1937 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
1938 PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscCtxDestroyFn *), (mat, f));
1939 PetscFunctionReturn(PETSC_SUCCESS);
1940 }
1941
1942 /*@C
1943 MatShellSetVecType - Sets the `VecType` of `Vec` returned by `MatCreateVecs()`
1944
1945 Logically Collective
1946
1947 Input Parameters:
1948 + mat - the `MATSHELL` shell matrix
1949 - vtype - type to use for creating vectors
1950
1951 Level: advanced
1952
1953 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateVecs()`
1954 @*/
MatShellSetVecType(Mat mat,VecType vtype)1955 PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype)
1956 {
1957 PetscFunctionBegin;
1958 PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype));
1959 PetscFunctionReturn(PETSC_SUCCESS);
1960 }
1961
1962 /*@
1963 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the `MATSHELL`. Must be called immediately
1964 after `MatCreateShell()`
1965
1966 Logically Collective
1967
1968 Input Parameter:
1969 . A - the `MATSHELL` shell matrix
1970
1971 Level: advanced
1972
1973 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()`
1974 @*/
MatShellSetManageScalingShifts(Mat A)1975 PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1976 {
1977 PetscFunctionBegin;
1978 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
1979 PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A));
1980 PetscFunctionReturn(PETSC_SUCCESS);
1981 }
1982
1983 /*@C
1984 MatShellGetScalingShifts - Gets members of a `MATSHELL` used internally for scaling and
1985 shifting the `Mat` or calling `MatAXPY()`, `MatZeroRows()`, or `MatZeroRowsColumns()` with it
1986
1987 Logically Collective
1988
1989 Input Parameter:
1990 . A - the `MATSHELL` shell matrix
1991
1992 Output Parameters:
1993 + vshift - `PetscScalar` pointer (can be `NULL`), see `MatShift()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be 0
1994 . vscale - `PetscScalar` pointer (can be `NULL`), see `MatScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be 1
1995 . dshift - `Vec` pointer (can be `NULL`), see `MatDiagonalSet()`, or `MAT_SHELL_NOT_ALLOWED` if the internal shift should be `NULL`
1996 . left - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
1997 . right - `Vec` pointer (can be `NULL`), see `MatDiagonalScale()`, or `MAT_SHELL_NOT_ALLOWED` if the internal scaling should be `NULL`
1998 . axpy - `Mat` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatAXPY()` should have not been called on `A`
1999 . zrows - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRows()` should have not been called on `A`
2000 - zcols - `Vec` pointer (can be `NULL`), or `MAT_SHELL_NOT_ALLOWED` if `MatZeroRowsColumns()` should have not been called on `A`
2001
2002 Level: advanced
2003
2004 Developer Notes:
2005 This is mostly useful to check for corner-cases in `MatType` deriving from
2006 `MATSHELL`, e.g, `MATCOMPOSITE` or `MATTRANSPOSEVIRTUAL`, since scaling and
2007 shifts often require extra work which is not always implemented.
2008
2009 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShift()`, `MatScale()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatZeroRows()`, `MatZeroRowsColumns()`, `MatShellSetManageScalingShifts()`
2010 @*/
MatShellGetScalingShifts(Mat A,PetscScalar * vshift,PetscScalar * vscale,Vec * dshift,Vec * left,Vec * right,Mat * axpy,IS * zrows,IS * zcols)2011 PetscErrorCode MatShellGetScalingShifts(Mat A, PetscScalar *vshift, PetscScalar *vscale, Vec *dshift, Vec *left, Vec *right, Mat *axpy, IS *zrows, IS *zcols)
2012 {
2013 PetscFunctionBegin;
2014 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
2015 PetscTryMethod(A, "MatShellGetScalingShifts_C", (Mat, PetscScalar *, PetscScalar *, Vec *, Vec *, Vec *, Mat *, IS *, IS *), (A, vshift, vscale, dshift, left, right, axpy, zrows, zcols));
2016 PetscFunctionReturn(PETSC_SUCCESS);
2017 }
2018
2019 /*@C
2020 MatShellTestMult - Compares the multiply routine provided to the `MATSHELL` with differencing on a given function.
2021
2022 Logically Collective; No Fortran Support
2023
2024 Input Parameters:
2025 + mat - the `MATSHELL` shell matrix
2026 . f - the function
2027 . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
2028 - ctx - an optional context for the function
2029
2030 Output Parameter:
2031 . flg - `PETSC_TRUE` if the multiply is likely correct
2032
2033 Options Database Key:
2034 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
2035
2036 Level: advanced
2037
2038 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()`
2039 @*/
MatShellTestMult(Mat mat,PetscErrorCode (* f)(void *,Vec,Vec),Vec base,PetscCtx ctx,PetscBool * flg)2040 PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, PetscCtx ctx, PetscBool *flg)
2041 {
2042 PetscInt m, n;
2043 Mat mf, Dmf, Dmat, Ddiff;
2044 PetscReal Diffnorm, Dmfnorm;
2045 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
2046
2047 PetscFunctionBegin;
2048 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
2049 PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v));
2050 PetscCall(MatGetLocalSize(mat, &m, &n));
2051 PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf));
2052 PetscCall(MatMFFDSetFunction(mf, f, ctx));
2053 PetscCall(MatMFFDSetBase(mf, base, NULL));
2054
2055 PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
2056 PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat));
2057
2058 PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
2059 PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
2060 PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
2061 PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
2062 if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
2063 flag = PETSC_FALSE;
2064 if (v) {
2065 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n", (double)(Diffnorm / Dmfnorm)));
2066 PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view"));
2067 PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view"));
2068 PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view"));
2069 }
2070 } else if (v) {
2071 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce the same results\n"));
2072 }
2073 if (flg) *flg = flag;
2074 PetscCall(MatDestroy(&Ddiff));
2075 PetscCall(MatDestroy(&mf));
2076 PetscCall(MatDestroy(&Dmf));
2077 PetscCall(MatDestroy(&Dmat));
2078 PetscFunctionReturn(PETSC_SUCCESS);
2079 }
2080
2081 /*@C
2082 MatShellTestMultTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on a given function.
2083
2084 Logically Collective; No Fortran Support
2085
2086 Input Parameters:
2087 + mat - the `MATSHELL` shell matrix
2088 . f - the function
2089 . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated
2090 - ctx - an optional context for the function
2091
2092 Output Parameter:
2093 . flg - `PETSC_TRUE` if the multiply is likely correct
2094
2095 Options Database Key:
2096 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
2097
2098 Level: advanced
2099
2100 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()`
2101 @*/
MatShellTestMultTranspose(Mat mat,PetscErrorCode (* f)(void *,Vec,Vec),Vec base,PetscCtx ctx,PetscBool * flg)2102 PetscErrorCode MatShellTestMultTranspose(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, PetscCtx ctx, PetscBool *flg)
2103 {
2104 Vec x, y, z;
2105 PetscInt m, n, M, N;
2106 Mat mf, Dmf, Dmat, Ddiff;
2107 PetscReal Diffnorm, Dmfnorm;
2108 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
2109
2110 PetscFunctionBegin;
2111 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
2112 PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_transpose_view", &v));
2113 PetscCall(MatCreateVecs(mat, &x, &y));
2114 PetscCall(VecDuplicate(y, &z));
2115 PetscCall(MatGetLocalSize(mat, &m, &n));
2116 PetscCall(MatGetSize(mat, &M, &N));
2117 PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, M, N, &mf));
2118 PetscCall(MatMFFDSetFunction(mf, f, ctx));
2119 PetscCall(MatMFFDSetBase(mf, base, NULL));
2120 PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf));
2121 PetscCall(MatTranspose(Dmf, MAT_INPLACE_MATRIX, &Dmf));
2122 PetscCall(MatComputeOperatorTranspose(mat, MATAIJ, &Dmat));
2123
2124 PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff));
2125 PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN));
2126 PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm));
2127 PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm));
2128 if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) {
2129 flag = PETSC_FALSE;
2130 if (v) {
2131 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n", (double)(Diffnorm / Dmfnorm)));
2132 PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2133 PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2134 PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_transpose_view"));
2135 }
2136 } else if (v) {
2137 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL transpose and matrix-free multiple appear to produce the same results\n"));
2138 }
2139 if (flg) *flg = flag;
2140 PetscCall(MatDestroy(&mf));
2141 PetscCall(MatDestroy(&Dmat));
2142 PetscCall(MatDestroy(&Ddiff));
2143 PetscCall(MatDestroy(&Dmf));
2144 PetscCall(VecDestroy(&x));
2145 PetscCall(VecDestroy(&y));
2146 PetscCall(VecDestroy(&z));
2147 PetscFunctionReturn(PETSC_SUCCESS);
2148 }
2149
2150 /*@C
2151 MatShellSetOperation - Allows user to set a matrix operation for a `MATSHELL` shell matrix.
2152
2153 Logically Collective
2154
2155 Input Parameters:
2156 + mat - the `MATSHELL` shell matrix
2157 . op - the name of the operation
2158 - g - the function that provides the operation.
2159
2160 Level: advanced
2161
2162 Example Usage:
2163 .vb
2164 extern PetscErrorCode usermult(Mat, Vec, Vec);
2165
2166 MatCreateShell(comm, m, n, M, N, ctx, &A);
2167 MatShellSetOperation(A, MATOP_MULT, (void(*)(void))usermult);
2168 .ve
2169
2170 Notes:
2171 See the file include/petscmat.h for a complete list of matrix
2172 operations, which all have the form MATOP_<OPERATION>, where
2173 <OPERATION> is the name (in all capital letters) of the
2174 user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).
2175
2176 All user-provided functions (except for `MATOP_DESTROY`) should have the same calling
2177 sequence as the usual matrix interface routines, since they
2178 are intended to be accessed via the usual matrix interface
2179 routines, e.g.,
2180 .vb
2181 MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)
2182 .ve
2183
2184 In particular each function MUST return an error code of `PETSC_SUCCESS` on success and
2185 nonzero on failure.
2186
2187 Within each user-defined routine, the user should call
2188 `MatShellGetContext()` to obtain the user-defined context that was
2189 set by `MatCreateShell()`.
2190
2191 Use `MatSetOperation()` to set an operation for any matrix type. For matrix product operations (i.e. `MatMatXXX()`, `MatTransposeMatXXX()` etc)
2192 use `MatShellSetMatProductOperation()`
2193
2194 Fortran Note:
2195 For `MatCreateVecs()` the user code should check if the input left or right vector is -1 and in that case not
2196 generate a vector. See `src/mat/tests/ex120f.F`
2197
2198 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()`
2199 @*/
MatShellSetOperation(Mat mat,MatOperation op,PetscErrorCodeFn * g)2200 PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, PetscErrorCodeFn *g)
2201 {
2202 PetscFunctionBegin;
2203 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
2204 PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, PetscErrorCodeFn *), (mat, op, g));
2205 PetscFunctionReturn(PETSC_SUCCESS);
2206 }
2207
2208 /*@C
2209 MatShellGetOperation - Gets a matrix function for a `MATSHELL` shell matrix.
2210
2211 Not Collective
2212
2213 Input Parameters:
2214 + mat - the `MATSHELL` shell matrix
2215 - op - the name of the operation
2216
2217 Output Parameter:
2218 . g - the function that provides the operation.
2219
2220 Level: advanced
2221
2222 Notes:
2223 See the file include/petscmat.h for a complete list of matrix
2224 operations, which all have the form MATOP_<OPERATION>, where
2225 <OPERATION> is the name (in all capital letters) of the
2226 user interface routine (e.g., `MatMult()` -> `MATOP_MULT`).
2227
2228 All user-provided functions have the same calling
2229 sequence as the usual matrix interface routines, since they
2230 are intended to be accessed via the usual matrix interface
2231 routines, e.g.,
2232 .vb
2233 MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec)
2234 .ve
2235
2236 Within each user-defined routine, the user should call
2237 `MatShellGetContext()` to obtain the user-defined context that was
2238 set by `MatCreateShell()`.
2239
2240 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()`
2241 @*/
MatShellGetOperation(Mat mat,MatOperation op,PetscErrorCodeFn ** g)2242 PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, PetscErrorCodeFn **g)
2243 {
2244 PetscFunctionBegin;
2245 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
2246 PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, PetscErrorCodeFn **), (mat, op, g));
2247 PetscFunctionReturn(PETSC_SUCCESS);
2248 }
2249
2250 /*@
2251 MatIsShell - Inquires if a matrix is derived from `MATSHELL`
2252
2253 Input Parameter:
2254 . mat - the matrix
2255
2256 Output Parameter:
2257 . flg - the Boolean value
2258
2259 Level: developer
2260
2261 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT`
2262 @*/
MatIsShell(Mat mat,PetscBool * flg)2263 PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2264 {
2265 PetscFunctionBegin;
2266 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1);
2267 PetscAssertPointer(flg, 2);
2268 *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2269 PetscFunctionReturn(PETSC_SUCCESS);
2270 }
2271