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 10 struct _MatShellOps { 11 /* 3 */ PetscErrorCode (*mult)(Mat,Vec,Vec); 12 /* 5 */ PetscErrorCode (*multtranspose)(Mat,Vec,Vec); 13 /* 17 */ PetscErrorCode (*getdiagonal)(Mat,Vec); 14 /* 43 */ PetscErrorCode (*copy)(Mat,Mat,MatStructure); 15 /* 60 */ PetscErrorCode (*destroy)(Mat); 16 }; 17 18 typedef struct { 19 struct _MatShellOps ops[1]; 20 21 PetscScalar vscale,vshift; 22 Vec dshift; 23 Vec left,right; 24 Vec left_work,right_work; 25 Vec left_add_work,right_add_work; 26 Mat axpy; 27 PetscScalar axpy_vscale; 28 PetscBool managescalingshifts; /* The user will manage the scaling and shifts for the MATSHELL, not the default */ 29 void *ctx; 30 } Mat_Shell; 31 32 /* 33 xx = diag(left)*x 34 */ 35 static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx) 36 { 37 Mat_Shell *shell = (Mat_Shell*)A->data; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 *xx = NULL; 42 if (!shell->left) { 43 *xx = x; 44 } else { 45 if (!shell->left_work) {ierr = VecDuplicate(shell->left,&shell->left_work);CHKERRQ(ierr);} 46 ierr = VecPointwiseMult(shell->left_work,x,shell->left);CHKERRQ(ierr); 47 *xx = shell->left_work; 48 } 49 PetscFunctionReturn(0); 50 } 51 52 /* 53 xx = diag(right)*x 54 */ 55 static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx) 56 { 57 Mat_Shell *shell = (Mat_Shell*)A->data; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 *xx = NULL; 62 if (!shell->right) { 63 *xx = x; 64 } else { 65 if (!shell->right_work) {ierr = VecDuplicate(shell->right,&shell->right_work);CHKERRQ(ierr);} 66 ierr = VecPointwiseMult(shell->right_work,x,shell->right);CHKERRQ(ierr); 67 *xx = shell->right_work; 68 } 69 PetscFunctionReturn(0); 70 } 71 72 /* 73 x = diag(left)*x 74 */ 75 static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x) 76 { 77 Mat_Shell *shell = (Mat_Shell*)A->data; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 if (shell->left) {ierr = VecPointwiseMult(x,x,shell->left);CHKERRQ(ierr);} 82 PetscFunctionReturn(0); 83 } 84 85 /* 86 x = diag(right)*x 87 */ 88 static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x) 89 { 90 Mat_Shell *shell = (Mat_Shell*)A->data; 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 if (shell->right) {ierr = VecPointwiseMult(x,x,shell->right);CHKERRQ(ierr);} 95 PetscFunctionReturn(0); 96 } 97 98 /* 99 Y = vscale*Y + diag(dshift)*X + vshift*X 100 101 On input Y already contains A*x 102 */ 103 static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y) 104 { 105 Mat_Shell *shell = (Mat_Shell*)A->data; 106 PetscErrorCode ierr; 107 108 PetscFunctionBegin; 109 if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */ 110 PetscInt i,m; 111 const PetscScalar *x,*d; 112 PetscScalar *y; 113 ierr = VecGetLocalSize(X,&m);CHKERRQ(ierr); 114 ierr = VecGetArrayRead(shell->dshift,&d);CHKERRQ(ierr); 115 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 116 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 117 for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i]; 118 ierr = VecRestoreArrayRead(shell->dshift,&d);CHKERRQ(ierr); 119 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 120 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 121 } else { 122 ierr = VecScale(Y,shell->vscale);CHKERRQ(ierr); 123 } 124 if (shell->vshift != 0.0) {ierr = VecAXPY(Y,shell->vshift,X);CHKERRQ(ierr);} /* if test is for non-square matrices */ 125 PetscFunctionReturn(0); 126 } 127 128 /*@ 129 MatShellGetContext - Returns the user-provided context associated with a shell matrix. 130 131 Not Collective 132 133 Input Parameter: 134 . mat - the matrix, should have been created with MatCreateShell() 135 136 Output Parameter: 137 . ctx - the user provided context 138 139 Level: advanced 140 141 Fortran Notes: 142 To use this from Fortran you must write a Fortran interface definition for this 143 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 144 145 .keywords: matrix, shell, get, context 146 147 .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext() 148 @*/ 149 PetscErrorCode MatShellGetContext(Mat mat,void *ctx) 150 { 151 PetscErrorCode ierr; 152 PetscBool flg; 153 154 PetscFunctionBegin; 155 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 156 PetscValidPointer(ctx,2); 157 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 158 if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx; 159 else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix"); 160 PetscFunctionReturn(0); 161 } 162 163 PetscErrorCode MatDestroy_Shell(Mat mat) 164 { 165 PetscErrorCode ierr; 166 Mat_Shell *shell = (Mat_Shell*)mat->data; 167 168 PetscFunctionBegin; 169 if (shell->ops->destroy) { 170 ierr = (*shell->ops->destroy)(mat);CHKERRQ(ierr); 171 } 172 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 173 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 174 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 175 ierr = VecDestroy(&shell->left_work);CHKERRQ(ierr); 176 ierr = VecDestroy(&shell->right_work);CHKERRQ(ierr); 177 ierr = VecDestroy(&shell->left_add_work);CHKERRQ(ierr); 178 ierr = VecDestroy(&shell->right_add_work);CHKERRQ(ierr); 179 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 180 ierr = PetscFree(mat->data);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str) 185 { 186 Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data; 187 PetscErrorCode ierr; 188 PetscBool matflg; 189 190 PetscFunctionBegin; 191 ierr = PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);CHKERRQ(ierr); 192 if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL"); 193 194 ierr = PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 195 ierr = PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));CHKERRQ(ierr); 196 197 if (shellA->ops->copy) { 198 ierr = (*shellA->ops->copy)(A,B,str);CHKERRQ(ierr); 199 } 200 shellB->vscale = shellA->vscale; 201 shellB->vshift = shellA->vshift; 202 if (shellA->dshift) { 203 if (!shellB->dshift) { 204 ierr = VecDuplicate(shellA->dshift,&shellB->dshift);CHKERRQ(ierr); 205 } 206 ierr = VecCopy(shellA->dshift,shellB->dshift);CHKERRQ(ierr); 207 } else { 208 ierr = VecDestroy(&shellB->dshift);CHKERRQ(ierr); 209 } 210 if (shellA->left) { 211 if (!shellB->left) { 212 ierr = VecDuplicate(shellA->left,&shellB->left);CHKERRQ(ierr); 213 } 214 ierr = VecCopy(shellA->left,shellB->left);CHKERRQ(ierr); 215 } else { 216 ierr = VecDestroy(&shellB->left);CHKERRQ(ierr); 217 } 218 if (shellA->right) { 219 if (!shellB->right) { 220 ierr = VecDuplicate(shellA->right,&shellB->right);CHKERRQ(ierr); 221 } 222 ierr = VecCopy(shellA->right,shellB->right);CHKERRQ(ierr); 223 } else { 224 ierr = VecDestroy(&shellB->right);CHKERRQ(ierr); 225 } 226 ierr = MatDestroy(&shellB->axpy);CHKERRQ(ierr); 227 if (shellA->axpy) { 228 ierr = PetscObjectReference((PetscObject)shellA->axpy);CHKERRQ(ierr); 229 shellB->axpy = shellA->axpy; 230 shellB->axpy_vscale = shellA->axpy_vscale; 231 } 232 PetscFunctionReturn(0); 233 } 234 235 PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M) 236 { 237 PetscErrorCode ierr; 238 void *ctx; 239 240 PetscFunctionBegin; 241 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 242 ierr = MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);CHKERRQ(ierr); 243 ierr = MatCopy(mat,*M,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 244 PetscFunctionReturn(0); 245 } 246 247 PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y) 248 { 249 Mat_Shell *shell = (Mat_Shell*)A->data; 250 PetscErrorCode ierr; 251 Vec xx; 252 PetscObjectState instate,outstate; 253 254 PetscFunctionBegin; 255 ierr = MatShellPreScaleRight(A,x,&xx);CHKERRQ(ierr); 256 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 257 if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL"); 258 ierr = (*shell->ops->mult)(A,xx,y);CHKERRQ(ierr); 259 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 260 if (instate == outstate) { 261 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 262 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 263 } 264 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 265 ierr = MatShellPostScaleLeft(A,y);CHKERRQ(ierr); 266 267 if (shell->axpy) { 268 if (!shell->left_work) {ierr = MatCreateVecs(A,&shell->left_work,NULL);CHKERRQ(ierr);} 269 ierr = MatMult(shell->axpy,x,shell->left_work);CHKERRQ(ierr); 270 ierr = VecAXPY(y,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr); 271 } 272 PetscFunctionReturn(0); 273 } 274 275 PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z) 276 { 277 Mat_Shell *shell = (Mat_Shell*)A->data; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 if (y == z) { 282 if (!shell->right_add_work) {ierr = VecDuplicate(z,&shell->right_add_work);CHKERRQ(ierr);} 283 ierr = MatMult(A,x,shell->right_add_work);CHKERRQ(ierr); 284 ierr = VecAXPY(z,1.0,shell->right_add_work);CHKERRQ(ierr); 285 } else { 286 ierr = MatMult(A,x,z);CHKERRQ(ierr); 287 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 288 } 289 PetscFunctionReturn(0); 290 } 291 292 PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y) 293 { 294 Mat_Shell *shell = (Mat_Shell*)A->data; 295 PetscErrorCode ierr; 296 Vec xx; 297 PetscObjectState instate,outstate; 298 299 PetscFunctionBegin; 300 ierr = MatShellPreScaleLeft(A,x,&xx);CHKERRQ(ierr); 301 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 302 if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL"); 303 ierr = (*shell->ops->multtranspose)(A,xx,y);CHKERRQ(ierr); 304 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 305 if (instate == outstate) { 306 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 307 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 308 } 309 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 310 ierr = MatShellPostScaleRight(A,y);CHKERRQ(ierr); 311 PetscFunctionReturn(0); 312 } 313 314 PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z) 315 { 316 Mat_Shell *shell = (Mat_Shell*)A->data; 317 PetscErrorCode ierr; 318 319 PetscFunctionBegin; 320 if (y == z) { 321 if (!shell->left_add_work) {ierr = VecDuplicate(z,&shell->left_add_work);CHKERRQ(ierr);} 322 ierr = MatMultTranspose(A,x,shell->left_add_work);CHKERRQ(ierr); 323 ierr = VecWAXPY(z,1.0,shell->left_add_work,y);CHKERRQ(ierr); 324 } else { 325 ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr); 326 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 327 } 328 PetscFunctionReturn(0); 329 } 330 331 /* 332 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 333 */ 334 PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v) 335 { 336 Mat_Shell *shell = (Mat_Shell*)A->data; 337 PetscErrorCode ierr; 338 339 PetscFunctionBegin; 340 if (shell->ops->getdiagonal) { 341 ierr = (*shell->ops->getdiagonal)(A,v);CHKERRQ(ierr); 342 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)"); 343 ierr = VecScale(v,shell->vscale);CHKERRQ(ierr); 344 if (shell->dshift) { 345 ierr = VecAXPY(v,1.0,shell->dshift);CHKERRQ(ierr); 346 } 347 ierr = VecShift(v,shell->vshift);CHKERRQ(ierr); 348 if (shell->left) {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);} 349 if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);} 350 if (shell->axpy) { 351 if (!shell->left_work) {ierr = VecDuplicate(v,&shell->left_work);CHKERRQ(ierr);} 352 ierr = MatGetDiagonal(shell->axpy,shell->left_work);CHKERRQ(ierr); 353 ierr = VecAXPY(v,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr); 354 } 355 PetscFunctionReturn(0); 356 } 357 358 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 359 { 360 Mat_Shell *shell = (Mat_Shell*)Y->data; 361 PetscErrorCode ierr; 362 363 PetscFunctionBegin; 364 if (shell->left || shell->right) { 365 if (!shell->dshift) { 366 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 367 ierr = VecSet(shell->dshift,a);CHKERRQ(ierr); 368 } else { 369 if (shell->left) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 370 if (shell->right) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 371 ierr = VecShift(shell->dshift,a);CHKERRQ(ierr); 372 } 373 if (shell->left) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 374 if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 375 } else shell->vshift += a; 376 PetscFunctionReturn(0); 377 } 378 379 PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins) 380 { 381 Mat_Shell *shell = (Mat_Shell*)A->data; 382 PetscErrorCode ierr; 383 384 PetscFunctionBegin; 385 if (ins == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported with INSERT_VALUES"); 386 if (!shell->dshift) {ierr = VecDuplicate(D,&shell->dshift);CHKERRQ(ierr);} 387 if (shell->left || shell->right) { 388 if (!shell->right_work) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);CHKERRQ(ierr);} 389 if (shell->left && shell->right) { 390 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 391 ierr = VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);CHKERRQ(ierr); 392 } else if (shell->left) { 393 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 394 } else { 395 ierr = VecPointwiseDivide(shell->right_work,D,shell->right);CHKERRQ(ierr); 396 } 397 if (!shell->dshift) { 398 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 399 ierr = VecCopy(shell->dshift,shell->right_work);CHKERRQ(ierr); 400 } else { 401 ierr = VecAXPY(shell->dshift,1.0,shell->right_work);CHKERRQ(ierr); 402 } 403 } else { 404 ierr = VecAXPY(shell->dshift,1.0,D);CHKERRQ(ierr); 405 } 406 PetscFunctionReturn(0); 407 } 408 409 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 410 { 411 Mat_Shell *shell = (Mat_Shell*)Y->data; 412 PetscErrorCode ierr; 413 414 PetscFunctionBegin; 415 shell->vscale *= a; 416 shell->vshift *= a; 417 if (shell->dshift) { 418 ierr = VecScale(shell->dshift,a);CHKERRQ(ierr); 419 } 420 shell->axpy_vscale *= a; 421 PetscFunctionReturn(0); 422 } 423 424 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right) 425 { 426 Mat_Shell *shell = (Mat_Shell*)Y->data; 427 PetscErrorCode ierr; 428 429 PetscFunctionBegin; 430 if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation"); 431 if (left) { 432 if (!shell->left) { 433 ierr = VecDuplicate(left,&shell->left);CHKERRQ(ierr); 434 ierr = VecCopy(left,shell->left);CHKERRQ(ierr); 435 } else { 436 ierr = VecPointwiseMult(shell->left,shell->left,left);CHKERRQ(ierr); 437 } 438 } 439 if (right) { 440 if (!shell->right) { 441 ierr = VecDuplicate(right,&shell->right);CHKERRQ(ierr); 442 ierr = VecCopy(right,shell->right);CHKERRQ(ierr); 443 } else { 444 ierr = VecPointwiseMult(shell->right,shell->right,right);CHKERRQ(ierr); 445 } 446 } 447 PetscFunctionReturn(0); 448 } 449 450 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 451 { 452 Mat_Shell *shell = (Mat_Shell*)Y->data; 453 PetscErrorCode ierr; 454 455 PetscFunctionBegin; 456 if (t == MAT_FINAL_ASSEMBLY) { 457 shell->vshift = 0.0; 458 shell->vscale = 1.0; 459 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 460 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 461 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 462 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 463 } 464 PetscFunctionReturn(0); 465 } 466 467 PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*); 468 469 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d) 470 { 471 PetscFunctionBegin; 472 *missing = PETSC_FALSE; 473 PetscFunctionReturn(0); 474 } 475 476 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str) 477 { 478 Mat_Shell *shell = (Mat_Shell*)Y->data; 479 PetscErrorCode ierr; 480 481 PetscFunctionBegin; 482 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 483 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 484 shell->axpy = X; 485 shell->axpy_vscale = a; 486 PetscFunctionReturn(0); 487 } 488 489 static struct _MatOps MatOps_Values = {0, 490 0, 491 0, 492 0, 493 /* 4*/ MatMultAdd_Shell, 494 0, 495 MatMultTransposeAdd_Shell, 496 0, 497 0, 498 0, 499 /*10*/ 0, 500 0, 501 0, 502 0, 503 0, 504 /*15*/ 0, 505 0, 506 0, 507 MatDiagonalScale_Shell, 508 0, 509 /*20*/ 0, 510 MatAssemblyEnd_Shell, 511 0, 512 0, 513 /*24*/ 0, 514 0, 515 0, 516 0, 517 0, 518 /*29*/ 0, 519 0, 520 0, 521 0, 522 0, 523 /*34*/ MatDuplicate_Shell, 524 0, 525 0, 526 0, 527 0, 528 /*39*/ MatAXPY_Shell, 529 0, 530 0, 531 0, 532 MatCopy_Shell, 533 /*44*/ 0, 534 MatScale_Shell, 535 MatShift_Shell, 536 MatDiagonalSet_Shell, 537 0, 538 /*49*/ 0, 539 0, 540 0, 541 0, 542 0, 543 /*54*/ 0, 544 0, 545 0, 546 0, 547 0, 548 /*59*/ 0, 549 MatDestroy_Shell, 550 0, 551 0, 552 0, 553 /*64*/ 0, 554 0, 555 0, 556 0, 557 0, 558 /*69*/ 0, 559 0, 560 MatConvert_Shell, 561 0, 562 0, 563 /*74*/ 0, 564 0, 565 0, 566 0, 567 0, 568 /*79*/ 0, 569 0, 570 0, 571 0, 572 0, 573 /*84*/ 0, 574 0, 575 0, 576 0, 577 0, 578 /*89*/ 0, 579 0, 580 0, 581 0, 582 0, 583 /*94*/ 0, 584 0, 585 0, 586 0, 587 0, 588 /*99*/ 0, 589 0, 590 0, 591 0, 592 0, 593 /*104*/ 0, 594 0, 595 0, 596 0, 597 0, 598 /*109*/ 0, 599 0, 600 0, 601 0, 602 MatMissingDiagonal_Shell, 603 /*114*/ 0, 604 0, 605 0, 606 0, 607 0, 608 /*119*/ 0, 609 0, 610 0, 611 0, 612 0, 613 /*124*/ 0, 614 0, 615 0, 616 0, 617 0, 618 /*129*/ 0, 619 0, 620 0, 621 0, 622 0, 623 /*134*/ 0, 624 0, 625 0, 626 0, 627 0, 628 /*139*/ 0, 629 0, 630 0 631 }; 632 633 /*MC 634 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 635 636 Level: advanced 637 638 .seealso: MatCreateShell() 639 M*/ 640 641 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 642 { 643 Mat_Shell *b; 644 PetscErrorCode ierr; 645 646 PetscFunctionBegin; 647 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 648 649 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 650 A->data = (void*)b; 651 652 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 653 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 654 655 b->ctx = 0; 656 b->vshift = 0.0; 657 b->vscale = 1.0; 658 b->managescalingshifts = PETSC_TRUE; 659 A->assembled = PETSC_TRUE; 660 A->preallocated = PETSC_FALSE; 661 662 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 663 PetscFunctionReturn(0); 664 } 665 666 /*@C 667 MatCreateShell - Creates a new matrix class for use with a user-defined 668 private data storage format. 669 670 Collective on MPI_Comm 671 672 Input Parameters: 673 + comm - MPI communicator 674 . m - number of local rows (must be given) 675 . n - number of local columns (must be given) 676 . M - number of global rows (may be PETSC_DETERMINE) 677 . N - number of global columns (may be PETSC_DETERMINE) 678 - ctx - pointer to data needed by the shell matrix routines 679 680 Output Parameter: 681 . A - the matrix 682 683 Level: advanced 684 685 Usage: 686 $ extern int mult(Mat,Vec,Vec); 687 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 688 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 689 $ [ Use matrix for operations that have been set ] 690 $ MatDestroy(mat); 691 692 Notes: 693 The shell matrix type is intended to provide a simple class to use 694 with KSP (such as, for use with matrix-free methods). You should not 695 use the shell type if you plan to define a complete matrix class. 696 697 Fortran Notes: 698 To use this from Fortran with a ctx you must write an interface definition for this 699 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 700 in as the ctx argument. 701 702 PETSc requires that matrices and vectors being used for certain 703 operations are partitioned accordingly. For example, when 704 creating a shell matrix, A, that supports parallel matrix-vector 705 products using MatMult(A,x,y) the user should set the number 706 of local matrix rows to be the number of local elements of the 707 corresponding result vector, y. Note that this is information is 708 required for use of the matrix interface routines, even though 709 the shell matrix may not actually be physically partitioned. 710 For example, 711 712 $ 713 $ Vec x, y 714 $ extern int mult(Mat,Vec,Vec); 715 $ Mat A 716 $ 717 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 718 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 719 $ VecGetLocalSize(y,&m); 720 $ VecGetLocalSize(x,&n); 721 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 722 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 723 $ MatMult(A,x,y); 724 $ MatDestroy(A); 725 $ VecDestroy(y); VecDestroy(x); 726 $ 727 728 729 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), and MatScale() internally so these 730 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 731 732 733 For rectangular matrices do all the scalings and shifts make sense? 734 735 Developers Notes: 736 Regarding shifting and scaling. The general form is 737 738 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 739 740 The order you apply the operations is important. For example if you have a dshift then 741 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 742 you get s*vscale*A + diag(shift) 743 744 A is the user provided function. 745 746 .keywords: matrix, shell, create 747 748 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts() 749 @*/ 750 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 751 { 752 PetscErrorCode ierr; 753 754 PetscFunctionBegin; 755 ierr = MatCreate(comm,A);CHKERRQ(ierr); 756 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 757 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 758 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 759 ierr = MatSetUp(*A);CHKERRQ(ierr); 760 PetscFunctionReturn(0); 761 } 762 763 /*@ 764 MatShellSetContext - sets the context for a shell matrix 765 766 Logically Collective on Mat 767 768 Input Parameters: 769 + mat - the shell matrix 770 - ctx - the context 771 772 Level: advanced 773 774 Fortran Notes: 775 To use this from Fortran you must write a Fortran interface definition for this 776 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 777 778 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 779 @*/ 780 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 781 { 782 Mat_Shell *shell = (Mat_Shell*)mat->data; 783 PetscErrorCode ierr; 784 PetscBool flg; 785 786 PetscFunctionBegin; 787 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 788 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 789 if (flg) { 790 shell->ctx = ctx; 791 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix"); 792 PetscFunctionReturn(0); 793 } 794 795 /*@ 796 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 797 after MatCreateShell() 798 799 Logically Collective on Mat 800 801 Input Parameter: 802 . mat - the shell matrix 803 804 Level: advanced 805 806 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation() 807 @*/ 808 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 809 { 810 PetscErrorCode ierr; 811 Mat_Shell *shell; 812 PetscBool flg; 813 814 PetscFunctionBegin; 815 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 816 ierr = PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);CHKERRQ(ierr); 817 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 818 shell = (Mat_Shell*)A->data; 819 shell->managescalingshifts = PETSC_FALSE; 820 A->ops->diagonalset = NULL; 821 A->ops->diagonalscale = NULL; 822 A->ops->scale = NULL; 823 A->ops->shift = NULL; 824 A->ops->axpy = NULL; 825 PetscFunctionReturn(0); 826 } 827 828 /*@C 829 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 830 831 Logically Collective on Mat 832 833 Input Parameters: 834 + mat - the shell matrix 835 . f - the function 836 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 837 - ctx - an optional context for the function 838 839 Output Parameter: 840 . flg - PETSC_TRUE if the multiply is likely correct 841 842 Options Database: 843 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 844 845 Level: advanced 846 847 Fortran Notes: 848 Not supported from Fortran 849 850 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose() 851 @*/ 852 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 853 { 854 PetscErrorCode ierr; 855 PetscInt m,n; 856 Mat mf,Dmf,Dmat,Ddiff; 857 PetscReal Diffnorm,Dmfnorm; 858 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 859 860 PetscFunctionBegin; 861 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 862 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr); 863 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 864 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr); 865 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 866 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 867 868 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 869 ierr = MatComputeExplicitOperator(mat,&Dmat);CHKERRQ(ierr); 870 871 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 872 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 873 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 874 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 875 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 876 flag = PETSC_FALSE; 877 if (v) { 878 ierr = 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));CHKERRQ(ierr); 879 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 880 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 881 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 882 } 883 } else if (v) { 884 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 885 } 886 if (flg) *flg = flag; 887 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 888 ierr = MatDestroy(&mf);CHKERRQ(ierr); 889 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 890 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 891 PetscFunctionReturn(0); 892 } 893 894 /*@C 895 MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 896 897 Logically Collective on Mat 898 899 Input Parameters: 900 + mat - the shell matrix 901 . f - the function 902 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 903 - ctx - an optional context for the function 904 905 Output Parameter: 906 . flg - PETSC_TRUE if the multiply is likely correct 907 908 Options Database: 909 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 910 911 Level: advanced 912 913 Fortran Notes: 914 Not supported from Fortran 915 916 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult() 917 @*/ 918 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 919 { 920 PetscErrorCode ierr; 921 Vec x,y,z; 922 PetscInt m,n,M,N; 923 Mat mf,Dmf,Dmat,Ddiff; 924 PetscReal Diffnorm,Dmfnorm; 925 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 926 927 PetscFunctionBegin; 928 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 929 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr); 930 ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr); 931 ierr = VecDuplicate(y,&z);CHKERRQ(ierr); 932 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 933 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 934 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr); 935 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 936 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 937 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 938 ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr); 939 ierr = MatComputeExplicitOperatorTranspose(mat,&Dmat);CHKERRQ(ierr); 940 941 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 942 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 943 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 944 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 945 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 946 flag = PETSC_FALSE; 947 if (v) { 948 ierr = 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));CHKERRQ(ierr); 949 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 950 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 951 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 952 } 953 } else if (v) { 954 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 955 } 956 if (flg) *flg = flag; 957 ierr = MatDestroy(&mf);CHKERRQ(ierr); 958 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 959 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 960 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 961 ierr = VecDestroy(&x);CHKERRQ(ierr); 962 ierr = VecDestroy(&y);CHKERRQ(ierr); 963 ierr = VecDestroy(&z);CHKERRQ(ierr); 964 PetscFunctionReturn(0); 965 } 966 967 /*@C 968 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 969 970 Logically Collective on Mat 971 972 Input Parameters: 973 + mat - the shell matrix 974 . op - the name of the operation 975 - f - the function that provides the operation. 976 977 Level: advanced 978 979 Usage: 980 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 981 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 982 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 983 984 Notes: 985 See the file include/petscmat.h for a complete list of matrix 986 operations, which all have the form MATOP_<OPERATION>, where 987 <OPERATION> is the name (in all capital letters) of the 988 user interface routine (e.g., MatMult() -> MATOP_MULT). 989 990 All user-provided functions (except for MATOP_DESTROY) should have the same calling 991 sequence as the usual matrix interface routines, since they 992 are intended to be accessed via the usual matrix interface 993 routines, e.g., 994 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 995 996 In particular each function MUST return an error code of 0 on success and 997 nonzero on failure. 998 999 Within each user-defined routine, the user should call 1000 MatShellGetContext() to obtain the user-defined context that was 1001 set by MatCreateShell(). 1002 1003 Fortran Notes: 1004 For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 1005 generate a matrix. See src/mat/examples/tests/ex120f.F 1006 1007 Use MatSetOperation() to set an operation for any matrix type 1008 1009 .keywords: matrix, shell, set, operation 1010 1011 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts() 1012 @*/ 1013 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 1014 { 1015 PetscBool flg; 1016 Mat_Shell *shell; 1017 PetscErrorCode ierr; 1018 1019 PetscFunctionBegin; 1020 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1021 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1022 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1023 shell = (Mat_Shell*)mat->data; 1024 1025 switch (op) { 1026 case MATOP_DESTROY: 1027 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1028 break; 1029 case MATOP_VIEW: 1030 if (!mat->ops->viewnative) { 1031 mat->ops->viewnative = mat->ops->view; 1032 } 1033 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1034 break; 1035 case MATOP_COPY: 1036 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1037 break; 1038 case MATOP_DIAGONAL_SET: 1039 case MATOP_DIAGONAL_SCALE: 1040 case MATOP_SHIFT: 1041 case MATOP_SCALE: 1042 case MATOP_AXPY: 1043 if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1044 (((void(**)(void))mat->ops)[op]) = f; 1045 break; 1046 case MATOP_GET_DIAGONAL: 1047 if (shell->managescalingshifts) { 1048 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1049 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1050 } else { 1051 shell->ops->getdiagonal = NULL; 1052 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1053 } 1054 break; 1055 case MATOP_MULT: 1056 if (shell->managescalingshifts) { 1057 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1058 mat->ops->mult = MatMult_Shell; 1059 } else { 1060 shell->ops->mult = NULL; 1061 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1062 } 1063 break; 1064 case MATOP_MULT_TRANSPOSE: 1065 if (shell->managescalingshifts) { 1066 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1067 mat->ops->multtranspose = MatMultTranspose_Shell; 1068 } else { 1069 shell->ops->multtranspose = NULL; 1070 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1071 } 1072 break; 1073 default: 1074 (((void(**)(void))mat->ops)[op]) = f; 1075 break; 1076 } 1077 PetscFunctionReturn(0); 1078 } 1079 1080 /*@C 1081 MatShellGetOperation - Gets a matrix function for a shell matrix. 1082 1083 Not Collective 1084 1085 Input Parameters: 1086 + mat - the shell matrix 1087 - op - the name of the operation 1088 1089 Output Parameter: 1090 . f - the function that provides the operation. 1091 1092 Level: advanced 1093 1094 Notes: 1095 See the file include/petscmat.h for a complete list of matrix 1096 operations, which all have the form MATOP_<OPERATION>, where 1097 <OPERATION> is the name (in all capital letters) of the 1098 user interface routine (e.g., MatMult() -> MATOP_MULT). 1099 1100 All user-provided functions have the same calling 1101 sequence as the usual matrix interface routines, since they 1102 are intended to be accessed via the usual matrix interface 1103 routines, e.g., 1104 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1105 1106 Within each user-defined routine, the user should call 1107 MatShellGetContext() to obtain the user-defined context that was 1108 set by MatCreateShell(). 1109 1110 .keywords: matrix, shell, set, operation 1111 1112 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 1113 @*/ 1114 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 1115 { 1116 PetscBool flg; 1117 Mat_Shell *shell; 1118 PetscErrorCode ierr; 1119 1120 PetscFunctionBegin; 1121 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1122 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1123 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1124 shell = (Mat_Shell*)mat->data; 1125 1126 switch (op) { 1127 case MATOP_DESTROY: 1128 *f = (void (*)(void))shell->ops->destroy; 1129 break; 1130 case MATOP_VIEW: 1131 *f = (void (*)(void))mat->ops->view; 1132 break; 1133 case MATOP_COPY: 1134 *f = (void (*)(void))shell->ops->copy; 1135 break; 1136 case MATOP_DIAGONAL_SET: 1137 case MATOP_DIAGONAL_SCALE: 1138 case MATOP_SHIFT: 1139 case MATOP_SCALE: 1140 case MATOP_AXPY: 1141 *f = (((void (**)(void))mat->ops)[op]); 1142 break; 1143 case MATOP_GET_DIAGONAL: 1144 if (shell->ops->getdiagonal) 1145 *f = (void (*)(void))shell->ops->getdiagonal; 1146 else 1147 *f = (((void (**)(void))mat->ops)[op]); 1148 break; 1149 case MATOP_MULT: 1150 if (shell->ops->mult) 1151 *f = (void (*)(void))shell->ops->mult; 1152 else 1153 *f = (((void (**)(void))mat->ops)[op]); 1154 break; 1155 case MATOP_MULT_TRANSPOSE: 1156 if (shell->ops->multtranspose) 1157 *f = (void (*)(void))shell->ops->multtranspose; 1158 else 1159 *f = (((void (**)(void))mat->ops)[op]); 1160 break; 1161 default: 1162 *f = (((void (**)(void))mat->ops)[op]); 1163 } 1164 PetscFunctionReturn(0); 1165 } 1166