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