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