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