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