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