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