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 PetscFunctionReturn(PETSC_SUCCESS); 1004 } 1005 1006 static PetscErrorCode MatMult_Shell(Mat A, Vec x, Vec y) 1007 { 1008 Mat_Shell *shell = (Mat_Shell *)A->data; 1009 Vec xx; 1010 PetscObjectState instate, outstate; 1011 1012 PetscFunctionBegin; 1013 PetscCall(MatShellPreZeroRight(A, x, &xx)); 1014 PetscCall(MatShellPreScaleRight(A, xx, &xx)); 1015 PetscCall(PetscObjectStateGet((PetscObject)y, &instate)); 1016 PetscCall((*shell->ops->mult)(A, xx, y)); 1017 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate)); 1018 if (instate == outstate) { 1019 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1020 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 1021 } 1022 PetscCall(MatShellShiftAndScale(A, xx, y)); 1023 PetscCall(MatShellPostScaleLeft(A, y)); 1024 PetscCall(MatShellPostZeroLeft(A, y)); 1025 1026 if (shell->axpy) { 1027 Mat X; 1028 PetscObjectState axpy_state; 1029 1030 PetscCall(MatShellGetContext(shell->axpy, &X)); 1031 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state)); 1032 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,...)"); 1033 1034 PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left)); 1035 PetscCall(VecCopy(x, shell->axpy_right)); 1036 PetscCall(MatMult(shell->axpy, shell->axpy_right, shell->axpy_left)); 1037 PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_left)); 1038 } 1039 PetscFunctionReturn(PETSC_SUCCESS); 1040 } 1041 1042 static PetscErrorCode MatMultAdd_Shell(Mat A, Vec x, Vec y, Vec z) 1043 { 1044 Mat_Shell *shell = (Mat_Shell *)A->data; 1045 1046 PetscFunctionBegin; 1047 if (y == z) { 1048 if (!shell->right_add_work) PetscCall(VecDuplicate(z, &shell->right_add_work)); 1049 PetscCall(MatMult(A, x, shell->right_add_work)); 1050 PetscCall(VecAXPY(z, 1.0, shell->right_add_work)); 1051 } else { 1052 PetscCall(MatMult(A, x, z)); 1053 PetscCall(VecAXPY(z, 1.0, y)); 1054 } 1055 PetscFunctionReturn(PETSC_SUCCESS); 1056 } 1057 1058 static PetscErrorCode MatMultTranspose_Shell(Mat A, Vec x, Vec y) 1059 { 1060 Mat_Shell *shell = (Mat_Shell *)A->data; 1061 Vec xx; 1062 PetscObjectState instate, outstate; 1063 1064 PetscFunctionBegin; 1065 PetscCall(MatShellPreZeroLeft(A, x, &xx)); 1066 PetscCall(MatShellPreScaleLeft(A, xx, &xx)); 1067 PetscCall(PetscObjectStateGet((PetscObject)y, &instate)); 1068 PetscCall((*shell->ops->multtranspose)(A, xx, y)); 1069 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate)); 1070 if (instate == outstate) { 1071 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1072 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 1073 } 1074 PetscCall(MatShellShiftAndScale(A, xx, y)); 1075 PetscCall(MatShellPostScaleRight(A, y)); 1076 PetscCall(MatShellPostZeroRight(A, y)); 1077 1078 if (shell->axpy) { 1079 Mat X; 1080 PetscObjectState axpy_state; 1081 1082 PetscCall(MatShellGetContext(shell->axpy, &X)); 1083 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state)); 1084 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,...)"); 1085 PetscCall(MatCreateVecs(shell->axpy, shell->axpy_right ? NULL : &shell->axpy_right, shell->axpy_left ? NULL : &shell->axpy_left)); 1086 PetscCall(VecCopy(x, shell->axpy_left)); 1087 PetscCall(MatMultTranspose(shell->axpy, shell->axpy_left, shell->axpy_right)); 1088 PetscCall(VecAXPY(y, shell->axpy_vscale, shell->axpy_right)); 1089 } 1090 PetscFunctionReturn(PETSC_SUCCESS); 1091 } 1092 1093 static PetscErrorCode MatMultTransposeAdd_Shell(Mat A, Vec x, Vec y, Vec z) 1094 { 1095 Mat_Shell *shell = (Mat_Shell *)A->data; 1096 1097 PetscFunctionBegin; 1098 if (y == z) { 1099 if (!shell->left_add_work) PetscCall(VecDuplicate(z, &shell->left_add_work)); 1100 PetscCall(MatMultTranspose(A, x, shell->left_add_work)); 1101 PetscCall(VecAXPY(z, 1.0, shell->left_add_work)); 1102 } else { 1103 PetscCall(MatMultTranspose(A, x, z)); 1104 PetscCall(VecAXPY(z, 1.0, y)); 1105 } 1106 PetscFunctionReturn(PETSC_SUCCESS); 1107 } 1108 1109 /* 1110 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1111 */ 1112 static PetscErrorCode MatGetDiagonal_Shell(Mat A, Vec v) 1113 { 1114 Mat_Shell *shell = (Mat_Shell *)A->data; 1115 1116 PetscFunctionBegin; 1117 if (shell->ops->getdiagonal) { 1118 PetscCall((*shell->ops->getdiagonal)(A, v)); 1119 } else SETERRQ(PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONGSTATE, "Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)"); 1120 PetscCall(VecScale(v, shell->vscale)); 1121 if (shell->dshift) PetscCall(VecAXPY(v, 1.0, shell->dshift)); 1122 PetscCall(VecShift(v, shell->vshift)); 1123 if (shell->left) PetscCall(VecPointwiseMult(v, v, shell->left)); 1124 if (shell->right) PetscCall(VecPointwiseMult(v, v, shell->right)); 1125 if (shell->zrows) { 1126 PetscCall(VecScatterBegin(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE)); 1127 PetscCall(VecScatterEnd(shell->zvals_sct_r, shell->zvals, v, INSERT_VALUES, SCATTER_REVERSE)); 1128 } 1129 if (shell->axpy) { 1130 Mat X; 1131 PetscObjectState axpy_state; 1132 1133 PetscCall(MatShellGetContext(shell->axpy, &X)); 1134 PetscCall(PetscObjectStateGet((PetscObject)X, &axpy_state)); 1135 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,...)"); 1136 PetscCall(MatCreateVecs(shell->axpy, NULL, shell->axpy_left ? NULL : &shell->axpy_left)); 1137 PetscCall(MatGetDiagonal(shell->axpy, shell->axpy_left)); 1138 PetscCall(VecAXPY(v, shell->axpy_vscale, shell->axpy_left)); 1139 } 1140 PetscFunctionReturn(PETSC_SUCCESS); 1141 } 1142 1143 static PetscErrorCode MatShift_Shell(Mat Y, PetscScalar a) 1144 { 1145 Mat_Shell *shell = (Mat_Shell *)Y->data; 1146 PetscBool flg; 1147 1148 PetscFunctionBegin; 1149 PetscCall(MatHasCongruentLayouts(Y, &flg)); 1150 PetscCheck(flg, PetscObjectComm((PetscObject)Y), PETSC_ERR_SUP, "Cannot shift shell matrix if it is not congruent"); 1151 if (shell->left || shell->right) { 1152 if (!shell->dshift) { 1153 PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift)); 1154 PetscCall(VecSet(shell->dshift, a)); 1155 } else { 1156 if (shell->left) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->left)); 1157 if (shell->right) PetscCall(VecPointwiseMult(shell->dshift, shell->dshift, shell->right)); 1158 PetscCall(VecShift(shell->dshift, a)); 1159 } 1160 if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->left)); 1161 if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift, shell->dshift, shell->right)); 1162 } else shell->vshift += a; 1163 if (shell->zrows) PetscCall(VecShift(shell->zvals, a)); 1164 PetscFunctionReturn(PETSC_SUCCESS); 1165 } 1166 1167 static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A, Vec D, PetscScalar s) 1168 { 1169 Mat_Shell *shell = (Mat_Shell *)A->data; 1170 1171 PetscFunctionBegin; 1172 if (!shell->dshift) PetscCall(VecDuplicate(D, &shell->dshift)); 1173 if (shell->left || shell->right) { 1174 if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work)); 1175 if (shell->left && shell->right) { 1176 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left)); 1177 PetscCall(VecPointwiseDivide(shell->right_work, shell->right_work, shell->right)); 1178 } else if (shell->left) { 1179 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->left)); 1180 } else { 1181 PetscCall(VecPointwiseDivide(shell->right_work, D, shell->right)); 1182 } 1183 PetscCall(VecAXPY(shell->dshift, s, shell->right_work)); 1184 } else { 1185 PetscCall(VecAXPY(shell->dshift, s, D)); 1186 } 1187 PetscFunctionReturn(PETSC_SUCCESS); 1188 } 1189 1190 static PetscErrorCode MatDiagonalSet_Shell(Mat A, Vec D, InsertMode ins) 1191 { 1192 Mat_Shell *shell = (Mat_Shell *)A->data; 1193 Vec d; 1194 PetscBool flg; 1195 1196 PetscFunctionBegin; 1197 PetscCall(MatHasCongruentLayouts(A, &flg)); 1198 PetscCheck(flg, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot diagonal set or shift shell matrix if it is not congruent"); 1199 if (ins == INSERT_VALUES) { 1200 PetscCall(VecDuplicate(D, &d)); 1201 PetscCall(MatGetDiagonal(A, d)); 1202 PetscCall(MatDiagonalSet_Shell_Private(A, d, -1.)); 1203 PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.)); 1204 PetscCall(VecDestroy(&d)); 1205 if (shell->zrows) PetscCall(VecCopy(D, shell->zvals)); 1206 } else { 1207 PetscCall(MatDiagonalSet_Shell_Private(A, D, 1.)); 1208 if (shell->zrows) PetscCall(VecAXPY(shell->zvals, 1.0, D)); 1209 } 1210 PetscFunctionReturn(PETSC_SUCCESS); 1211 } 1212 1213 static PetscErrorCode MatScale_Shell(Mat Y, PetscScalar a) 1214 { 1215 Mat_Shell *shell = (Mat_Shell *)Y->data; 1216 1217 PetscFunctionBegin; 1218 shell->vscale *= a; 1219 shell->vshift *= a; 1220 if (shell->dshift) PetscCall(VecScale(shell->dshift, a)); 1221 shell->axpy_vscale *= a; 1222 if (shell->zrows) PetscCall(VecScale(shell->zvals, a)); 1223 PetscFunctionReturn(PETSC_SUCCESS); 1224 } 1225 1226 static PetscErrorCode MatDiagonalScale_Shell(Mat Y, Vec left, Vec right) 1227 { 1228 Mat_Shell *shell = (Mat_Shell *)Y->data; 1229 1230 PetscFunctionBegin; 1231 if (left) { 1232 if (!shell->left) { 1233 PetscCall(VecDuplicate(left, &shell->left)); 1234 PetscCall(VecCopy(left, shell->left)); 1235 } else { 1236 PetscCall(VecPointwiseMult(shell->left, shell->left, left)); 1237 } 1238 if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, left)); 1239 } 1240 if (right) { 1241 if (!shell->right) { 1242 PetscCall(VecDuplicate(right, &shell->right)); 1243 PetscCall(VecCopy(right, shell->right)); 1244 } else { 1245 PetscCall(VecPointwiseMult(shell->right, shell->right, right)); 1246 } 1247 if (shell->zrows) { 1248 if (!shell->left_work) PetscCall(MatCreateVecs(Y, NULL, &shell->left_work)); 1249 PetscCall(VecSet(shell->zvals_w, 1.0)); 1250 PetscCall(VecScatterBegin(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD)); 1251 PetscCall(VecScatterEnd(shell->zvals_sct_c, right, shell->zvals_w, INSERT_VALUES, SCATTER_FORWARD)); 1252 PetscCall(VecPointwiseMult(shell->zvals, shell->zvals, shell->zvals_w)); 1253 } 1254 } 1255 if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy, left, right)); 1256 PetscFunctionReturn(PETSC_SUCCESS); 1257 } 1258 1259 static PetscErrorCode MatAssemblyEnd_Shell(Mat Y, MatAssemblyType t) 1260 { 1261 Mat_Shell *shell = (Mat_Shell *)Y->data; 1262 1263 PetscFunctionBegin; 1264 if (t == MAT_FINAL_ASSEMBLY) { 1265 shell->vshift = 0.0; 1266 shell->vscale = 1.0; 1267 shell->axpy_vscale = 0.0; 1268 shell->axpy_state = 0; 1269 PetscCall(VecDestroy(&shell->dshift)); 1270 PetscCall(VecDestroy(&shell->left)); 1271 PetscCall(VecDestroy(&shell->right)); 1272 PetscCall(MatDestroy(&shell->axpy)); 1273 PetscCall(VecDestroy(&shell->axpy_left)); 1274 PetscCall(VecDestroy(&shell->axpy_right)); 1275 PetscCall(VecScatterDestroy(&shell->zvals_sct_c)); 1276 PetscCall(VecScatterDestroy(&shell->zvals_sct_r)); 1277 PetscCall(ISDestroy(&shell->zrows)); 1278 PetscCall(ISDestroy(&shell->zcols)); 1279 } 1280 PetscFunctionReturn(PETSC_SUCCESS); 1281 } 1282 1283 static PetscErrorCode MatMissingDiagonal_Shell(Mat A, PetscBool *missing, PetscInt *d) 1284 { 1285 PetscFunctionBegin; 1286 *missing = PETSC_FALSE; 1287 PetscFunctionReturn(PETSC_SUCCESS); 1288 } 1289 1290 static PetscErrorCode MatAXPY_Shell(Mat Y, PetscScalar a, Mat X, MatStructure str) 1291 { 1292 Mat_Shell *shell = (Mat_Shell *)Y->data; 1293 1294 PetscFunctionBegin; 1295 if (X == Y) { 1296 PetscCall(MatScale(Y, 1.0 + a)); 1297 PetscFunctionReturn(PETSC_SUCCESS); 1298 } 1299 if (!shell->axpy) { 1300 PetscCall(MatConvertFrom_Shell(X, MATSHELL, MAT_INITIAL_MATRIX, &shell->axpy)); 1301 shell->axpy_vscale = a; 1302 PetscCall(PetscObjectStateGet((PetscObject)X, &shell->axpy_state)); 1303 } else { 1304 PetscCall(MatAXPY(shell->axpy, a / shell->axpy_vscale, X, str)); 1305 } 1306 PetscFunctionReturn(PETSC_SUCCESS); 1307 } 1308 1309 static struct _MatOps MatOps_Values = {NULL, 1310 NULL, 1311 NULL, 1312 NULL, 1313 /* 4*/ MatMultAdd_Shell, 1314 NULL, 1315 MatMultTransposeAdd_Shell, 1316 NULL, 1317 NULL, 1318 NULL, 1319 /*10*/ NULL, 1320 NULL, 1321 NULL, 1322 NULL, 1323 NULL, 1324 /*15*/ NULL, 1325 NULL, 1326 NULL, 1327 MatDiagonalScale_Shell, 1328 NULL, 1329 /*20*/ NULL, 1330 MatAssemblyEnd_Shell, 1331 NULL, 1332 NULL, 1333 /*24*/ MatZeroRows_Shell, 1334 NULL, 1335 NULL, 1336 NULL, 1337 NULL, 1338 /*29*/ NULL, 1339 NULL, 1340 NULL, 1341 NULL, 1342 NULL, 1343 /*34*/ MatDuplicate_Shell, 1344 NULL, 1345 NULL, 1346 NULL, 1347 NULL, 1348 /*39*/ MatAXPY_Shell, 1349 NULL, 1350 NULL, 1351 NULL, 1352 MatCopy_Shell, 1353 /*44*/ NULL, 1354 MatScale_Shell, 1355 MatShift_Shell, 1356 MatDiagonalSet_Shell, 1357 MatZeroRowsColumns_Shell, 1358 /*49*/ NULL, 1359 NULL, 1360 NULL, 1361 NULL, 1362 NULL, 1363 /*54*/ NULL, 1364 NULL, 1365 NULL, 1366 NULL, 1367 NULL, 1368 /*59*/ NULL, 1369 MatDestroy_Shell, 1370 NULL, 1371 MatConvertFrom_Shell, 1372 NULL, 1373 /*64*/ NULL, 1374 NULL, 1375 NULL, 1376 NULL, 1377 NULL, 1378 /*69*/ NULL, 1379 NULL, 1380 MatConvert_Shell, 1381 NULL, 1382 NULL, 1383 /*74*/ NULL, 1384 NULL, 1385 NULL, 1386 NULL, 1387 NULL, 1388 /*79*/ NULL, 1389 NULL, 1390 NULL, 1391 NULL, 1392 NULL, 1393 /*84*/ NULL, 1394 NULL, 1395 NULL, 1396 NULL, 1397 NULL, 1398 /*89*/ NULL, 1399 NULL, 1400 NULL, 1401 NULL, 1402 NULL, 1403 /*94*/ NULL, 1404 NULL, 1405 NULL, 1406 NULL, 1407 NULL, 1408 /*99*/ NULL, 1409 NULL, 1410 NULL, 1411 NULL, 1412 NULL, 1413 /*104*/ NULL, 1414 NULL, 1415 NULL, 1416 NULL, 1417 NULL, 1418 /*109*/ NULL, 1419 NULL, 1420 NULL, 1421 NULL, 1422 MatMissingDiagonal_Shell, 1423 /*114*/ NULL, 1424 NULL, 1425 NULL, 1426 NULL, 1427 NULL, 1428 /*119*/ NULL, 1429 NULL, 1430 NULL, 1431 NULL, 1432 NULL, 1433 /*124*/ NULL, 1434 NULL, 1435 NULL, 1436 NULL, 1437 NULL, 1438 /*129*/ NULL, 1439 NULL, 1440 NULL, 1441 NULL, 1442 NULL, 1443 /*134*/ NULL, 1444 NULL, 1445 NULL, 1446 NULL, 1447 NULL, 1448 /*139*/ NULL, 1449 NULL, 1450 NULL, 1451 NULL, 1452 NULL, 1453 /*144*/ NULL, 1454 NULL, 1455 NULL, 1456 NULL, 1457 NULL, 1458 NULL, 1459 /*150*/ NULL, 1460 NULL}; 1461 1462 static PetscErrorCode MatShellSetContext_Shell(Mat mat, void *ctx) 1463 { 1464 Mat_Shell *shell = (Mat_Shell *)mat->data; 1465 1466 PetscFunctionBegin; 1467 if (ctx) { 1468 PetscContainer ctxcontainer; 1469 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat), &ctxcontainer)); 1470 PetscCall(PetscContainerSetPointer(ctxcontainer, ctx)); 1471 PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", (PetscObject)ctxcontainer)); 1472 shell->ctxcontainer = ctxcontainer; 1473 PetscCall(PetscContainerDestroy(&ctxcontainer)); 1474 } else { 1475 PetscCall(PetscObjectCompose((PetscObject)mat, "MatShell ctx", NULL)); 1476 shell->ctxcontainer = NULL; 1477 } 1478 1479 PetscFunctionReturn(PETSC_SUCCESS); 1480 } 1481 1482 static PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat, PetscErrorCode (*f)(void *)) 1483 { 1484 Mat_Shell *shell = (Mat_Shell *)mat->data; 1485 1486 PetscFunctionBegin; 1487 if (shell->ctxcontainer) PetscCall(PetscContainerSetUserDestroy(shell->ctxcontainer, f)); 1488 PetscFunctionReturn(PETSC_SUCCESS); 1489 } 1490 1491 static PetscErrorCode MatShellSetVecType_Shell(Mat mat, VecType vtype) 1492 { 1493 PetscFunctionBegin; 1494 PetscCall(PetscFree(mat->defaultvectype)); 1495 PetscCall(PetscStrallocpy(vtype, (char **)&mat->defaultvectype)); 1496 PetscFunctionReturn(PETSC_SUCCESS); 1497 } 1498 1499 static PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A) 1500 { 1501 Mat_Shell *shell = (Mat_Shell *)A->data; 1502 1503 PetscFunctionBegin; 1504 shell->managescalingshifts = PETSC_FALSE; 1505 A->ops->diagonalset = NULL; 1506 A->ops->diagonalscale = NULL; 1507 A->ops->scale = NULL; 1508 A->ops->shift = NULL; 1509 A->ops->axpy = NULL; 1510 PetscFunctionReturn(PETSC_SUCCESS); 1511 } 1512 1513 static PetscErrorCode MatShellSetOperation_Shell(Mat mat, MatOperation op, void (*f)(void)) 1514 { 1515 Mat_Shell *shell = (Mat_Shell *)mat->data; 1516 1517 PetscFunctionBegin; 1518 switch (op) { 1519 case MATOP_DESTROY: 1520 shell->ops->destroy = (PetscErrorCode(*)(Mat))f; 1521 break; 1522 case MATOP_VIEW: 1523 if (!mat->ops->viewnative) mat->ops->viewnative = mat->ops->view; 1524 mat->ops->view = (PetscErrorCode(*)(Mat, PetscViewer))f; 1525 break; 1526 case MATOP_COPY: 1527 shell->ops->copy = (PetscErrorCode(*)(Mat, Mat, MatStructure))f; 1528 break; 1529 case MATOP_DIAGONAL_SET: 1530 case MATOP_DIAGONAL_SCALE: 1531 case MATOP_SHIFT: 1532 case MATOP_SCALE: 1533 case MATOP_AXPY: 1534 case MATOP_ZERO_ROWS: 1535 case MATOP_ZERO_ROWS_COLUMNS: 1536 PetscCheck(!shell->managescalingshifts, PetscObjectComm((PetscObject)mat), PETSC_ERR_ARG_WRONGSTATE, "MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1537 (((void (**)(void))mat->ops)[op]) = f; 1538 break; 1539 case MATOP_GET_DIAGONAL: 1540 if (shell->managescalingshifts) { 1541 shell->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f; 1542 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1543 } else { 1544 shell->ops->getdiagonal = NULL; 1545 mat->ops->getdiagonal = (PetscErrorCode(*)(Mat, Vec))f; 1546 } 1547 break; 1548 case MATOP_MULT: 1549 if (shell->managescalingshifts) { 1550 shell->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f; 1551 mat->ops->mult = MatMult_Shell; 1552 } else { 1553 shell->ops->mult = NULL; 1554 mat->ops->mult = (PetscErrorCode(*)(Mat, Vec, Vec))f; 1555 } 1556 break; 1557 case MATOP_MULT_TRANSPOSE: 1558 if (shell->managescalingshifts) { 1559 shell->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f; 1560 mat->ops->multtranspose = MatMultTranspose_Shell; 1561 } else { 1562 shell->ops->multtranspose = NULL; 1563 mat->ops->multtranspose = (PetscErrorCode(*)(Mat, Vec, Vec))f; 1564 } 1565 break; 1566 default: 1567 (((void (**)(void))mat->ops)[op]) = f; 1568 break; 1569 } 1570 PetscFunctionReturn(PETSC_SUCCESS); 1571 } 1572 1573 static PetscErrorCode MatShellGetOperation_Shell(Mat mat, MatOperation op, void (**f)(void)) 1574 { 1575 Mat_Shell *shell = (Mat_Shell *)mat->data; 1576 1577 PetscFunctionBegin; 1578 switch (op) { 1579 case MATOP_DESTROY: 1580 *f = (void (*)(void))shell->ops->destroy; 1581 break; 1582 case MATOP_VIEW: 1583 *f = (void (*)(void))mat->ops->view; 1584 break; 1585 case MATOP_COPY: 1586 *f = (void (*)(void))shell->ops->copy; 1587 break; 1588 case MATOP_DIAGONAL_SET: 1589 case MATOP_DIAGONAL_SCALE: 1590 case MATOP_SHIFT: 1591 case MATOP_SCALE: 1592 case MATOP_AXPY: 1593 case MATOP_ZERO_ROWS: 1594 case MATOP_ZERO_ROWS_COLUMNS: 1595 *f = (((void (**)(void))mat->ops)[op]); 1596 break; 1597 case MATOP_GET_DIAGONAL: 1598 if (shell->ops->getdiagonal) *f = (void (*)(void))shell->ops->getdiagonal; 1599 else *f = (((void (**)(void))mat->ops)[op]); 1600 break; 1601 case MATOP_MULT: 1602 if (shell->ops->mult) *f = (void (*)(void))shell->ops->mult; 1603 else *f = (((void (**)(void))mat->ops)[op]); 1604 break; 1605 case MATOP_MULT_TRANSPOSE: 1606 if (shell->ops->multtranspose) *f = (void (*)(void))shell->ops->multtranspose; 1607 else *f = (((void (**)(void))mat->ops)[op]); 1608 break; 1609 default: 1610 *f = (((void (**)(void))mat->ops)[op]); 1611 } 1612 PetscFunctionReturn(PETSC_SUCCESS); 1613 } 1614 1615 /*MC 1616 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix-free. 1617 1618 Level: advanced 1619 1620 .seealso: [](ch_matrices), `Mat`, `MatCreateShell()` 1621 M*/ 1622 1623 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 1624 { 1625 Mat_Shell *b; 1626 1627 PetscFunctionBegin; 1628 PetscCall(PetscNew(&b)); 1629 A->data = (void *)b; 1630 A->ops[0] = MatOps_Values; 1631 1632 b->ctxcontainer = NULL; 1633 b->vshift = 0.0; 1634 b->vscale = 1.0; 1635 b->managescalingshifts = PETSC_TRUE; 1636 A->assembled = PETSC_TRUE; 1637 A->preallocated = PETSC_FALSE; 1638 1639 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetContext_C", MatShellGetContext_Shell)); 1640 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContext_C", MatShellSetContext_Shell)); 1641 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Shell)); 1642 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetVecType_C", MatShellSetVecType_Shell)); 1643 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Shell)); 1644 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetOperation_C", MatShellSetOperation_Shell)); 1645 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellGetOperation_C", MatShellGetOperation_Shell)); 1646 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatShellSetMatProductOperation_C", MatShellSetMatProductOperation_Shell)); 1647 PetscCall(PetscObjectChangeTypeName((PetscObject)A, MATSHELL)); 1648 PetscFunctionReturn(PETSC_SUCCESS); 1649 } 1650 1651 /*@C 1652 MatCreateShell - Creates a new matrix of `MatType` `MATSHELL` for use with a user-defined 1653 private data storage format. 1654 1655 Collective 1656 1657 Input Parameters: 1658 + comm - MPI communicator 1659 . m - number of local rows (must be given) 1660 . n - number of local columns (must be given) 1661 . M - number of global rows (may be `PETSC_DETERMINE`) 1662 . N - number of global columns (may be `PETSC_DETERMINE`) 1663 - ctx - pointer to data needed by the shell matrix routines 1664 1665 Output Parameter: 1666 . A - the matrix 1667 1668 Level: advanced 1669 1670 Example Usage: 1671 .vb 1672 extern PetscErrorCode mult(Mat, Vec, Vec); 1673 1674 MatCreateShell(comm, m, n, M, N, ctx, &mat); 1675 MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))mult); 1676 // Use matrix for operations that have been set 1677 MatDestroy(mat); 1678 .ve 1679 1680 Notes: 1681 The shell matrix type is intended to provide a simple class to use 1682 with `KSP` (such as, for use with matrix-free methods). You should not 1683 use the shell type if you plan to define a complete matrix class. 1684 1685 PETSc requires that matrices and vectors being used for certain 1686 operations are partitioned accordingly. For example, when 1687 creating a shell matrix, `A`, that supports parallel matrix-vector 1688 products using `MatMult`(A,x,y) the user should set the number 1689 of local matrix rows to be the number of local elements of the 1690 corresponding result vector, y. Note that this is information is 1691 required for use of the matrix interface routines, even though 1692 the shell matrix may not actually be physically partitioned. 1693 For example, 1694 1695 .vb 1696 Vec x, y 1697 extern PetscErrorCode mult(Mat,Vec,Vec); 1698 Mat A 1699 1700 VecCreateMPI(comm,PETSC_DECIDE,M,&y); 1701 VecCreateMPI(comm,PETSC_DECIDE,N,&x); 1702 VecGetLocalSize(y,&m); 1703 VecGetLocalSize(x,&n); 1704 MatCreateShell(comm,m,n,M,N,ctx,&A); 1705 MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1706 MatMult(A,x,y); 1707 MatDestroy(&A); 1708 VecDestroy(&y); 1709 VecDestroy(&x); 1710 .ve 1711 1712 `MATSHELL` handles `MatShift()`, `MatDiagonalSet()`, `MatDiagonalScale()`, `MatAXPY()`, `MatScale()`, `MatZeroRows()` and `MatZeroRowsColumns()` 1713 internally, so these operations cannot be overwritten unless `MatShellSetManageScalingShifts()` is called. 1714 1715 Developer Notes: 1716 For rectangular matrices do all the scalings and shifts make sense? 1717 1718 Regarding shifting and scaling. The general form is 1719 1720 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1721 1722 The order you apply the operations is important. For example if you have a dshift then 1723 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 1724 you get s*vscale*A + diag(shift) 1725 1726 A is the user provided function. 1727 1728 `KSP`/`PC` uses changes in the `Mat`'s "state" to decide if preconditioners need to be rebuilt `PCSetUp()` only calls the setup() for 1729 for the `PC` implementation if the `Mat` state has increased from the previous call. Thus to get changes in a `MATSHELL` to trigger 1730 an update in the preconditioner you must call `MatAssemblyBegin()` and `MatAssemblyEnd()` or `PetscObjectStateIncrease`((`PetscObject`)mat); 1731 each time the `MATSHELL` matrix has changed. 1732 1733 Matrix product operations (i.e. `MatMat()`, `MatTransposeMat()` etc) can be specified using `MatShellSetMatProductOperation()` 1734 1735 Calling `MatAssemblyBegin()`/`MatAssemblyEnd()` on a `MATSHELL` removes any previously supplied shift and scales that were provided 1736 with `MatDiagonalSet()`, `MatShift()`, `MatScale()`, or `MatDiagonalScale()`. 1737 1738 Fortran Notes: 1739 To use this from Fortran with a `ctx` you must write an interface definition for this 1740 function and for `MatShellGetContext()` that tells Fortran the Fortran derived data type you are passing 1741 in as the `ctx` argument. 1742 1743 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()` 1744 @*/ 1745 PetscErrorCode MatCreateShell(MPI_Comm comm, PetscInt m, PetscInt n, PetscInt M, PetscInt N, void *ctx, Mat *A) 1746 { 1747 PetscFunctionBegin; 1748 PetscCall(MatCreate(comm, A)); 1749 PetscCall(MatSetSizes(*A, m, n, M, N)); 1750 PetscCall(MatSetType(*A, MATSHELL)); 1751 PetscCall(MatShellSetContext(*A, ctx)); 1752 PetscCall(MatSetUp(*A)); 1753 PetscFunctionReturn(PETSC_SUCCESS); 1754 } 1755 1756 /*@ 1757 MatShellSetContext - sets the context for a `MATSHELL` shell matrix 1758 1759 Logically Collective 1760 1761 Input Parameters: 1762 + mat - the `MATSHELL` shell matrix 1763 - ctx - the context 1764 1765 Level: advanced 1766 1767 Fortran Notes: 1768 You must write a Fortran interface definition for this 1769 function that tells Fortran the Fortran derived data type that you are passing in as the `ctx` argument. 1770 1771 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()` 1772 @*/ 1773 PetscErrorCode MatShellSetContext(Mat mat, void *ctx) 1774 { 1775 PetscFunctionBegin; 1776 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 1777 PetscTryMethod(mat, "MatShellSetContext_C", (Mat, void *), (mat, ctx)); 1778 PetscFunctionReturn(PETSC_SUCCESS); 1779 } 1780 1781 /*@C 1782 MatShellSetContextDestroy - sets the destroy function for a `MATSHELL` shell matrix context 1783 1784 Logically Collective 1785 1786 Input Parameters: 1787 + mat - the shell matrix 1788 - f - the context destroy function 1789 1790 Level: advanced 1791 1792 Note: 1793 If the `MatShell` is never duplicated, the behavior of this function is equivalent 1794 to `MatShellSetOperation`(`Mat`,`MATOP_DESTROY`,f). However, `MatShellSetContextDestroy()` 1795 ensures proper reference counting for the user provided context data in the case that 1796 the `MATSHELL` is duplicated. 1797 1798 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellSetContext()` 1799 @*/ 1800 PetscErrorCode MatShellSetContextDestroy(Mat mat, PetscErrorCode (*f)(void *)) 1801 { 1802 PetscFunctionBegin; 1803 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 1804 PetscTryMethod(mat, "MatShellSetContextDestroy_C", (Mat, PetscErrorCode(*)(void *)), (mat, f)); 1805 PetscFunctionReturn(PETSC_SUCCESS); 1806 } 1807 1808 /*@C 1809 MatShellSetVecType - Sets the `VecType` of `Vec` returned by `MatCreateVecs()` 1810 1811 Logically Collective 1812 1813 Input Parameters: 1814 + mat - the `MATSHELL` shell matrix 1815 - vtype - type to use for creating vectors 1816 1817 Level: advanced 1818 1819 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateVecs()` 1820 @*/ 1821 PetscErrorCode MatShellSetVecType(Mat mat, VecType vtype) 1822 { 1823 PetscFunctionBegin; 1824 PetscTryMethod(mat, "MatShellSetVecType_C", (Mat, VecType), (mat, vtype)); 1825 PetscFunctionReturn(PETSC_SUCCESS); 1826 } 1827 1828 /*@ 1829 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the `MATSHELL`. Must be called immediately 1830 after `MatCreateShell()` 1831 1832 Logically Collective 1833 1834 Input Parameter: 1835 . A - the `MATSHELL` shell matrix 1836 1837 Level: advanced 1838 1839 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()` 1840 @*/ 1841 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 1842 { 1843 PetscFunctionBegin; 1844 PetscValidHeaderSpecific(A, MAT_CLASSID, 1); 1845 PetscTryMethod(A, "MatShellSetManageScalingShifts_C", (Mat), (A)); 1846 PetscFunctionReturn(PETSC_SUCCESS); 1847 } 1848 1849 /*@C 1850 MatShellTestMult - Compares the multiply routine provided to the `MATSHELL` with differencing on a given function. 1851 1852 Logically Collective; No Fortran Support 1853 1854 Input Parameters: 1855 + mat - the `MATSHELL` shell matrix 1856 . f - the function 1857 . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated 1858 - ctx - an optional context for the function 1859 1860 Output Parameter: 1861 . flg - `PETSC_TRUE` if the multiply is likely correct 1862 1863 Options Database Key: 1864 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1865 1866 Level: advanced 1867 1868 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()` 1869 @*/ 1870 PetscErrorCode MatShellTestMult(Mat mat, PetscErrorCode (*f)(void *, Vec, Vec), Vec base, void *ctx, PetscBool *flg) 1871 { 1872 PetscInt m, n; 1873 Mat mf, Dmf, Dmat, Ddiff; 1874 PetscReal Diffnorm, Dmfnorm; 1875 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1876 1877 PetscFunctionBegin; 1878 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 1879 PetscCall(PetscOptionsHasName(NULL, ((PetscObject)mat)->prefix, "-mat_shell_test_mult_view", &v)); 1880 PetscCall(MatGetLocalSize(mat, &m, &n)); 1881 PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat), m, n, PETSC_DECIDE, PETSC_DECIDE, &mf)); 1882 PetscCall(MatMFFDSetFunction(mf, f, ctx)); 1883 PetscCall(MatMFFDSetBase(mf, base, NULL)); 1884 1885 PetscCall(MatComputeOperator(mf, MATAIJ, &Dmf)); 1886 PetscCall(MatComputeOperator(mat, MATAIJ, &Dmat)); 1887 1888 PetscCall(MatDuplicate(Dmat, MAT_COPY_VALUES, &Ddiff)); 1889 PetscCall(MatAXPY(Ddiff, -1.0, Dmf, DIFFERENT_NONZERO_PATTERN)); 1890 PetscCall(MatNorm(Ddiff, NORM_FROBENIUS, &Diffnorm)); 1891 PetscCall(MatNorm(Dmf, NORM_FROBENIUS, &Dmfnorm)); 1892 if (Diffnorm / Dmfnorm > 10 * PETSC_SQRT_MACHINE_EPSILON) { 1893 flag = PETSC_FALSE; 1894 if (v) { 1895 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))); 1896 PetscCall(MatViewFromOptions(Ddiff, (PetscObject)mat, "-mat_shell_test_mult_view")); 1897 PetscCall(MatViewFromOptions(Dmf, (PetscObject)mat, "-mat_shell_test_mult_view")); 1898 PetscCall(MatViewFromOptions(Dmat, (PetscObject)mat, "-mat_shell_test_mult_view")); 1899 } 1900 } else if (v) { 1901 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat), "MATSHELL and matrix-free multiple appear to produce the same results\n")); 1902 } 1903 if (flg) *flg = flag; 1904 PetscCall(MatDestroy(&Ddiff)); 1905 PetscCall(MatDestroy(&mf)); 1906 PetscCall(MatDestroy(&Dmf)); 1907 PetscCall(MatDestroy(&Dmat)); 1908 PetscFunctionReturn(PETSC_SUCCESS); 1909 } 1910 1911 /*@C 1912 MatShellTestMultTranspose - Compares the multiply transpose routine provided to the `MATSHELL` with differencing on a given function. 1913 1914 Logically Collective; No Fortran Support 1915 1916 Input Parameters: 1917 + mat - the `MATSHELL` shell matrix 1918 . f - the function 1919 . base - differences are computed around this vector, see `MatMFFDSetBase()`, for Jacobians this is the point at which the Jacobian is being evaluated 1920 - ctx - an optional context for the function 1921 1922 Output Parameter: 1923 . flg - `PETSC_TRUE` if the multiply is likely correct 1924 1925 Options Database Key: 1926 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1927 1928 Level: advanced 1929 1930 .seealso: [](ch_matrices), `Mat`, `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(PETSC_SUCCESS); 1978 } 1979 1980 /*@C 1981 MatShellSetOperation - Allows user to set a matrix operation for a `MATSHELL` shell matrix. 1982 1983 Logically Collective 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 Example Usage: 1993 .vb 1994 extern PetscErrorCode usermult(Mat, Vec, Vec); 1995 1996 MatCreateShell(comm, m, n, M, N, ctx, &A); 1997 MatShellSetOperation(A, MATOP_MULT, (void(*)(void))usermult); 1998 .ve 1999 2000 Notes: 2001 See the file include/petscmat.h for a complete list of matrix 2002 operations, which all have the form MATOP_<OPERATION>, where 2003 <OPERATION> is the name (in all capital letters) of the 2004 user interface routine (e.g., `MatMult()` -> `MATOP_MULT`). 2005 2006 All user-provided functions (except for `MATOP_DESTROY`) should have the same calling 2007 sequence as the usual matrix interface routines, since they 2008 are intended to be accessed via the usual matrix interface 2009 routines, e.g., 2010 $ MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec) 2011 2012 In particular each function MUST return an error code of 0 on success and 2013 nonzero on failure. 2014 2015 Within each user-defined routine, the user should call 2016 `MatShellGetContext()` to obtain the user-defined context that was 2017 set by `MatCreateShell()`. 2018 2019 Use `MatSetOperation()` to set an operation for any matrix type. For matrix product operations (i.e. `MatMatXXX()`, `MatTransposeMatXXX()` etc) 2020 use `MatShellSetMatProductOperation()` 2021 2022 Fortran Notes: 2023 For `MatCreateVecs()` the user code should check if the input left or right matrix is -1 and in that case not 2024 generate a matrix. See src/mat/tests/ex120f.F 2025 2026 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()` 2027 @*/ 2028 PetscErrorCode MatShellSetOperation(Mat mat, MatOperation op, void (*g)(void)) 2029 { 2030 PetscFunctionBegin; 2031 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2032 PetscTryMethod(mat, "MatShellSetOperation_C", (Mat, MatOperation, void (*)(void)), (mat, op, g)); 2033 PetscFunctionReturn(PETSC_SUCCESS); 2034 } 2035 2036 /*@C 2037 MatShellGetOperation - Gets a matrix function for a `MATSHELL` shell matrix. 2038 2039 Not Collective 2040 2041 Input Parameters: 2042 + mat - the `MATSHELL` shell matrix 2043 - op - the name of the operation 2044 2045 Output Parameter: 2046 . g - the function that provides the operation. 2047 2048 Level: advanced 2049 2050 Notes: 2051 See the file include/petscmat.h for a complete list of matrix 2052 operations, which all have the form MATOP_<OPERATION>, where 2053 <OPERATION> is the name (in all capital letters) of the 2054 user interface routine (e.g., `MatMult()` -> `MATOP_MULT`). 2055 2056 All user-provided functions have the same calling 2057 sequence as the usual matrix interface routines, since they 2058 are intended to be accessed via the usual matrix interface 2059 routines, e.g., 2060 $ MatMult(Mat, Vec, Vec) -> usermult(Mat, Vec, Vec) 2061 2062 Within each user-defined routine, the user should call 2063 `MatShellGetContext()` to obtain the user-defined context that was 2064 set by `MatCreateShell()`. 2065 2066 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()` 2067 @*/ 2068 PetscErrorCode MatShellGetOperation(Mat mat, MatOperation op, void (**g)(void)) 2069 { 2070 PetscFunctionBegin; 2071 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2072 PetscUseMethod(mat, "MatShellGetOperation_C", (Mat, MatOperation, void (**)(void)), (mat, op, g)); 2073 PetscFunctionReturn(PETSC_SUCCESS); 2074 } 2075 2076 /*@ 2077 MatIsShell - Inquires if a matrix is derived from `MATSHELL` 2078 2079 Input Parameter: 2080 . mat - the matrix 2081 2082 Output Parameter: 2083 . flg - the boolean value 2084 2085 Level: developer 2086 2087 Developer Notes: 2088 In the future, we should allow the object type name to be changed still using the `MATSHELL` data structure for other matrices 2089 (i.e. `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT` etc) 2090 2091 .seealso: [](ch_matrices), `Mat`, `MATSHELL`, `MATMFFD`, `MatCreateShell()`, `MATTRANSPOSEVIRTUAL`, `MATSCHURCOMPLEMENT` 2092 @*/ 2093 PetscErrorCode MatIsShell(Mat mat, PetscBool *flg) 2094 { 2095 PetscFunctionBegin; 2096 PetscValidHeaderSpecific(mat, MAT_CLASSID, 1); 2097 PetscAssertPointer(flg, 2); 2098 *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell); 2099 PetscFunctionReturn(PETSC_SUCCESS); 2100 } 2101