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 }; 1466 1467 PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx) 1468 { 1469 Mat_Shell *shell = (Mat_Shell*)mat->data; 1470 1471 PetscFunctionBegin; 1472 shell->ctx = ctx; 1473 PetscFunctionReturn(0); 1474 } 1475 1476 static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype) 1477 { 1478 PetscFunctionBegin; 1479 PetscCall(PetscFree(mat->defaultvectype)); 1480 PetscCall(PetscStrallocpy(vtype,(char**)&mat->defaultvectype)); 1481 PetscFunctionReturn(0); 1482 } 1483 1484 PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A) 1485 { 1486 Mat_Shell *shell = (Mat_Shell*)A->data; 1487 1488 PetscFunctionBegin; 1489 shell->managescalingshifts = PETSC_FALSE; 1490 A->ops->diagonalset = NULL; 1491 A->ops->diagonalscale = NULL; 1492 A->ops->scale = NULL; 1493 A->ops->shift = NULL; 1494 A->ops->axpy = NULL; 1495 PetscFunctionReturn(0); 1496 } 1497 1498 PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void)) 1499 { 1500 Mat_Shell *shell = (Mat_Shell*)mat->data; 1501 1502 PetscFunctionBegin; 1503 switch (op) { 1504 case MATOP_DESTROY: 1505 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1506 break; 1507 case MATOP_VIEW: 1508 if (!mat->ops->viewnative) { 1509 mat->ops->viewnative = mat->ops->view; 1510 } 1511 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1512 break; 1513 case MATOP_COPY: 1514 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1515 break; 1516 case MATOP_DIAGONAL_SET: 1517 case MATOP_DIAGONAL_SCALE: 1518 case MATOP_SHIFT: 1519 case MATOP_SCALE: 1520 case MATOP_AXPY: 1521 case MATOP_ZERO_ROWS: 1522 case MATOP_ZERO_ROWS_COLUMNS: 1523 PetscCheck(!shell->managescalingshifts,PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1524 (((void(**)(void))mat->ops)[op]) = f; 1525 break; 1526 case MATOP_GET_DIAGONAL: 1527 if (shell->managescalingshifts) { 1528 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1529 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1530 } else { 1531 shell->ops->getdiagonal = NULL; 1532 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1533 } 1534 break; 1535 case MATOP_MULT: 1536 if (shell->managescalingshifts) { 1537 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1538 mat->ops->mult = MatMult_Shell; 1539 } else { 1540 shell->ops->mult = NULL; 1541 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1542 } 1543 break; 1544 case MATOP_MULT_TRANSPOSE: 1545 if (shell->managescalingshifts) { 1546 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1547 mat->ops->multtranspose = MatMultTranspose_Shell; 1548 } else { 1549 shell->ops->multtranspose = NULL; 1550 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1551 } 1552 break; 1553 default: 1554 (((void(**)(void))mat->ops)[op]) = f; 1555 break; 1556 } 1557 PetscFunctionReturn(0); 1558 } 1559 1560 PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void)) 1561 { 1562 Mat_Shell *shell = (Mat_Shell*)mat->data; 1563 1564 PetscFunctionBegin; 1565 switch (op) { 1566 case MATOP_DESTROY: 1567 *f = (void (*)(void))shell->ops->destroy; 1568 break; 1569 case MATOP_VIEW: 1570 *f = (void (*)(void))mat->ops->view; 1571 break; 1572 case MATOP_COPY: 1573 *f = (void (*)(void))shell->ops->copy; 1574 break; 1575 case MATOP_DIAGONAL_SET: 1576 case MATOP_DIAGONAL_SCALE: 1577 case MATOP_SHIFT: 1578 case MATOP_SCALE: 1579 case MATOP_AXPY: 1580 case MATOP_ZERO_ROWS: 1581 case MATOP_ZERO_ROWS_COLUMNS: 1582 *f = (((void (**)(void))mat->ops)[op]); 1583 break; 1584 case MATOP_GET_DIAGONAL: 1585 if (shell->ops->getdiagonal) 1586 *f = (void (*)(void))shell->ops->getdiagonal; 1587 else 1588 *f = (((void (**)(void))mat->ops)[op]); 1589 break; 1590 case MATOP_MULT: 1591 if (shell->ops->mult) 1592 *f = (void (*)(void))shell->ops->mult; 1593 else 1594 *f = (((void (**)(void))mat->ops)[op]); 1595 break; 1596 case MATOP_MULT_TRANSPOSE: 1597 if (shell->ops->multtranspose) 1598 *f = (void (*)(void))shell->ops->multtranspose; 1599 else 1600 *f = (((void (**)(void))mat->ops)[op]); 1601 break; 1602 default: 1603 *f = (((void (**)(void))mat->ops)[op]); 1604 } 1605 PetscFunctionReturn(0); 1606 } 1607 1608 /*MC 1609 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 1610 1611 Level: advanced 1612 1613 .seealso: `MatCreateShell()` 1614 M*/ 1615 1616 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 1617 { 1618 Mat_Shell *b; 1619 1620 PetscFunctionBegin; 1621 PetscCall(PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps))); 1622 1623 PetscCall(PetscNewLog(A,&b)); 1624 A->data = (void*)b; 1625 1626 b->ctx = NULL; 1627 b->vshift = 0.0; 1628 b->vscale = 1.0; 1629 b->managescalingshifts = PETSC_TRUE; 1630 A->assembled = PETSC_TRUE; 1631 A->preallocated = PETSC_FALSE; 1632 1633 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell)); 1634 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell)); 1635 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell)); 1636 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell)); 1637 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell)); 1638 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell)); 1639 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetMatProductOperation_C",MatShellSetMatProductOperation_Shell)); 1640 PetscCall(PetscObjectChangeTypeName((PetscObject)A,MATSHELL)); 1641 PetscFunctionReturn(0); 1642 } 1643 1644 /*@C 1645 MatCreateShell - Creates a new matrix class for use with a user-defined 1646 private data storage format. 1647 1648 Collective 1649 1650 Input Parameters: 1651 + comm - MPI communicator 1652 . m - number of local rows (must be given) 1653 . n - number of local columns (must be given) 1654 . M - number of global rows (may be PETSC_DETERMINE) 1655 . N - number of global columns (may be PETSC_DETERMINE) 1656 - ctx - pointer to data needed by the shell matrix routines 1657 1658 Output Parameter: 1659 . A - the matrix 1660 1661 Level: advanced 1662 1663 Usage: 1664 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1665 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 1666 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1667 $ [ Use matrix for operations that have been set ] 1668 $ MatDestroy(mat); 1669 1670 Notes: 1671 The shell matrix type is intended to provide a simple class to use 1672 with KSP (such as, for use with matrix-free methods). You should not 1673 use the shell type if you plan to define a complete matrix class. 1674 1675 Fortran Notes: 1676 To use this from Fortran with a ctx you must write an interface definition for this 1677 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 1678 in as the ctx argument. 1679 1680 PETSc requires that matrices and vectors being used for certain 1681 operations are partitioned accordingly. For example, when 1682 creating a shell matrix, A, that supports parallel matrix-vector 1683 products using MatMult(A,x,y) the user should set the number 1684 of local matrix rows to be the number of local elements of the 1685 corresponding result vector, y. Note that this is information is 1686 required for use of the matrix interface routines, even though 1687 the shell matrix may not actually be physically partitioned. 1688 For example, 1689 1690 $ 1691 $ Vec x, y 1692 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1693 $ Mat A 1694 $ 1695 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 1696 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 1697 $ VecGetLocalSize(y,&m); 1698 $ VecGetLocalSize(x,&n); 1699 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 1700 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1701 $ MatMult(A,x,y); 1702 $ MatDestroy(&A); 1703 $ VecDestroy(&y); 1704 $ VecDestroy(&x); 1705 $ 1706 1707 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these 1708 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 1709 1710 For rectangular matrices do all the scalings and shifts make sense? 1711 1712 Developers Notes: 1713 Regarding shifting and scaling. The general form is 1714 1715 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1716 1717 The order you apply the operations is important. For example if you have a dshift then 1718 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 1719 you get s*vscale*A + diag(shift) 1720 1721 A is the user provided function. 1722 1723 KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for 1724 for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger 1725 an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat); 1726 each time the MATSHELL matrix has changed. 1727 1728 Matrix product operations (i.e. MatMat, MatTransposeMat etc) can be specified using MatShellSetMatProductOperation() 1729 1730 Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided 1731 with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale(). 1732 1733 .seealso: `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MATSHELL`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()` 1734 @*/ 1735 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 1736 { 1737 PetscFunctionBegin; 1738 PetscCall(MatCreate(comm,A)); 1739 PetscCall(MatSetSizes(*A,m,n,M,N)); 1740 PetscCall(MatSetType(*A,MATSHELL)); 1741 PetscCall(MatShellSetContext(*A,ctx)); 1742 PetscCall(MatSetUp(*A)); 1743 PetscFunctionReturn(0); 1744 } 1745 1746 /*@ 1747 MatShellSetContext - sets the context for a shell matrix 1748 1749 Logically Collective on Mat 1750 1751 Input Parameters: 1752 + mat - the shell matrix 1753 - ctx - the context 1754 1755 Level: advanced 1756 1757 Fortran Notes: 1758 To use this from Fortran you must write a Fortran interface definition for this 1759 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 1760 1761 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()` 1762 @*/ 1763 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 1764 { 1765 PetscFunctionBegin; 1766 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1767 PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx)); 1768 PetscFunctionReturn(0); 1769 } 1770 1771 /*@C 1772 MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs() 1773 1774 Logically collective 1775 1776 Input Parameters: 1777 + mat - the shell matrix 1778 - vtype - type to use for creating vectors 1779 1780 Notes: 1781 1782 Level: advanced 1783 1784 .seealso: `MatCreateVecs()` 1785 @*/ 1786 PetscErrorCode MatShellSetVecType(Mat mat,VecType vtype) 1787 { 1788 PetscFunctionBegin; 1789 PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype)); 1790 PetscFunctionReturn(0); 1791 } 1792 1793 /*@ 1794 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 1795 after MatCreateShell() 1796 1797 Logically Collective on Mat 1798 1799 Input Parameter: 1800 . mat - the shell matrix 1801 1802 Level: advanced 1803 1804 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()` 1805 @*/ 1806 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 1807 { 1808 PetscFunctionBegin; 1809 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1810 PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A)); 1811 PetscFunctionReturn(0); 1812 } 1813 1814 /*@C 1815 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 1816 1817 Logically Collective on Mat 1818 1819 Input Parameters: 1820 + mat - the shell matrix 1821 . f - the function 1822 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1823 - ctx - an optional context for the function 1824 1825 Output Parameter: 1826 . flg - PETSC_TRUE if the multiply is likely correct 1827 1828 Options Database: 1829 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1830 1831 Level: advanced 1832 1833 Fortran Notes: 1834 Not supported from Fortran 1835 1836 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()` 1837 @*/ 1838 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1839 { 1840 PetscInt m,n; 1841 Mat mf,Dmf,Dmat,Ddiff; 1842 PetscReal Diffnorm,Dmfnorm; 1843 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1844 1845 PetscFunctionBegin; 1846 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1847 PetscCall(PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v)); 1848 PetscCall(MatGetLocalSize(mat,&m,&n)); 1849 PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf)); 1850 PetscCall(MatMFFDSetFunction(mf,f,ctx)); 1851 PetscCall(MatMFFDSetBase(mf,base,NULL)); 1852 1853 PetscCall(MatComputeOperator(mf,MATAIJ,&Dmf)); 1854 PetscCall(MatComputeOperator(mat,MATAIJ,&Dmat)); 1855 1856 PetscCall(MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff)); 1857 PetscCall(MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN)); 1858 PetscCall(MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm)); 1859 PetscCall(MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm)); 1860 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 1861 flag = PETSC_FALSE; 1862 if (v) { 1863 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))); 1864 PetscCall(MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view")); 1865 PetscCall(MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view")); 1866 PetscCall(MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view")); 1867 } 1868 } else if (v) { 1869 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n")); 1870 } 1871 if (flg) *flg = flag; 1872 PetscCall(MatDestroy(&Ddiff)); 1873 PetscCall(MatDestroy(&mf)); 1874 PetscCall(MatDestroy(&Dmf)); 1875 PetscCall(MatDestroy(&Dmat)); 1876 PetscFunctionReturn(0); 1877 } 1878 1879 /*@C 1880 MatShellTestMultTranspose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 1881 1882 Logically Collective on Mat 1883 1884 Input Parameters: 1885 + mat - the shell matrix 1886 . f - the function 1887 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1888 - ctx - an optional context for the function 1889 1890 Output Parameter: 1891 . flg - PETSC_TRUE if the multiply is likely correct 1892 1893 Options Database: 1894 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1895 1896 Level: advanced 1897 1898 Fortran Notes: 1899 Not supported from Fortran 1900 1901 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()` 1902 @*/ 1903 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1904 { 1905 Vec x,y,z; 1906 PetscInt m,n,M,N; 1907 Mat mf,Dmf,Dmat,Ddiff; 1908 PetscReal Diffnorm,Dmfnorm; 1909 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1910 1911 PetscFunctionBegin; 1912 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1913 PetscCall(PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v)); 1914 PetscCall(MatCreateVecs(mat,&x,&y)); 1915 PetscCall(VecDuplicate(y,&z)); 1916 PetscCall(MatGetLocalSize(mat,&m,&n)); 1917 PetscCall(MatGetSize(mat,&M,&N)); 1918 PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf)); 1919 PetscCall(MatMFFDSetFunction(mf,f,ctx)); 1920 PetscCall(MatMFFDSetBase(mf,base,NULL)); 1921 PetscCall(MatComputeOperator(mf,MATAIJ,&Dmf)); 1922 PetscCall(MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf)); 1923 PetscCall(MatComputeOperatorTranspose(mat,MATAIJ,&Dmat)); 1924 1925 PetscCall(MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff)); 1926 PetscCall(MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN)); 1927 PetscCall(MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm)); 1928 PetscCall(MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm)); 1929 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 1930 flag = PETSC_FALSE; 1931 if (v) { 1932 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))); 1933 PetscCall(MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view")); 1934 PetscCall(MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view")); 1935 PetscCall(MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view")); 1936 } 1937 } else if (v) { 1938 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n")); 1939 } 1940 if (flg) *flg = flag; 1941 PetscCall(MatDestroy(&mf)); 1942 PetscCall(MatDestroy(&Dmat)); 1943 PetscCall(MatDestroy(&Ddiff)); 1944 PetscCall(MatDestroy(&Dmf)); 1945 PetscCall(VecDestroy(&x)); 1946 PetscCall(VecDestroy(&y)); 1947 PetscCall(VecDestroy(&z)); 1948 PetscFunctionReturn(0); 1949 } 1950 1951 /*@C 1952 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 1953 1954 Logically Collective on Mat 1955 1956 Input Parameters: 1957 + mat - the shell matrix 1958 . op - the name of the operation 1959 - g - the function that provides the operation. 1960 1961 Level: advanced 1962 1963 Usage: 1964 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 1965 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 1966 $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 1967 1968 Notes: 1969 See the file include/petscmat.h for a complete list of matrix 1970 operations, which all have the form MATOP_<OPERATION>, where 1971 <OPERATION> is the name (in all capital letters) of the 1972 user interface routine (e.g., MatMult() -> MATOP_MULT). 1973 1974 All user-provided functions (except for MATOP_DESTROY) should have the same calling 1975 sequence as the usual matrix interface routines, since they 1976 are intended to be accessed via the usual matrix interface 1977 routines, e.g., 1978 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1979 1980 In particular each function MUST return an error code of 0 on success and 1981 nonzero on failure. 1982 1983 Within each user-defined routine, the user should call 1984 MatShellGetContext() to obtain the user-defined context that was 1985 set by MatCreateShell(). 1986 1987 Use MatSetOperation() to set an operation for any matrix type. For matrix product operations (i.e. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation() 1988 1989 Fortran Notes: 1990 For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 1991 generate a matrix. See src/mat/tests/ex120f.F 1992 1993 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()` 1994 @*/ 1995 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void)) 1996 { 1997 PetscFunctionBegin; 1998 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1999 PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g)); 2000 PetscFunctionReturn(0); 2001 } 2002 2003 /*@C 2004 MatShellGetOperation - Gets a matrix function for a shell matrix. 2005 2006 Not Collective 2007 2008 Input Parameters: 2009 + mat - the shell matrix 2010 - op - the name of the operation 2011 2012 Output Parameter: 2013 . g - the function that provides the operation. 2014 2015 Level: advanced 2016 2017 Notes: 2018 See the file include/petscmat.h for a complete list of matrix 2019 operations, which all have the form MATOP_<OPERATION>, where 2020 <OPERATION> is the name (in all capital letters) of the 2021 user interface routine (e.g., MatMult() -> MATOP_MULT). 2022 2023 All user-provided functions have the same calling 2024 sequence as the usual matrix interface routines, since they 2025 are intended to be accessed via the usual matrix interface 2026 routines, e.g., 2027 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 2028 2029 Within each user-defined routine, the user should call 2030 MatShellGetContext() to obtain the user-defined context that was 2031 set by MatCreateShell(). 2032 2033 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()` 2034 @*/ 2035 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void)) 2036 { 2037 PetscFunctionBegin; 2038 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2039 PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g)); 2040 PetscFunctionReturn(0); 2041 } 2042 2043 /*@ 2044 MatIsShell - Inquires if a matrix is derived from MATSHELL 2045 2046 Input Parameter: 2047 . mat - the matrix 2048 2049 Output Parameter: 2050 . flg - the boolean value 2051 2052 Level: developer 2053 2054 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) 2055 2056 .seealso: `MatCreateShell()` 2057 @*/ 2058 PetscErrorCode MatIsShell(Mat mat, PetscBool *flg) 2059 { 2060 PetscFunctionBegin; 2061 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2062 PetscValidBoolPointer(flg,2); 2063 *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell); 2064 PetscFunctionReturn(0); 2065 } 2066