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