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 PetscCheck(shell->ops->mult,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL"); 1019 PetscCall(MatShellPreZeroRight(A,x,&xx)); 1020 PetscCall(MatShellPreScaleRight(A,xx,&xx)); 1021 PetscCall(PetscObjectStateGet((PetscObject)y, &instate)); 1022 PetscCall((*shell->ops->mult)(A,xx,y)); 1023 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate)); 1024 if (instate == outstate) { 1025 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1026 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 1027 } 1028 PetscCall(MatShellShiftAndScale(A,xx,y)); 1029 PetscCall(MatShellPostScaleLeft(A,y)); 1030 PetscCall(MatShellPostZeroLeft(A,y)); 1031 1032 if (shell->axpy) { 1033 Mat X; 1034 PetscObjectState axpy_state; 1035 1036 PetscCall(MatShellGetContext(shell->axpy,&X)); 1037 PetscCall(PetscObjectStateGet((PetscObject)X,&axpy_state)); 1038 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,...)"); 1039 1040 PetscCall(MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left)); 1041 PetscCall(VecCopy(x,shell->axpy_right)); 1042 PetscCall(MatMult(shell->axpy,shell->axpy_right,shell->axpy_left)); 1043 PetscCall(VecAXPY(y,shell->axpy_vscale,shell->axpy_left)); 1044 } 1045 PetscFunctionReturn(0); 1046 } 1047 1048 static PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z) 1049 { 1050 Mat_Shell *shell = (Mat_Shell*)A->data; 1051 1052 PetscFunctionBegin; 1053 if (y == z) { 1054 if (!shell->right_add_work) PetscCall(VecDuplicate(z,&shell->right_add_work)); 1055 PetscCall(MatMult(A,x,shell->right_add_work)); 1056 PetscCall(VecAXPY(z,1.0,shell->right_add_work)); 1057 } else { 1058 PetscCall(MatMult(A,x,z)); 1059 PetscCall(VecAXPY(z,1.0,y)); 1060 } 1061 PetscFunctionReturn(0); 1062 } 1063 1064 static PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y) 1065 { 1066 Mat_Shell *shell = (Mat_Shell*)A->data; 1067 Vec xx; 1068 PetscObjectState instate,outstate; 1069 1070 PetscFunctionBegin; 1071 PetscCheck(shell->ops->multtranspose,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL"); 1072 PetscCall(MatShellPreZeroLeft(A,x,&xx)); 1073 PetscCall(MatShellPreScaleLeft(A,xx,&xx)); 1074 PetscCall(PetscObjectStateGet((PetscObject)y, &instate)); 1075 PetscCall((*shell->ops->multtranspose)(A,xx,y)); 1076 PetscCall(PetscObjectStateGet((PetscObject)y, &outstate)); 1077 if (instate == outstate) { 1078 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 1079 PetscCall(PetscObjectStateIncrease((PetscObject)y)); 1080 } 1081 PetscCall(MatShellShiftAndScale(A,xx,y)); 1082 PetscCall(MatShellPostScaleRight(A,y)); 1083 PetscCall(MatShellPostZeroRight(A,y)); 1084 1085 if (shell->axpy) { 1086 Mat X; 1087 PetscObjectState axpy_state; 1088 1089 PetscCall(MatShellGetContext(shell->axpy,&X)); 1090 PetscCall(PetscObjectStateGet((PetscObject)X,&axpy_state)); 1091 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,...)"); 1092 PetscCall(MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left)); 1093 PetscCall(VecCopy(x,shell->axpy_left)); 1094 PetscCall(MatMultTranspose(shell->axpy,shell->axpy_left,shell->axpy_right)); 1095 PetscCall(VecAXPY(y,shell->axpy_vscale,shell->axpy_right)); 1096 } 1097 PetscFunctionReturn(0); 1098 } 1099 1100 static PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z) 1101 { 1102 Mat_Shell *shell = (Mat_Shell*)A->data; 1103 1104 PetscFunctionBegin; 1105 if (y == z) { 1106 if (!shell->left_add_work) PetscCall(VecDuplicate(z,&shell->left_add_work)); 1107 PetscCall(MatMultTranspose(A,x,shell->left_add_work)); 1108 PetscCall(VecAXPY(z,1.0,shell->left_add_work)); 1109 } else { 1110 PetscCall(MatMultTranspose(A,x,z)); 1111 PetscCall(VecAXPY(z,1.0,y)); 1112 } 1113 PetscFunctionReturn(0); 1114 } 1115 1116 /* 1117 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1118 */ 1119 static PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v) 1120 { 1121 Mat_Shell *shell = (Mat_Shell*)A->data; 1122 1123 PetscFunctionBegin; 1124 if (shell->ops->getdiagonal) { 1125 PetscCall((*shell->ops->getdiagonal)(A,v)); 1126 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)"); 1127 PetscCall(VecScale(v,shell->vscale)); 1128 if (shell->dshift) PetscCall(VecAXPY(v,1.0,shell->dshift)); 1129 PetscCall(VecShift(v,shell->vshift)); 1130 if (shell->left) PetscCall(VecPointwiseMult(v,v,shell->left)); 1131 if (shell->right) PetscCall(VecPointwiseMult(v,v,shell->right)); 1132 if (shell->zrows) { 1133 PetscCall(VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE)); 1134 PetscCall(VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE)); 1135 } 1136 if (shell->axpy) { 1137 Mat X; 1138 PetscObjectState axpy_state; 1139 1140 PetscCall(MatShellGetContext(shell->axpy,&X)); 1141 PetscCall(PetscObjectStateGet((PetscObject)X,&axpy_state)); 1142 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,...)"); 1143 PetscCall(MatCreateVecs(shell->axpy,NULL,shell->axpy_left ? NULL : &shell->axpy_left)); 1144 PetscCall(MatGetDiagonal(shell->axpy,shell->axpy_left)); 1145 PetscCall(VecAXPY(v,shell->axpy_vscale,shell->axpy_left)); 1146 } 1147 PetscFunctionReturn(0); 1148 } 1149 1150 static PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 1151 { 1152 Mat_Shell *shell = (Mat_Shell*)Y->data; 1153 PetscBool flg; 1154 1155 PetscFunctionBegin; 1156 PetscCall(MatHasCongruentLayouts(Y,&flg)); 1157 PetscCheck(flg,PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent"); 1158 if (shell->left || shell->right) { 1159 if (!shell->dshift) { 1160 PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift)); 1161 PetscCall(VecSet(shell->dshift,a)); 1162 } else { 1163 if (shell->left) PetscCall(VecPointwiseMult(shell->dshift,shell->dshift,shell->left)); 1164 if (shell->right) PetscCall(VecPointwiseMult(shell->dshift,shell->dshift,shell->right)); 1165 PetscCall(VecShift(shell->dshift,a)); 1166 } 1167 if (shell->left) PetscCall(VecPointwiseDivide(shell->dshift,shell->dshift,shell->left)); 1168 if (shell->right) PetscCall(VecPointwiseDivide(shell->dshift,shell->dshift,shell->right)); 1169 } else shell->vshift += a; 1170 if (shell->zrows) PetscCall(VecShift(shell->zvals,a)); 1171 PetscFunctionReturn(0); 1172 } 1173 1174 static PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s) 1175 { 1176 Mat_Shell *shell = (Mat_Shell*)A->data; 1177 1178 PetscFunctionBegin; 1179 if (!shell->dshift) PetscCall(VecDuplicate(D,&shell->dshift)); 1180 if (shell->left || shell->right) { 1181 if (!shell->right_work) PetscCall(VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work)); 1182 if (shell->left && shell->right) { 1183 PetscCall(VecPointwiseDivide(shell->right_work,D,shell->left)); 1184 PetscCall(VecPointwiseDivide(shell->right_work,shell->right_work,shell->right)); 1185 } else if (shell->left) { 1186 PetscCall(VecPointwiseDivide(shell->right_work,D,shell->left)); 1187 } else { 1188 PetscCall(VecPointwiseDivide(shell->right_work,D,shell->right)); 1189 } 1190 PetscCall(VecAXPY(shell->dshift,s,shell->right_work)); 1191 } else { 1192 PetscCall(VecAXPY(shell->dshift,s,D)); 1193 } 1194 PetscFunctionReturn(0); 1195 } 1196 1197 PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins) 1198 { 1199 Mat_Shell *shell = (Mat_Shell*)A->data; 1200 Vec d; 1201 PetscBool flg; 1202 1203 PetscFunctionBegin; 1204 PetscCall(MatHasCongruentLayouts(A,&flg)); 1205 PetscCheck(flg,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent"); 1206 if (ins == INSERT_VALUES) { 1207 PetscCheck(A->ops->getdiagonal,PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first"); 1208 PetscCall(VecDuplicate(D,&d)); 1209 PetscCall(MatGetDiagonal(A,d)); 1210 PetscCall(MatDiagonalSet_Shell_Private(A,d,-1.)); 1211 PetscCall(MatDiagonalSet_Shell_Private(A,D,1.)); 1212 PetscCall(VecDestroy(&d)); 1213 if (shell->zrows) PetscCall(VecCopy(D,shell->zvals)); 1214 } else { 1215 PetscCall(MatDiagonalSet_Shell_Private(A,D,1.)); 1216 if (shell->zrows) PetscCall(VecAXPY(shell->zvals,1.0,D)); 1217 } 1218 PetscFunctionReturn(0); 1219 } 1220 1221 static PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 1222 { 1223 Mat_Shell *shell = (Mat_Shell*)Y->data; 1224 1225 PetscFunctionBegin; 1226 shell->vscale *= a; 1227 shell->vshift *= a; 1228 if (shell->dshift) PetscCall(VecScale(shell->dshift,a)); 1229 shell->axpy_vscale *= a; 1230 if (shell->zrows) PetscCall(VecScale(shell->zvals,a)); 1231 PetscFunctionReturn(0); 1232 } 1233 1234 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right) 1235 { 1236 Mat_Shell *shell = (Mat_Shell*)Y->data; 1237 1238 PetscFunctionBegin; 1239 if (left) { 1240 if (!shell->left) { 1241 PetscCall(VecDuplicate(left,&shell->left)); 1242 PetscCall(VecCopy(left,shell->left)); 1243 } else { 1244 PetscCall(VecPointwiseMult(shell->left,shell->left,left)); 1245 } 1246 if (shell->zrows) PetscCall(VecPointwiseMult(shell->zvals,shell->zvals,left)); 1247 } 1248 if (right) { 1249 if (!shell->right) { 1250 PetscCall(VecDuplicate(right,&shell->right)); 1251 PetscCall(VecCopy(right,shell->right)); 1252 } else { 1253 PetscCall(VecPointwiseMult(shell->right,shell->right,right)); 1254 } 1255 if (shell->zrows) { 1256 if (!shell->left_work) { 1257 PetscCall(MatCreateVecs(Y,NULL,&shell->left_work)); 1258 } 1259 PetscCall(VecSet(shell->zvals_w,1.0)); 1260 PetscCall(VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD)); 1261 PetscCall(VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD)); 1262 PetscCall(VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w)); 1263 } 1264 } 1265 if (shell->axpy) PetscCall(MatDiagonalScale(shell->axpy,left,right)); 1266 PetscFunctionReturn(0); 1267 } 1268 1269 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 1270 { 1271 Mat_Shell *shell = (Mat_Shell*)Y->data; 1272 1273 PetscFunctionBegin; 1274 if (t == MAT_FINAL_ASSEMBLY) { 1275 shell->vshift = 0.0; 1276 shell->vscale = 1.0; 1277 shell->axpy_vscale = 0.0; 1278 shell->axpy_state = 0; 1279 PetscCall(VecDestroy(&shell->dshift)); 1280 PetscCall(VecDestroy(&shell->left)); 1281 PetscCall(VecDestroy(&shell->right)); 1282 PetscCall(MatDestroy(&shell->axpy)); 1283 PetscCall(VecDestroy(&shell->axpy_left)); 1284 PetscCall(VecDestroy(&shell->axpy_right)); 1285 PetscCall(VecScatterDestroy(&shell->zvals_sct_c)); 1286 PetscCall(VecScatterDestroy(&shell->zvals_sct_r)); 1287 PetscCall(ISDestroy(&shell->zrows)); 1288 PetscCall(ISDestroy(&shell->zcols)); 1289 } 1290 PetscFunctionReturn(0); 1291 } 1292 1293 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d) 1294 { 1295 PetscFunctionBegin; 1296 *missing = PETSC_FALSE; 1297 PetscFunctionReturn(0); 1298 } 1299 1300 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str) 1301 { 1302 Mat_Shell *shell = (Mat_Shell*)Y->data; 1303 1304 PetscFunctionBegin; 1305 if (X == Y) { 1306 PetscCall(MatScale(Y,1.0 + a)); 1307 PetscFunctionReturn(0); 1308 } 1309 if (!shell->axpy) { 1310 PetscCall(MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy)); 1311 shell->axpy_vscale = a; 1312 PetscCall(PetscObjectStateGet((PetscObject)X,&shell->axpy_state)); 1313 } else { 1314 PetscCall(MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str)); 1315 } 1316 PetscFunctionReturn(0); 1317 } 1318 1319 static struct _MatOps MatOps_Values = {NULL, 1320 NULL, 1321 NULL, 1322 NULL, 1323 /* 4*/ MatMultAdd_Shell, 1324 NULL, 1325 MatMultTransposeAdd_Shell, 1326 NULL, 1327 NULL, 1328 NULL, 1329 /*10*/ NULL, 1330 NULL, 1331 NULL, 1332 NULL, 1333 NULL, 1334 /*15*/ NULL, 1335 NULL, 1336 NULL, 1337 MatDiagonalScale_Shell, 1338 NULL, 1339 /*20*/ NULL, 1340 MatAssemblyEnd_Shell, 1341 NULL, 1342 NULL, 1343 /*24*/ MatZeroRows_Shell, 1344 NULL, 1345 NULL, 1346 NULL, 1347 NULL, 1348 /*29*/ NULL, 1349 NULL, 1350 NULL, 1351 NULL, 1352 NULL, 1353 /*34*/ MatDuplicate_Shell, 1354 NULL, 1355 NULL, 1356 NULL, 1357 NULL, 1358 /*39*/ MatAXPY_Shell, 1359 NULL, 1360 NULL, 1361 NULL, 1362 MatCopy_Shell, 1363 /*44*/ NULL, 1364 MatScale_Shell, 1365 MatShift_Shell, 1366 MatDiagonalSet_Shell, 1367 MatZeroRowsColumns_Shell, 1368 /*49*/ NULL, 1369 NULL, 1370 NULL, 1371 NULL, 1372 NULL, 1373 /*54*/ NULL, 1374 NULL, 1375 NULL, 1376 NULL, 1377 NULL, 1378 /*59*/ NULL, 1379 MatDestroy_Shell, 1380 NULL, 1381 MatConvertFrom_Shell, 1382 NULL, 1383 /*64*/ NULL, 1384 NULL, 1385 NULL, 1386 NULL, 1387 NULL, 1388 /*69*/ NULL, 1389 NULL, 1390 MatConvert_Shell, 1391 NULL, 1392 NULL, 1393 /*74*/ NULL, 1394 NULL, 1395 NULL, 1396 NULL, 1397 NULL, 1398 /*79*/ NULL, 1399 NULL, 1400 NULL, 1401 NULL, 1402 NULL, 1403 /*84*/ NULL, 1404 NULL, 1405 NULL, 1406 NULL, 1407 NULL, 1408 /*89*/ NULL, 1409 NULL, 1410 NULL, 1411 NULL, 1412 NULL, 1413 /*94*/ NULL, 1414 NULL, 1415 NULL, 1416 NULL, 1417 NULL, 1418 /*99*/ NULL, 1419 NULL, 1420 NULL, 1421 NULL, 1422 NULL, 1423 /*104*/ NULL, 1424 NULL, 1425 NULL, 1426 NULL, 1427 NULL, 1428 /*109*/ NULL, 1429 NULL, 1430 NULL, 1431 NULL, 1432 MatMissingDiagonal_Shell, 1433 /*114*/ NULL, 1434 NULL, 1435 NULL, 1436 NULL, 1437 NULL, 1438 /*119*/ NULL, 1439 NULL, 1440 NULL, 1441 NULL, 1442 NULL, 1443 /*124*/ NULL, 1444 NULL, 1445 NULL, 1446 NULL, 1447 NULL, 1448 /*129*/ NULL, 1449 NULL, 1450 NULL, 1451 NULL, 1452 NULL, 1453 /*134*/ NULL, 1454 NULL, 1455 NULL, 1456 NULL, 1457 NULL, 1458 /*139*/ NULL, 1459 NULL, 1460 NULL, 1461 NULL, 1462 NULL, 1463 /*144*/ NULL, 1464 NULL, 1465 NULL, 1466 NULL, 1467 NULL, 1468 NULL, 1469 /*150*/ NULL 1470 }; 1471 1472 static PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx) 1473 { 1474 Mat_Shell *shell = (Mat_Shell*)mat->data; 1475 1476 PetscFunctionBegin; 1477 if (ctx) { 1478 PetscContainer ctxcontainer; 1479 PetscCall(PetscContainerCreate(PetscObjectComm((PetscObject)mat),&ctxcontainer)); 1480 PetscCall(PetscContainerSetPointer(ctxcontainer,ctx)); 1481 PetscCall(PetscObjectCompose((PetscObject)mat,"MatShell ctx",(PetscObject)ctxcontainer)); 1482 shell->ctxcontainer = ctxcontainer; 1483 PetscCall(PetscContainerDestroy(&ctxcontainer)); 1484 } else { 1485 PetscCall(PetscObjectCompose((PetscObject)mat,"MatShell ctx",NULL)); 1486 shell->ctxcontainer = NULL; 1487 } 1488 1489 PetscFunctionReturn(0); 1490 } 1491 1492 PetscErrorCode MatShellSetContextDestroy_Shell(Mat mat,PetscErrorCode (*f)(void*)) 1493 { 1494 Mat_Shell *shell = (Mat_Shell*)mat->data; 1495 1496 PetscFunctionBegin; 1497 if (shell->ctxcontainer) PetscCall(PetscContainerSetUserDestroy(shell->ctxcontainer,f)); 1498 PetscFunctionReturn(0); 1499 } 1500 1501 static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype) 1502 { 1503 PetscFunctionBegin; 1504 PetscCall(PetscFree(mat->defaultvectype)); 1505 PetscCall(PetscStrallocpy(vtype,(char**)&mat->defaultvectype)); 1506 PetscFunctionReturn(0); 1507 } 1508 1509 PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A) 1510 { 1511 Mat_Shell *shell = (Mat_Shell*)A->data; 1512 1513 PetscFunctionBegin; 1514 shell->managescalingshifts = PETSC_FALSE; 1515 A->ops->diagonalset = NULL; 1516 A->ops->diagonalscale = NULL; 1517 A->ops->scale = NULL; 1518 A->ops->shift = NULL; 1519 A->ops->axpy = NULL; 1520 PetscFunctionReturn(0); 1521 } 1522 1523 static PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void)) 1524 { 1525 Mat_Shell *shell = (Mat_Shell*)mat->data; 1526 1527 PetscFunctionBegin; 1528 switch (op) { 1529 case MATOP_DESTROY: 1530 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1531 break; 1532 case MATOP_VIEW: 1533 if (!mat->ops->viewnative) { 1534 mat->ops->viewnative = mat->ops->view; 1535 } 1536 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1537 break; 1538 case MATOP_COPY: 1539 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1540 break; 1541 case MATOP_DIAGONAL_SET: 1542 case MATOP_DIAGONAL_SCALE: 1543 case MATOP_SHIFT: 1544 case MATOP_SCALE: 1545 case MATOP_AXPY: 1546 case MATOP_ZERO_ROWS: 1547 case MATOP_ZERO_ROWS_COLUMNS: 1548 PetscCheck(!shell->managescalingshifts,PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1549 (((void(**)(void))mat->ops)[op]) = f; 1550 break; 1551 case MATOP_GET_DIAGONAL: 1552 if (shell->managescalingshifts) { 1553 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1554 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1555 } else { 1556 shell->ops->getdiagonal = NULL; 1557 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1558 } 1559 break; 1560 case MATOP_MULT: 1561 if (shell->managescalingshifts) { 1562 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1563 mat->ops->mult = MatMult_Shell; 1564 } else { 1565 shell->ops->mult = NULL; 1566 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1567 } 1568 break; 1569 case MATOP_MULT_TRANSPOSE: 1570 if (shell->managescalingshifts) { 1571 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1572 mat->ops->multtranspose = MatMultTranspose_Shell; 1573 } else { 1574 shell->ops->multtranspose = NULL; 1575 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1576 } 1577 break; 1578 default: 1579 (((void(**)(void))mat->ops)[op]) = f; 1580 break; 1581 } 1582 PetscFunctionReturn(0); 1583 } 1584 1585 static PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void)) 1586 { 1587 Mat_Shell *shell = (Mat_Shell*)mat->data; 1588 1589 PetscFunctionBegin; 1590 switch (op) { 1591 case MATOP_DESTROY: 1592 *f = (void (*)(void))shell->ops->destroy; 1593 break; 1594 case MATOP_VIEW: 1595 *f = (void (*)(void))mat->ops->view; 1596 break; 1597 case MATOP_COPY: 1598 *f = (void (*)(void))shell->ops->copy; 1599 break; 1600 case MATOP_DIAGONAL_SET: 1601 case MATOP_DIAGONAL_SCALE: 1602 case MATOP_SHIFT: 1603 case MATOP_SCALE: 1604 case MATOP_AXPY: 1605 case MATOP_ZERO_ROWS: 1606 case MATOP_ZERO_ROWS_COLUMNS: 1607 *f = (((void (**)(void))mat->ops)[op]); 1608 break; 1609 case MATOP_GET_DIAGONAL: 1610 if (shell->ops->getdiagonal) 1611 *f = (void (*)(void))shell->ops->getdiagonal; 1612 else 1613 *f = (((void (**)(void))mat->ops)[op]); 1614 break; 1615 case MATOP_MULT: 1616 if (shell->ops->mult) 1617 *f = (void (*)(void))shell->ops->mult; 1618 else 1619 *f = (((void (**)(void))mat->ops)[op]); 1620 break; 1621 case MATOP_MULT_TRANSPOSE: 1622 if (shell->ops->multtranspose) 1623 *f = (void (*)(void))shell->ops->multtranspose; 1624 else 1625 *f = (((void (**)(void))mat->ops)[op]); 1626 break; 1627 default: 1628 *f = (((void (**)(void))mat->ops)[op]); 1629 } 1630 PetscFunctionReturn(0); 1631 } 1632 1633 /*MC 1634 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 1635 1636 Level: advanced 1637 1638 .seealso: `MatCreateShell()` 1639 M*/ 1640 1641 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 1642 { 1643 Mat_Shell *b; 1644 1645 PetscFunctionBegin; 1646 PetscCall(PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps))); 1647 1648 PetscCall(PetscNewLog(A,&b)); 1649 A->data = (void*)b; 1650 1651 b->ctxcontainer = NULL; 1652 b->vshift = 0.0; 1653 b->vscale = 1.0; 1654 b->managescalingshifts = PETSC_TRUE; 1655 A->assembled = PETSC_TRUE; 1656 A->preallocated = PETSC_FALSE; 1657 1658 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell)); 1659 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell)); 1660 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetContextDestroy_C",MatShellSetContextDestroy_Shell)); 1661 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell)); 1662 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell)); 1663 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell)); 1664 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell)); 1665 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatShellSetMatProductOperation_C",MatShellSetMatProductOperation_Shell)); 1666 PetscCall(PetscObjectChangeTypeName((PetscObject)A,MATSHELL)); 1667 PetscFunctionReturn(0); 1668 } 1669 1670 /*@C 1671 MatCreateShell - Creates a new matrix class for use with a user-defined 1672 private data storage format. 1673 1674 Collective 1675 1676 Input Parameters: 1677 + comm - MPI communicator 1678 . m - number of local rows (must be given) 1679 . n - number of local columns (must be given) 1680 . M - number of global rows (may be PETSC_DETERMINE) 1681 . N - number of global columns (may be PETSC_DETERMINE) 1682 - ctx - pointer to data needed by the shell matrix routines 1683 1684 Output Parameter: 1685 . A - the matrix 1686 1687 Level: advanced 1688 1689 Usage: 1690 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1691 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 1692 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1693 $ [ Use matrix for operations that have been set ] 1694 $ MatDestroy(mat); 1695 1696 Notes: 1697 The shell matrix type is intended to provide a simple class to use 1698 with KSP (such as, for use with matrix-free methods). You should not 1699 use the shell type if you plan to define a complete matrix class. 1700 1701 Fortran Notes: 1702 To use this from Fortran with a ctx you must write an interface definition for this 1703 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 1704 in as the ctx argument. 1705 1706 PETSc requires that matrices and vectors being used for certain 1707 operations are partitioned accordingly. For example, when 1708 creating a shell matrix, A, that supports parallel matrix-vector 1709 products using MatMult(A,x,y) the user should set the number 1710 of local matrix rows to be the number of local elements of the 1711 corresponding result vector, y. Note that this is information is 1712 required for use of the matrix interface routines, even though 1713 the shell matrix may not actually be physically partitioned. 1714 For example, 1715 1716 $ 1717 $ Vec x, y 1718 $ extern PetscErrorCode mult(Mat,Vec,Vec); 1719 $ Mat A 1720 $ 1721 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 1722 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 1723 $ VecGetLocalSize(y,&m); 1724 $ VecGetLocalSize(x,&n); 1725 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 1726 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 1727 $ MatMult(A,x,y); 1728 $ MatDestroy(&A); 1729 $ VecDestroy(&y); 1730 $ VecDestroy(&x); 1731 $ 1732 1733 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these 1734 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 1735 1736 For rectangular matrices do all the scalings and shifts make sense? 1737 1738 Developers Notes: 1739 Regarding shifting and scaling. The general form is 1740 1741 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 1742 1743 The order you apply the operations is important. For example if you have a dshift then 1744 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 1745 you get s*vscale*A + diag(shift) 1746 1747 A is the user provided function. 1748 1749 KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for 1750 for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger 1751 an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat); 1752 each time the MATSHELL matrix has changed. 1753 1754 Matrix product operations (i.e. MatMat, MatTransposeMat etc) can be specified using MatShellSetMatProductOperation() 1755 1756 Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided 1757 with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale(). 1758 1759 .seealso: `MatShellSetOperation()`, `MatHasOperation()`, `MatShellGetContext()`, `MatShellSetContext()`, `MATSHELL`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()` 1760 @*/ 1761 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 1762 { 1763 PetscFunctionBegin; 1764 PetscCall(MatCreate(comm,A)); 1765 PetscCall(MatSetSizes(*A,m,n,M,N)); 1766 PetscCall(MatSetType(*A,MATSHELL)); 1767 PetscCall(MatShellSetContext(*A,ctx)); 1768 PetscCall(MatSetUp(*A)); 1769 PetscFunctionReturn(0); 1770 } 1771 1772 /*@ 1773 MatShellSetContext - sets the context for a shell matrix 1774 1775 Logically Collective on Mat 1776 1777 Input Parameters: 1778 + mat - the shell matrix 1779 - ctx - the context 1780 1781 Level: advanced 1782 1783 Fortran Notes: 1784 To use this from Fortran you must write a Fortran interface definition for this 1785 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 1786 1787 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()` 1788 @*/ 1789 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 1790 { 1791 PetscFunctionBegin; 1792 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1793 PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx)); 1794 PetscFunctionReturn(0); 1795 } 1796 1797 /*@ 1798 MatShellSetContextDestroy - sets the destroy function for a shell matrix context 1799 1800 Logically Collective on Mat 1801 1802 Input Parameters: 1803 + mat - the shell matrix 1804 - f - the context destroy function 1805 1806 Note: 1807 If the `MatShell` is never duplicated, the behavior of this function is equivalent 1808 to `MatShellSetOperation(Mat,MATOP_DESTROY,f)`. However, `MatShellSetContextDestroy()` 1809 ensures proper reference counting for the user provided context data in the case that 1810 the `MatShell` is duplicated. 1811 1812 Level: advanced 1813 1814 .seealso: `MatCreateShell()`, `MatShellSetContext()` 1815 @*/ 1816 PetscErrorCode MatShellSetContextDestroy(Mat mat,PetscErrorCode (*f)(void*)) 1817 { 1818 PetscFunctionBegin; 1819 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1820 PetscTryMethod(mat,"MatShellSetContextDestroy_C",(Mat,PetscErrorCode (*)(void*)),(mat,f)); 1821 PetscFunctionReturn(0); 1822 } 1823 1824 /*@C 1825 MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs() 1826 1827 Logically collective 1828 1829 Input Parameters: 1830 + mat - the shell matrix 1831 - vtype - type to use for creating vectors 1832 1833 Notes: 1834 1835 Level: advanced 1836 1837 .seealso: `MatCreateVecs()` 1838 @*/ 1839 PetscErrorCode MatShellSetVecType(Mat mat,VecType vtype) 1840 { 1841 PetscFunctionBegin; 1842 PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype)); 1843 PetscFunctionReturn(0); 1844 } 1845 1846 /*@ 1847 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 1848 after MatCreateShell() 1849 1850 Logically Collective on Mat 1851 1852 Input Parameter: 1853 . mat - the shell matrix 1854 1855 Level: advanced 1856 1857 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatShellSetOperation()` 1858 @*/ 1859 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 1860 { 1861 PetscFunctionBegin; 1862 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 1863 PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A)); 1864 PetscFunctionReturn(0); 1865 } 1866 1867 /*@C 1868 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 1869 1870 Logically Collective on Mat 1871 1872 Input Parameters: 1873 + mat - the shell matrix 1874 . f - the function 1875 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1876 - ctx - an optional context for the function 1877 1878 Output Parameter: 1879 . flg - PETSC_TRUE if the multiply is likely correct 1880 1881 Options Database: 1882 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1883 1884 Level: advanced 1885 1886 Fortran Notes: 1887 Not supported from Fortran 1888 1889 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMultTranspose()` 1890 @*/ 1891 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1892 { 1893 PetscInt m,n; 1894 Mat mf,Dmf,Dmat,Ddiff; 1895 PetscReal Diffnorm,Dmfnorm; 1896 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1897 1898 PetscFunctionBegin; 1899 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1900 PetscCall(PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v)); 1901 PetscCall(MatGetLocalSize(mat,&m,&n)); 1902 PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf)); 1903 PetscCall(MatMFFDSetFunction(mf,f,ctx)); 1904 PetscCall(MatMFFDSetBase(mf,base,NULL)); 1905 1906 PetscCall(MatComputeOperator(mf,MATAIJ,&Dmf)); 1907 PetscCall(MatComputeOperator(mat,MATAIJ,&Dmat)); 1908 1909 PetscCall(MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff)); 1910 PetscCall(MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN)); 1911 PetscCall(MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm)); 1912 PetscCall(MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm)); 1913 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 1914 flag = PETSC_FALSE; 1915 if (v) { 1916 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))); 1917 PetscCall(MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view")); 1918 PetscCall(MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view")); 1919 PetscCall(MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view")); 1920 } 1921 } else if (v) { 1922 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n")); 1923 } 1924 if (flg) *flg = flag; 1925 PetscCall(MatDestroy(&Ddiff)); 1926 PetscCall(MatDestroy(&mf)); 1927 PetscCall(MatDestroy(&Dmf)); 1928 PetscCall(MatDestroy(&Dmat)); 1929 PetscFunctionReturn(0); 1930 } 1931 1932 /*@C 1933 MatShellTestMultTranspose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 1934 1935 Logically Collective on Mat 1936 1937 Input Parameters: 1938 + mat - the shell matrix 1939 . f - the function 1940 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 1941 - ctx - an optional context for the function 1942 1943 Output Parameter: 1944 . flg - PETSC_TRUE if the multiply is likely correct 1945 1946 Options Database: 1947 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 1948 1949 Level: advanced 1950 1951 Fortran Notes: 1952 Not supported from Fortran 1953 1954 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellTestMult()` 1955 @*/ 1956 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 1957 { 1958 Vec x,y,z; 1959 PetscInt m,n,M,N; 1960 Mat mf,Dmf,Dmat,Ddiff; 1961 PetscReal Diffnorm,Dmfnorm; 1962 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 1963 1964 PetscFunctionBegin; 1965 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1966 PetscCall(PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v)); 1967 PetscCall(MatCreateVecs(mat,&x,&y)); 1968 PetscCall(VecDuplicate(y,&z)); 1969 PetscCall(MatGetLocalSize(mat,&m,&n)); 1970 PetscCall(MatGetSize(mat,&M,&N)); 1971 PetscCall(MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf)); 1972 PetscCall(MatMFFDSetFunction(mf,f,ctx)); 1973 PetscCall(MatMFFDSetBase(mf,base,NULL)); 1974 PetscCall(MatComputeOperator(mf,MATAIJ,&Dmf)); 1975 PetscCall(MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf)); 1976 PetscCall(MatComputeOperatorTranspose(mat,MATAIJ,&Dmat)); 1977 1978 PetscCall(MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff)); 1979 PetscCall(MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN)); 1980 PetscCall(MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm)); 1981 PetscCall(MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm)); 1982 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 1983 flag = PETSC_FALSE; 1984 if (v) { 1985 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))); 1986 PetscCall(MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view")); 1987 PetscCall(MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view")); 1988 PetscCall(MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view")); 1989 } 1990 } else if (v) { 1991 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n")); 1992 } 1993 if (flg) *flg = flag; 1994 PetscCall(MatDestroy(&mf)); 1995 PetscCall(MatDestroy(&Dmat)); 1996 PetscCall(MatDestroy(&Ddiff)); 1997 PetscCall(MatDestroy(&Dmf)); 1998 PetscCall(VecDestroy(&x)); 1999 PetscCall(VecDestroy(&y)); 2000 PetscCall(VecDestroy(&z)); 2001 PetscFunctionReturn(0); 2002 } 2003 2004 /*@C 2005 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 2006 2007 Logically Collective on Mat 2008 2009 Input Parameters: 2010 + mat - the shell matrix 2011 . op - the name of the operation 2012 - g - the function that provides the operation. 2013 2014 Level: advanced 2015 2016 Usage: 2017 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 2018 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 2019 $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 2020 2021 Notes: 2022 See the file include/petscmat.h for a complete list of matrix 2023 operations, which all have the form MATOP_<OPERATION>, where 2024 <OPERATION> is the name (in all capital letters) of the 2025 user interface routine (e.g., MatMult() -> MATOP_MULT). 2026 2027 All user-provided functions (except for MATOP_DESTROY) should have the same calling 2028 sequence as the usual matrix interface routines, since they 2029 are intended to be accessed via the usual matrix interface 2030 routines, e.g., 2031 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 2032 2033 In particular each function MUST return an error code of 0 on success and 2034 nonzero on failure. 2035 2036 Within each user-defined routine, the user should call 2037 MatShellGetContext() to obtain the user-defined context that was 2038 set by MatCreateShell(). 2039 2040 Use MatSetOperation() to set an operation for any matrix type. For matrix product operations (i.e. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation() 2041 2042 Fortran Notes: 2043 For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 2044 generate a matrix. See src/mat/tests/ex120f.F 2045 2046 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellGetOperation()`, `MatShellSetContext()`, `MatSetOperation()`, `MatShellSetManageScalingShifts()`, `MatShellSetMatProductOperation()` 2047 @*/ 2048 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void)) 2049 { 2050 PetscFunctionBegin; 2051 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2052 PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g)); 2053 PetscFunctionReturn(0); 2054 } 2055 2056 /*@C 2057 MatShellGetOperation - Gets a matrix function for a shell matrix. 2058 2059 Not Collective 2060 2061 Input Parameters: 2062 + mat - the shell matrix 2063 - op - the name of the operation 2064 2065 Output Parameter: 2066 . g - the function that provides the operation. 2067 2068 Level: advanced 2069 2070 Notes: 2071 See the file include/petscmat.h for a complete list of matrix 2072 operations, which all have the form MATOP_<OPERATION>, where 2073 <OPERATION> is the name (in all capital letters) of the 2074 user interface routine (e.g., MatMult() -> MATOP_MULT). 2075 2076 All user-provided functions have the same calling 2077 sequence as the usual matrix interface routines, since they 2078 are intended to be accessed via the usual matrix interface 2079 routines, e.g., 2080 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 2081 2082 Within each user-defined routine, the user should call 2083 MatShellGetContext() to obtain the user-defined context that was 2084 set by MatCreateShell(). 2085 2086 .seealso: `MatCreateShell()`, `MatShellGetContext()`, `MatShellSetOperation()`, `MatShellSetContext()` 2087 @*/ 2088 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void)) 2089 { 2090 PetscFunctionBegin; 2091 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2092 PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g)); 2093 PetscFunctionReturn(0); 2094 } 2095 2096 /*@ 2097 MatIsShell - Inquires if a matrix is derived from MATSHELL 2098 2099 Input Parameter: 2100 . mat - the matrix 2101 2102 Output Parameter: 2103 . flg - the boolean value 2104 2105 Level: developer 2106 2107 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) 2108 2109 .seealso: `MatCreateShell()` 2110 @*/ 2111 PetscErrorCode MatIsShell(Mat mat, PetscBool *flg) 2112 { 2113 PetscFunctionBegin; 2114 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 2115 PetscValidBoolPointer(flg,2); 2116 *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell); 2117 PetscFunctionReturn(0); 2118 } 2119