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 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 Notes: 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: `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 shell matrix. 813 814 Logically Collective on Mat 815 816 Input Parameters: 817 + A - the 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 If the symbolic phase is not specified, MatSetUp() is called on the result matrix that must have its type set if Ctype is NULL. 844 Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback. 845 PETSc will take care of calling the user-defined callbacks. 846 It is allowed to specify the same callbacks for different Btype matrix types. 847 The couple (Btype,ptype) uniquely identifies the operation: the last specified callbacks takes precedence. 848 849 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatProductType`, `MatType`, `MatSetUp()` 850 @*/ 851 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) { 852 PetscFunctionBegin; 853 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 854 PetscValidLogicalCollectiveEnum(A, ptype, 2); 855 PetscCheck(ptype != MATPRODUCT_ABC, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not for product type %s", MatProductTypes[ptype]); 856 PetscCheck(numeric, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "Missing numeric routine, argument 4"); 857 PetscValidCharPointer(Btype, 6); 858 if (Ctype) PetscValidCharPointer(Ctype, 7); 859 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)); 860 PetscFunctionReturn(0); 861 } 862 863 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) { 864 PetscBool flg; 865 char composedname[256]; 866 MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList; 867 PetscMPIInt size; 868 869 PetscFunctionBegin; 870 PetscValidType(A, 1); 871 while (Bnames) { /* user passed in the root name */ 872 PetscCall(PetscStrcmp(Btype, Bnames->rname, &flg)); 873 if (flg) break; 874 Bnames = Bnames->next; 875 } 876 while (Cnames) { /* user passed in the root name */ 877 PetscCall(PetscStrcmp(Ctype, Cnames->rname, &flg)); 878 if (flg) break; 879 Cnames = Cnames->next; 880 } 881 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size)); 882 Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype; 883 Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype; 884 PetscCall(PetscSNPrintf(composedname, sizeof(composedname), "MatProductSetFromOptions_%s_%s_C", ((PetscObject)A)->type_name, Btype)); 885 PetscCall(MatShellSetMatProductOperation_Private(A, ptype, symbolic, numeric, destroy, composedname, Ctype)); 886 PetscFunctionReturn(0); 887 } 888 889 static PetscErrorCode MatCopy_Shell(Mat A, Mat B, MatStructure str) { 890 Mat_Shell *shellA = (Mat_Shell *)A->data, *shellB = (Mat_Shell *)B->data; 891 PetscBool matflg; 892 MatShellMatFunctionList matmatA; 893 894 PetscFunctionBegin; 895 PetscCall(MatIsShell(B, &matflg)); 896 PetscCheck(matflg, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix %s not derived from MATSHELL", ((PetscObject)B)->type_name); 897 898 PetscCall(PetscMemcpy(B->ops, A->ops, sizeof(struct _MatOps))); 899 PetscCall(PetscMemcpy(shellB->ops, shellA->ops, sizeof(struct _MatShellOps))); 900 901 if (shellA->ops->copy) PetscCall((*shellA->ops->copy)(A, B, str)); 902 shellB->vscale = shellA->vscale; 903 shellB->vshift = shellA->vshift; 904 if (shellA->dshift) { 905 if (!shellB->dshift) PetscCall(VecDuplicate(shellA->dshift, &shellB->dshift)); 906 PetscCall(VecCopy(shellA->dshift, shellB->dshift)); 907 } else { 908 PetscCall(VecDestroy(&shellB->dshift)); 909 } 910 if (shellA->left) { 911 if (!shellB->left) PetscCall(VecDuplicate(shellA->left, &shellB->left)); 912 PetscCall(VecCopy(shellA->left, shellB->left)); 913 } else { 914 PetscCall(VecDestroy(&shellB->left)); 915 } 916 if (shellA->right) { 917 if (!shellB->right) PetscCall(VecDuplicate(shellA->right, &shellB->right)); 918 PetscCall(VecCopy(shellA->right, shellB->right)); 919 } else { 920 PetscCall(VecDestroy(&shellB->right)); 921 } 922 PetscCall(MatDestroy(&shellB->axpy)); 923 shellB->axpy_vscale = 0.0; 924 shellB->axpy_state = 0; 925 if (shellA->axpy) { 926 PetscCall(PetscObjectReference((PetscObject)shellA->axpy)); 927 shellB->axpy = shellA->axpy; 928 shellB->axpy_vscale = shellA->axpy_vscale; 929 shellB->axpy_state = shellA->axpy_state; 930 } 931 if (shellA->zrows) { 932 PetscCall(ISDuplicate(shellA->zrows, &shellB->zrows)); 933 if (shellA->zcols) PetscCall(ISDuplicate(shellA->zcols, &shellB->zcols)); 934 PetscCall(VecDuplicate(shellA->zvals, &shellB->zvals)); 935 PetscCall(VecCopy(shellA->zvals, shellB->zvals)); 936 PetscCall(VecDuplicate(shellA->zvals_w, &shellB->zvals_w)); 937 PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_r)); 938 PetscCall(PetscObjectReference((PetscObject)shellA->zvals_sct_c)); 939 shellB->zvals_sct_r = shellA->zvals_sct_r; 940 shellB->zvals_sct_c = shellA->zvals_sct_c; 941 } 942 943 matmatA = shellA->matmat; 944 if (matmatA) { 945 while (matmatA->next) { 946 PetscCall(MatShellSetMatProductOperation_Private(B, matmatA->ptype, matmatA->symbolic, matmatA->numeric, matmatA->destroy, matmatA->composedname, matmatA->resultname)); 947 matmatA = matmatA->next; 948 } 949 } 950 PetscFunctionReturn(0); 951 } 952 953 static PetscErrorCode MatDuplicate_Shell(Mat mat, MatDuplicateOption op, Mat *M) { 954 PetscFunctionBegin; 955 PetscCall(MatCreateShell(PetscObjectComm((PetscObject)mat), mat->rmap->n, mat->cmap->n, mat->rmap->N, mat->cmap->N, NULL, M)); 956 ((Mat_Shell *)(*M)->data)->ctxcontainer = ((Mat_Shell *)mat->data)->ctxcontainer; 957 PetscCall(PetscObjectCompose((PetscObject)(*M), "MatShell ctx", (PetscObject)((Mat_Shell *)(*M)->data)->ctxcontainer)); 958 PetscCall(PetscObjectChangeTypeName((PetscObject)(*M), ((PetscObject)mat)->type_name)); 959 if (op != MAT_DO_NOT_COPY_VALUES) PetscCall(MatCopy(mat, *M, SAME_NONZERO_PATTERN)); 960 PetscFunctionReturn(0); 961 } 962 963 PetscErrorCode MatMult_Shell(Mat A, Vec x, Vec y) { 964 Mat_Shell *shell = (Mat_Shell *)A->data; 965 Vec xx; 966 PetscObjectState instate, outstate; 967 968 PetscFunctionBegin; 969 PetscCall(MatShellPreZeroRight(A, x, &xx)); 970 PetscCall(MatShellPreScaleRight(A, xx, &xx)); 971 PetscCall(PetscObjectStateGet((PetscObject)y, &instate)); 972 PetscCall((*shell->ops->mult)(A, xx, y)); 973 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate)); 974 if (instate == outstate) { 975 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 976 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 977 } 978 PetscCall(MatShellShiftAndScale(A, xx, y)); 979 PetscCall(MatShellPostScaleLeft(A, y)); 980 PetscCall(MatShellPostZeroLeft(A, y)); 981 982 if (shell->axpy) { 983 Mat X; 984 PetscObjectState axpy_state; 985 986 PetscCall(MatShellGetContext(shell->axpy, &X)); 987 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state)); 988 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,...)"); 989 990 PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left)); 991 PetscCall(VecCopy(x, shell->axpy_right)); 992 PetscCall(MatMult(shell->axpy, shell->axpy_right, shell->axpy_left)); 993 PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_left)); 994 } 995 PetscFunctionReturn(0); 996 } 997 998 static PetscErrorCode MatMultAdd_Shell(Mat A, Vec x, Vec y, Vec z) { 999 Mat_Shell *shell = (Mat_Shell *)A->data; 1000 1001 PetscFunctionBegin; 1002 if (y == z) { 1003 if (!shell->right_add_work) PetscCall(VecDuplicate(z, &shell->right_add_work)); 1004 PetscCall(MatMult(A, x, shell->right_add_work)); 1005 PetscCall(VecAXPY(z, 1.0, shell->right_add_work)); 1006 } else { 1007 PetscCall(MatMult(A, x, z)); 1008 PetscCall(VecAXPY(z, 1.0, y)); 1009 } 1010 PetscFunctionReturn(0); 1011 } 1012 1013 static PetscErrorCode MatMultTranspose_Shell(Mat A, Vec x, Vec y) { 1014 Mat_Shell *shell = (Mat_Shell *)A->data; 1015 Vec xx; 1016 PetscObjectState instate, outstate; 1017 1018 PetscFunctionBegin; 1019 PetscCall(MatShellPreZeroLeft(A, x, &xx)); 1020 PetscCall(MatShellPreScaleLeft(A, xx, &xx)); 1021 PetscCall(PetscObjectStateGet((PetscObject)y, &instate)); 1022 PetscCall((*shell->ops->multtranspose)(A, xx, y)); 1023 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate)); 1024 if (instate == outstate) { 1025 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1026 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 1027 } 1028 PetscCall(MatShellShiftAndScale(A, xx, y)); 1029 PetscCall(MatShellPostScaleRight(A, y)); 1030 PetscCall(MatShellPostZeroRight(A, y)); 1031 1032 if (shell->axpy) { 1033 Mat X; 1034 PetscObjectState axpy_state; 1035 1036 PetscCall(MatShellGetContext(shell->axpy, &X)); 1037 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state)); 1038 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,...)"); 1039 PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left)); 1040 PetscCall(VecCopy(x, shell->axpy_left)); 1041 PetscCall(MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right)); 1042 PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_right)); 1043 } 1044 PetscFunctionReturn(0); 1045 } 1046 1047 static PetscErrorCode MatMultTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z) { 1048 Mat_Shell *shell = (Mat_Shell *)A->data; 1049 1050 PetscFunctionBegin; 1051 if (y == z) { 1052 if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work)); 1053 PetscCall(MatMultTranspose(A, x, shell->left_add_work)); 1054 PetscCall(VecAXPY(z, 1.0, shell->left_add_work)); 1055 } else { 1056 PetscCall(MatMultTranspose(A, x, z)); 1057 PetscCall(VecAXPY(z, 1.0, y)); 1058 } 1059 PetscFunctionReturn(0); 1060 } 1061 1062 /* 1063 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1064 */ 1065 static PetscErrorCode MatGetDiagonal_Shell(Mat A, Vec v) { 1066 Mat_Shell *shell = (Mat_Shell *)A->data; 1067 1068 PetscFunctionBegin; 1069 if (shell->ops->getdiagonal) { 1070 PetscCall((*shell->ops->getdiagonal)(A, v)); 1071 } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)"); 1072 PetscCall(VecScale(v, shell->vscale)); 1073 if (shell->dshift) PetscCall(VecAXPY(v, 1.0, shell->dshift)); 1074 PetscCall(VecShift(v, shell->vshift)); 1075 if (shell->left) PetscCall(VecPointwiseMult(v, v, shell->left)); 1076 if (shell->right) PetscCall(VecPointwiseMult(v, v, shell->right)); 1077 if (shell->zrows) { 1078 PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE)); 1079 PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE)); 1080 } 1081 if (shell->axpy) { 1082 Mat X; 1083 PetscObjectState axpy_state; 1084 1085 PetscCall(MatShellGetContext(shell->axpy, &X)); 1086 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state)); 1087 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,...)"); 1088 PetscCall(MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left)); 1089 PetscCall(MatGetDiagonal(shell->axpy, shell->axpy_left)); 1090 PetscCall(VecAXPY(v, shell->axpy_vscale, shell->axpy_left)); 1091 } 1092 PetscFunctionReturn(0); 1093 } 1094 1095 static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a) { 1096 Mat_Shell *shell = (Mat_Shell *)Y->data; 1097 PetscBool flg; 1098 1099 PetscFunctionBegin; 1100 PetscCall(MatHasCongruentLayouts(Y, &flg)); 1101 PetscCheck(flg, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot shift shell matrix if it is not congruent"); 1102 if (shell->left || shell->right) { 1103 if (!shell->dshift) { 1104 PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift)); 1105 PetscCall(VecSet(shell->dshift, a)); 1106 } else { 1107 if (shell->left) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->left)); 1108 if (shell->right) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->right)); 1109 PetscCall(VecShift(shell->dshift, a)); 1110 } 1111 if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->left)); 1112 if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->right)); 1113 } else shell->vshift += a; 1114 if (shell->zrows) PetscCall(VecShift(shell->zvals, a)); 1115 PetscFunctionReturn(0); 1116 } 1117 1118 static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s) { 1119 Mat_Shell *shell = (Mat_Shell *)A->data; 1120 1121 PetscFunctionBegin; 1122 if (!shell->dshift) PetscCall(VecDuplicate(D, &shell->dshift)); 1123 if (shell->left || shell->right) { 1124 if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work)); 1125 if (shell->left && shell->right) { 1126 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left)); 1127 PetscCall(VecPointwiseDivide(shell->right_work, shell->right_work, shell->right)); 1128 } else if (shell->left) { 1129 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left)); 1130 } else { 1131 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->right)); 1132 } 1133 PetscCall(VecAXPY(shell->dshift, s, shell->right_work)); 1134 } else { 1135 PetscCall(VecAXPY(shell->dshift, s, D)); 1136 } 1137 PetscFunctionReturn(0); 1138 } 1139 1140 PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins) { 1141 Mat_Shell *shell = (Mat_Shell *)A->data; 1142 Vec d; 1143 PetscBool flg; 1144 1145 PetscFunctionBegin; 1146 PetscCall(MatHasCongruentLayouts(A, &flg)); 1147 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot diagonal set or shift shell matrix if it is not congruent"); 1148 if (ins == INSERT_VALUES) { 1149 PetscCall(VecDuplicate(D, &d)); 1150 PetscCall(MatGetDiagonal(A, d)); 1151 PetscCall(MatDiagonalSet_Shell_Private(A, d, -1.)); 1152 PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.)); 1153 PetscCall(VecDestroy(&d)); 1154 if (shell->zrows) PetscCall(VecCopy(D, shell->zvals)); 1155 } else { 1156 PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.)); 1157 if (shell->zrows) PetscCall(VecAXPY(shell->zvals, 1.0, D)); 1158 } 1159 PetscFunctionReturn(0); 1160 } 1161 1162 static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a) { 1163 Mat_Shell *shell = (Mat_Shell *)Y->data; 1164 1165 PetscFunctionBegin; 1166 shell->vscale *= a; 1167 shell->vshift *= a; 1168 if (shell->dshift) PetscCall(VecScale(shell->dshift, a)); 1169 shell->axpy_vscale *= a; 1170 if (shell->zrows) PetscCall(VecScale(shell->zvals, a)); 1171 PetscFunctionReturn(0); 1172 } 1173 1174 static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right) { 1175 Mat_Shell *shell = (Mat_Shell *)Y->data; 1176 1177 PetscFunctionBegin; 1178 if (left) { 1179 if (!shell->left) { 1180 PetscCall(VecDuplicate(left, &shell->left)); 1181 PetscCall(VecCopy(left, shell->left)); 1182 } else { 1183 PetscCall(VecPointwiseMult(shell->left, shell->left, left)); 1184 } 1185 if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, left)); 1186 } 1187 if (right) { 1188 if (!shell->right) { 1189 PetscCall(VecDuplicate(right, &shell->right)); 1190 PetscCall(VecCopy(right, shell->right)); 1191 } else { 1192 PetscCall(VecPointwiseMult(shell->right, shell->right, right)); 1193 } 1194 if (shell->zrows) { 1195 if (!shell->left_work) PetscCall(MatCreateVecs(Y, NULL, &shell->left_work)); 1196 PetscCall(VecSet(shell->zvals_w, 1.0)); 1197 PetscCall(VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD)); 1198 PetscCall(VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD)); 1199 PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w)); 1200 } 1201 } 1202 if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy, left, right)); 1203 PetscFunctionReturn(0); 1204 } 1205 1206 PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t) { 1207 Mat_Shell *shell = (Mat_Shell *)Y->data; 1208 1209 PetscFunctionBegin; 1210 if (t == MAT_FINAL_ASSEMBLY) { 1211 shell->vshift = 0.0; 1212 shell->vscale = 1.0; 1213 shell->axpy_vscale = 0.0; 1214 shell->axpy_state = 0; 1215 PetscCall(VecDestroy(&shell->dshift)); 1216 PetscCall(VecDestroy(&shell->left)); 1217 PetscCall(VecDestroy(&shell->right)); 1218 PetscCall(MatDestroy(&shell->axpy)); 1219 PetscCall(VecDestroy(&shell->axpy_left)); 1220 PetscCall(VecDestroy(&shell->axpy_right)); 1221 PetscCall(VecScatterDestroy(&shell->zvals_sct_c)); 1222 PetscCall(VecScatterDestroy(&shell->zvals_sct_r)); 1223 PetscCall(ISDestroy(&shell->zrows)); 1224 PetscCall(ISDestroy(&shell->zcols)); 1225 } 1226 PetscFunctionReturn(0); 1227 } 1228 1229 static PetscErrorCode MatMissingDiagonal_Shell(Mat A, PetscBool *missing, PetscInt *d) { 1230 PetscFunctionBegin; 1231 *missing = PETSC_FALSE; 1232 PetscFunctionReturn(0); 1233 } 1234 1235 PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str) { 1236 Mat_Shell *shell = (Mat_Shell *)Y->data; 1237 1238 PetscFunctionBegin; 1239 if (X == Y) { 1240 PetscCall(MatScale(Y, 1.0 + a)); 1241 PetscFunctionReturn(0); 1242 } 1243 if (!shell->axpy) { 1244 PetscCall(MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy)); 1245 shell->axpy_vscale = a; 1246 PetscCall(PetscObjectStateGet((PetscObject)X, &shell->axpy_state)); 1247 } else { 1248 PetscCall(MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str)); 1249 } 1250 PetscFunctionReturn(0); 1251 } 1252 1253 static struct _MatOps MatOps_Values = {NULL, 1254 NULL, 1255 NULL, 1256 NULL, 1257 /* 4*/ MatMultAdd_Shell, 1258 NULL, 1259 MatMultTransposeAdd_Shell, 1260 NULL, 1261 NULL, 1262 NULL, 1263 /*10*/ NULL, 1264 NULL, 1265 NULL, 1266 NULL, 1267 NULL, 1268 /*15*/ NULL, 1269 NULL, 1270 NULL, 1271 MatDiagonalScale_Shell, 1272 NULL, 1273 /*20*/ NULL, 1274 MatAssemblyEnd_Shell, 1275 NULL, 1276 NULL, 1277 /*24*/ MatZeroRows_Shell, 1278 NULL, 1279 NULL, 1280 NULL, 1281 NULL, 1282 /*29*/ NULL, 1283 NULL, 1284 NULL, 1285 NULL, 1286 NULL, 1287 /*34*/ MatDuplicate_Shell, 1288 NULL, 1289 NULL, 1290 NULL, 1291 NULL, 1292 /*39*/ MatAXPY_Shell, 1293 NULL, 1294 NULL, 1295 NULL, 1296 MatCopy_Shell, 1297 /*44*/ NULL, 1298 MatScale_Shell, 1299 MatShift_Shell, 1300 MatDiagonalSet_Shell, 1301 MatZeroRowsColumns_Shell, 1302 /*49*/ NULL, 1303 NULL, 1304 NULL, 1305 NULL, 1306 NULL, 1307 /*54*/ NULL, 1308 NULL, 1309 NULL, 1310 NULL, 1311 NULL, 1312 /*59*/ NULL, 1313 MatDestroy_Shell, 1314 NULL, 1315 MatConvertFrom_Shell, 1316 NULL, 1317 /*64*/ NULL, 1318 NULL, 1319 NULL, 1320 NULL, 1321 NULL, 1322 /*69*/ NULL, 1323 NULL, 1324 MatConvert_Shell, 1325 NULL, 1326 NULL, 1327 /*74*/ NULL, 1328 NULL, 1329 NULL, 1330 NULL, 1331 NULL, 1332 /*79*/ NULL, 1333 NULL, 1334 NULL, 1335 NULL, 1336 NULL, 1337 /*84*/ NULL, 1338 NULL, 1339 NULL, 1340 NULL, 1341 NULL, 1342 /*89*/ NULL, 1343 NULL, 1344 NULL, 1345 NULL, 1346 NULL, 1347 /*94*/ NULL, 1348 NULL, 1349 NULL, 1350 NULL, 1351 NULL, 1352 /*99*/ NULL, 1353 NULL, 1354 NULL, 1355 NULL, 1356 NULL, 1357 /*104*/ NULL, 1358 NULL, 1359 NULL, 1360 NULL, 1361 NULL, 1362 /*109*/ NULL, 1363 NULL, 1364 NULL, 1365 NULL, 1366 MatMissingDiagonal_Shell, 1367 /*114*/ NULL, 1368 NULL, 1369 NULL, 1370 NULL, 1371 NULL, 1372 /*119*/ NULL, 1373 NULL, 1374 NULL, 1375 NULL, 1376 NULL, 1377 /*124*/ NULL, 1378 NULL, 1379 NULL, 1380 NULL, 1381 NULL, 1382 /*129*/ NULL, 1383 NULL, 1384 NULL, 1385 NULL, 1386 NULL, 1387 /*134*/ NULL, 1388 NULL, 1389 NULL, 1390 NULL, 1391 NULL, 1392 /*139*/ NULL, 1393 NULL, 1394 NULL, 1395 NULL, 1396 NULL, 1397 /*144*/ NULL, 1398 NULL, 1399 NULL, 1400 NULL, 1401 NULL, 1402 NULL, 1403 /*150*/ NULL}; 1404 1405 static PetscErrorCode MatShellSetContext_Shell(Mat mat, void *ctx) { 1406 Mat_Shell *shell = (Mat_Shell *)mat->data; 1407 1408 PetscFunctionBegin; 1409 if (ctx) { 1410 PetscContainer ctxcontainer; 1411 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer)); 1412 PetscCall(PetscContainerSetPointer(ctxcontainer, ctx)); 1413 PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer)); 1414 shell->ctxcontainer = ctxcontainer; 1415 PetscCall(PetscContainerDestroy(&ctxcontainer)); 1416 } else { 1417 PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL)); 1418 shell->ctxcontainer = NULL; 1419 } 1420 1421 PetscFunctionReturn(0); 1422 } 1423 1424 PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscErrorCode (*f)(void *)) { 1425 Mat_Shell *shell = (Mat_Shell *)mat->data; 1426 1427 PetscFunctionBegin; 1428 if (shell->ctxcontainer) PetscCall(PetscContainerSetUserDestroy(shell->ctxcontainer, f)); 1429 PetscFunctionReturn(0); 1430 } 1431 1432 static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype) { 1433 PetscFunctionBegin; 1434 PetscCall(PetscFree(mat->defaultvectype)); 1435 PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype)); 1436 PetscFunctionReturn(0); 1437 } 1438 1439 PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A) { 1440 Mat_Shell *shell = (Mat_Shell *)A->data; 1441 1442 PetscFunctionBegin; 1443 shell->managescalingshifts = PETSC_FALSE; 1444 A->ops->diagonalset = NULL; 1445 A->ops->diagonalscale = NULL; 1446 A->ops->scale = NULL; 1447 A->ops->shift = NULL; 1448 A->ops->axpy = NULL; 1449 PetscFunctionReturn(0); 1450 } 1451 1452 static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, void (*f)(void)) { 1453 Mat_Shell *shell = (Mat_Shell *)mat->data; 1454 1455 PetscFunctionBegin; 1456 switch (op) { 1457 case MATOP_DESTROY: shell->ops->destroy = (PetscErrorCode(*)(Mat))f; break; 1458 case MATOP_VIEW: 1459 if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view; 1460 mat->ops->view = (PetscErrorCode(*)(Mat, PetscViewer))f; 1461 break; 1462 case MATOP_COPY: shell->ops->copy = (PetscErrorCode(*)(Mat, Mat, MatStructure))f; break; 1463 case MATOP_DIAGONAL_SET: 1464 case MATOP_DIAGONAL_SCALE: 1465 case MATOP_SHIFT: 1466 case MATOP_SCALE: 1467 case MATOP_AXPY: 1468 case MATOP_ZERO_ROWS: 1469 case MATOP_ZERO_ROWS_COLUMNS: 1470 PetscCheck(!shell->managescalingshifts, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1471 (((void (**)(void))mat->ops)[op]) = f; 1472 break; 1473 case MATOP_GET_DIAGONAL: 1474 if (shell->managescalingshifts) { 1475 shell->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f; 1476 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1477 } else { 1478 shell->ops->getdiagonal = NULL; 1479 mat->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f; 1480 } 1481 break; 1482 case MATOP_MULT: 1483 if (shell->managescalingshifts) { 1484 shell->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f; 1485 mat->ops->mult = MatMult_Shell; 1486 } else { 1487 shell->ops->mult = NULL; 1488 mat->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f; 1489 } 1490 break; 1491 case MATOP_MULT_TRANSPOSE: 1492 if (shell->managescalingshifts) { 1493 shell->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f; 1494 mat->ops->multtranspose = MatMultTranspose_Shell; 1495 } else { 1496 shell->ops->multtranspose = NULL; 1497 mat->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f; 1498 } 1499 break; 1500 default: (((void (**)(void))mat->ops)[op]) = f; break; 1501 } 1502 PetscFunctionReturn(0); 1503 } 1504 1505 static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, void (**f)(void)) { 1506 Mat_Shell *shell = (Mat_Shell *)mat->data; 1507 1508 PetscFunctionBegin; 1509 switch (op) { 1510 case MATOP_DESTROY: *f = (void (*)(void))shell->ops->destroy; break; 1511 case MATOP_VIEW: *f = (void (*)(void))mat->ops->view; break; 1512 case MATOP_COPY: *f = (void (*)(void))shell->ops->copy; break; 1513 case MATOP_DIAGONAL_SET: 1514 case MATOP_DIAGONAL_SCALE: 1515 case MATOP_SHIFT: 1516 case MATOP_SCALE: 1517 case MATOP_AXPY: 1518 case MATOP_ZERO_ROWS: 1519 case MATOP_ZERO_ROWS_COLUMNS: *f = (((void (**)(void))mat->ops)[op]); break; 1520 case MATOP_GET_DIAGONAL: 1521 if (shell->ops->getdiagonal) *f = (void (*)(void))shell->ops->getdiagonal; 1522 else *f = (((void (**)(void))mat->ops)[op]); 1523 break; 1524 case MATOP_MULT: 1525 if (shell->ops->mult) *f = (void (*)(void))shell->ops->mult; 1526 else *f = (((void (**)(void))mat->ops)[op]); 1527 break; 1528 case MATOP_MULT_TRANSPOSE: 1529 if (shell->ops->multtranspose) *f = (void (*)(void))shell->ops->multtranspose; 1530 else *f = (((void (**)(void))mat->ops)[op]); 1531 break; 1532 default: *f = (((void (**)(void))mat->ops)[op]); 1533 } 1534 PetscFunctionReturn(0); 1535 } 1536 1537 /*MC 1538 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 1539 1540 Level: advanced 1541 1542 .seealso: `MatCreateShell()` 1543 M*/ 1544 1545 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) { 1546 Mat_Shell *b; 1547 1548 PetscFunctionBegin; 1549 PetscCall(PetscMemcpy(A->ops, &MatOps_Values, sizeof(struct _MatOps))); 1550 1551 PetscCall(PetscNewLog(A, &b)); 1552 A->data = (void *)b; 1553 1554 b->ctxcontainer = NULL; 1555 b->vshift = 0.0; 1556 b->vscale = 1.0; 1557 b->managescalingshifts = PETSC_TRUE; 1558 A->assembled = PETSC_TRUE; 1559 A->preallocated = PETSC_FALSE; 1560 1561 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell)); 1562 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell)); 1563 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell)); 1564 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell)); 1565 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell)); 1566 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell)); 1567 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell)); 1568 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell)); 1569 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL)); 1570 PetscFunctionReturn(0); 1571 } 1572 1573 /*@C 1574 MatCreateShell - Creates a new matrix class for use with a user-defined 1575 private data storage format. 1576 1577 Collective 1578 1579 Input Parameters: 1580 + comm - MPI communicator 1581 . m - number of local rows (must be given) 1582 . n - number of local columns (must be given) 1583 . M - number of global rows (may be PETSC_DETERMINE) 1584 . N - number of global columns (may be PETSC_DETERMINE) 1585 - ctx - pointer to data needed by the shell matrix routines 1586 1587 Output Parameter: 1588 . A - the matrix 1589 1590 Level: advanced 1591 1592 Usage: 1593 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1594 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 1595 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1596 $ [ Use matrix for operations that have been set ] 1597 $ MatDestroy(mat); 1598 1599 Notes: 1600 The shell matrix type is intended to provide a simple class to use 1601 with KSP (such as, for use with matrix-free methods). You should not 1602 use the shell type if you plan to define a complete matrix class. 1603 1604 Fortran Notes: 1605 To use this from Fortran with a ctx you must write an interface definition for this 1606 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 1607 in as the ctx argument. 1608 1609 PETSc requires that matrices and vectors being used for certain 1610 operations are partitioned accordingly. For example, when 1611 creating a shell matrix, A, that supports parallel matrix-vector 1612 products using MatMult(A,x,y) the user should set the number 1613 of local matrix rows to be the number of local elements of the 1614 corresponding result vector, y. Note that this is information is 1615 required for use of the matrix interface routines, even though 1616 the shell matrix may not actually be physically partitioned. 1617 For example, 1618 1619 $ 1620 $ Vec x, y 1621 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1622 $ Mat A 1623 $ 1624 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 1625 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 1626 $ VecGetLocalSize(y,&m); 1627 $ VecGetLocalSize(x,&n); 1628 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 1629 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1630 $ MatMult(A,x,y); 1631 $ MatDestroy(&A); 1632 $ VecDestroy(&y); 1633 $ VecDestroy(&x); 1634 $ 1635 1636 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these 1637 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 1638 1639 For rectangular matrices do all the scalings and shifts make sense? 1640 1641 Developers Notes: 1642 Regarding shifting and scaling. The general form is 1643 1644 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1645 1646 The order you apply the operations is important. For example if you have a dshift then 1647 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 1648 you get s*vscale*A + diag(shift) 1649 1650 A is the user provided function. 1651 1652 KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for 1653 for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger 1654 an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat); 1655 each time the MATSHELL matrix has changed. 1656 1657 Matrix product operations (i.e. MatMat, MatTransposeMat etc) can be specified using MatShellSetMatProductOperation() 1658 1659 Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided 1660 with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale(). 1661 1662 .seealso: `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MATSHELL`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()` 1663 @*/ 1664 PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, void *ctx, Mat *A) { 1665 PetscFunctionBegin; 1666 PetscCall(MatCreate(comm, A)); 1667 PetscCall(MatSetSizes(*A, m, n, M, N)); 1668 PetscCall(MatSetType(*A, MATSHELL)); 1669 PetscCall(MatShellSetContext(*A, ctx)); 1670 PetscCall(MatSetUp(*A)); 1671 PetscFunctionReturn(0); 1672 } 1673 1674 /*@ 1675 MatShellSetContext - sets the context for a shell matrix 1676 1677 Logically Collective on Mat 1678 1679 Input Parameters: 1680 + mat - the shell matrix 1681 - ctx - the context 1682 1683 Level: advanced 1684 1685 Fortran Notes: 1686 To use this from Fortran you must write a Fortran interface definition for this 1687 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 1688 1689 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()` 1690 @*/ 1691 PetscErrorCode MatShellSetContext(Mat mat, void *ctx) { 1692 PetscFunctionBegin; 1693 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 1694 PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx)); 1695 PetscFunctionReturn(0); 1696 } 1697 1698 /*@ 1699 MatShellSetContextDestroy - sets the destroy function for a shell matrix context 1700 1701 Logically Collective on Mat 1702 1703 Input Parameters: 1704 + mat - the shell matrix 1705 - f - the context destroy function 1706 1707 Note: 1708 If the `MatShell` is never duplicated, the behavior of this function is equivalent 1709 to `MatShellSetOperation(Mat,MATOP_DESTROY,f)`. However, `MatShellSetContextDestroy()` 1710 ensures proper reference counting for the user provided context data in the case that 1711 the `MatShell` is duplicated. 1712 1713 Level: advanced 1714 1715 .seealso: `MatCreateShell()`, `MatShellSetContext()` 1716 @*/ 1717 PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscErrorCode (*f)(void *)) { 1718 PetscFunctionBegin; 1719 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 1720 PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscErrorCode(*)(void *)), (mat, f)); 1721 PetscFunctionReturn(0); 1722 } 1723 1724 /*@C 1725 MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs() 1726 1727 Logically collective 1728 1729 Input Parameters: 1730 + mat - the shell matrix 1731 - vtype - type to use for creating vectors 1732 1733 Notes: 1734 1735 Level: advanced 1736 1737 .seealso: `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 Mat 1750 1751 Input Parameter: 1752 . mat - the shell matrix 1753 1754 Level: advanced 1755 1756 .seealso: `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 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 Notes: 1785 Not supported from Fortran 1786 1787 .seealso: `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 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 Notes: 1849 Not supported from Fortran 1850 1851 .seealso: `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 shell matrix. 1902 1903 Logically Collective on Mat 1904 1905 Input Parameters: 1906 + mat - the 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. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation() 1937 1938 Fortran Notes: 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: `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 shell matrix. 1953 1954 Not Collective 1955 1956 Input Parameters: 1957 + mat - the 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 Notes: 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: `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 Notes: in the future, we should allow the object type name to be changed still using the MatShell data structure for other matrices (i.e. MATTRANSPOSEMAT, MATSCHURCOMPLEMENT etc) 2002 2003 .seealso: `MatCreateShell()` 2004 @*/ 2005 PetscErrorCode MatIsShell(Mat mat, PetscBool *flg) { 2006 PetscFunctionBegin; 2007 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2008 PetscValidBoolPointer(flg, 2); 2009 *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell); 2010 PetscFunctionReturn(0); 2011 } 2012