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