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