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