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