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