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