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