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 = VecAXPY(z,1.0,shell->left_add_work);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 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d) 468 { 469 PetscFunctionBegin; 470 *missing = PETSC_FALSE; 471 PetscFunctionReturn(0); 472 } 473 474 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str) 475 { 476 Mat_Shell *shell = (Mat_Shell*)Y->data; 477 PetscErrorCode ierr; 478 479 PetscFunctionBegin; 480 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 481 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 482 shell->axpy = X; 483 shell->axpy_vscale = a; 484 PetscFunctionReturn(0); 485 } 486 487 static struct _MatOps MatOps_Values = {0, 488 0, 489 0, 490 0, 491 /* 4*/ MatMultAdd_Shell, 492 0, 493 MatMultTransposeAdd_Shell, 494 0, 495 0, 496 0, 497 /*10*/ 0, 498 0, 499 0, 500 0, 501 0, 502 /*15*/ 0, 503 0, 504 0, 505 MatDiagonalScale_Shell, 506 0, 507 /*20*/ 0, 508 MatAssemblyEnd_Shell, 509 0, 510 0, 511 /*24*/ 0, 512 0, 513 0, 514 0, 515 0, 516 /*29*/ 0, 517 0, 518 0, 519 0, 520 0, 521 /*34*/ MatDuplicate_Shell, 522 0, 523 0, 524 0, 525 0, 526 /*39*/ MatAXPY_Shell, 527 0, 528 0, 529 0, 530 MatCopy_Shell, 531 /*44*/ 0, 532 MatScale_Shell, 533 MatShift_Shell, 534 MatDiagonalSet_Shell, 535 0, 536 /*49*/ 0, 537 0, 538 0, 539 0, 540 0, 541 /*54*/ 0, 542 0, 543 0, 544 0, 545 0, 546 /*59*/ 0, 547 MatDestroy_Shell, 548 0, 549 0, 550 0, 551 /*64*/ 0, 552 0, 553 0, 554 0, 555 0, 556 /*69*/ 0, 557 0, 558 MatConvert_Shell, 559 0, 560 0, 561 /*74*/ 0, 562 0, 563 0, 564 0, 565 0, 566 /*79*/ 0, 567 0, 568 0, 569 0, 570 0, 571 /*84*/ 0, 572 0, 573 0, 574 0, 575 0, 576 /*89*/ 0, 577 0, 578 0, 579 0, 580 0, 581 /*94*/ 0, 582 0, 583 0, 584 0, 585 0, 586 /*99*/ 0, 587 0, 588 0, 589 0, 590 0, 591 /*104*/ 0, 592 0, 593 0, 594 0, 595 0, 596 /*109*/ 0, 597 0, 598 0, 599 0, 600 MatMissingDiagonal_Shell, 601 /*114*/ 0, 602 0, 603 0, 604 0, 605 0, 606 /*119*/ 0, 607 0, 608 0, 609 0, 610 0, 611 /*124*/ 0, 612 0, 613 0, 614 0, 615 0, 616 /*129*/ 0, 617 0, 618 0, 619 0, 620 0, 621 /*134*/ 0, 622 0, 623 0, 624 0, 625 0, 626 /*139*/ 0, 627 0, 628 0 629 }; 630 631 /*MC 632 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 633 634 Level: advanced 635 636 .seealso: MatCreateShell() 637 M*/ 638 639 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 640 { 641 Mat_Shell *b; 642 PetscErrorCode ierr; 643 644 PetscFunctionBegin; 645 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 646 647 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 648 A->data = (void*)b; 649 650 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 651 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 652 653 b->ctx = 0; 654 b->vshift = 0.0; 655 b->vscale = 1.0; 656 b->managescalingshifts = PETSC_TRUE; 657 A->assembled = PETSC_TRUE; 658 A->preallocated = PETSC_FALSE; 659 660 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 661 PetscFunctionReturn(0); 662 } 663 664 /*@C 665 MatCreateShell - Creates a new matrix class for use with a user-defined 666 private data storage format. 667 668 Collective on MPI_Comm 669 670 Input Parameters: 671 + comm - MPI communicator 672 . m - number of local rows (must be given) 673 . n - number of local columns (must be given) 674 . M - number of global rows (may be PETSC_DETERMINE) 675 . N - number of global columns (may be PETSC_DETERMINE) 676 - ctx - pointer to data needed by the shell matrix routines 677 678 Output Parameter: 679 . A - the matrix 680 681 Level: advanced 682 683 Usage: 684 $ extern int mult(Mat,Vec,Vec); 685 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 686 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 687 $ [ Use matrix for operations that have been set ] 688 $ MatDestroy(mat); 689 690 Notes: 691 The shell matrix type is intended to provide a simple class to use 692 with KSP (such as, for use with matrix-free methods). You should not 693 use the shell type if you plan to define a complete matrix class. 694 695 Fortran Notes: 696 To use this from Fortran with a ctx you must write an interface definition for this 697 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 698 in as the ctx argument. 699 700 PETSc requires that matrices and vectors being used for certain 701 operations are partitioned accordingly. For example, when 702 creating a shell matrix, A, that supports parallel matrix-vector 703 products using MatMult(A,x,y) the user should set the number 704 of local matrix rows to be the number of local elements of the 705 corresponding result vector, y. Note that this is information is 706 required for use of the matrix interface routines, even though 707 the shell matrix may not actually be physically partitioned. 708 For example, 709 710 $ 711 $ Vec x, y 712 $ extern int mult(Mat,Vec,Vec); 713 $ Mat A 714 $ 715 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 716 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 717 $ VecGetLocalSize(y,&m); 718 $ VecGetLocalSize(x,&n); 719 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 720 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 721 $ MatMult(A,x,y); 722 $ MatDestroy(A); 723 $ VecDestroy(y); VecDestroy(x); 724 $ 725 726 727 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), and MatScale() internally so these 728 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 729 730 731 For rectangular matrices do all the scalings and shifts make sense? 732 733 Developers Notes: 734 Regarding shifting and scaling. The general form is 735 736 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 737 738 The order you apply the operations is important. For example if you have a dshift then 739 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 740 you get s*vscale*A + diag(shift) 741 742 A is the user provided function. 743 744 .keywords: matrix, shell, create 745 746 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts() 747 @*/ 748 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 749 { 750 PetscErrorCode ierr; 751 752 PetscFunctionBegin; 753 ierr = MatCreate(comm,A);CHKERRQ(ierr); 754 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 755 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 756 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 757 ierr = MatSetUp(*A);CHKERRQ(ierr); 758 PetscFunctionReturn(0); 759 } 760 761 /*@ 762 MatShellSetContext - sets the context for a shell matrix 763 764 Logically Collective on Mat 765 766 Input Parameters: 767 + mat - the shell matrix 768 - ctx - the context 769 770 Level: advanced 771 772 Fortran Notes: 773 To use this from Fortran you must write a Fortran interface definition for this 774 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 775 776 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 777 @*/ 778 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 779 { 780 Mat_Shell *shell = (Mat_Shell*)mat->data; 781 PetscErrorCode ierr; 782 PetscBool flg; 783 784 PetscFunctionBegin; 785 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 786 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 787 if (flg) { 788 shell->ctx = ctx; 789 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix"); 790 PetscFunctionReturn(0); 791 } 792 793 /*@ 794 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 795 after MatCreateShell() 796 797 Logically Collective on Mat 798 799 Input Parameter: 800 . mat - the shell matrix 801 802 Level: advanced 803 804 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation() 805 @*/ 806 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 807 { 808 PetscErrorCode ierr; 809 Mat_Shell *shell; 810 PetscBool flg; 811 812 PetscFunctionBegin; 813 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 814 ierr = PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);CHKERRQ(ierr); 815 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 816 shell = (Mat_Shell*)A->data; 817 shell->managescalingshifts = PETSC_FALSE; 818 A->ops->diagonalset = NULL; 819 A->ops->diagonalscale = NULL; 820 A->ops->scale = NULL; 821 A->ops->shift = NULL; 822 A->ops->axpy = NULL; 823 PetscFunctionReturn(0); 824 } 825 826 /*@C 827 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 828 829 Logically Collective on Mat 830 831 Input Parameters: 832 + mat - the shell matrix 833 . f - the function 834 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 835 - ctx - an optional context for the function 836 837 Output Parameter: 838 . flg - PETSC_TRUE if the multiply is likely correct 839 840 Options Database: 841 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 842 843 Level: advanced 844 845 Fortran Notes: 846 Not supported from Fortran 847 848 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose() 849 @*/ 850 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 851 { 852 PetscErrorCode ierr; 853 PetscInt m,n; 854 Mat mf,Dmf,Dmat,Ddiff; 855 PetscReal Diffnorm,Dmfnorm; 856 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 857 858 PetscFunctionBegin; 859 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 860 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr); 861 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 862 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr); 863 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 864 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 865 866 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 867 ierr = MatComputeExplicitOperator(mat,&Dmat);CHKERRQ(ierr); 868 869 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 870 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 871 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 872 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 873 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 874 flag = PETSC_FALSE; 875 if (v) { 876 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); 877 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 878 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 879 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 880 } 881 } else if (v) { 882 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 883 } 884 if (flg) *flg = flag; 885 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 886 ierr = MatDestroy(&mf);CHKERRQ(ierr); 887 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 888 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 889 PetscFunctionReturn(0); 890 } 891 892 /*@C 893 MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 894 895 Logically Collective on Mat 896 897 Input Parameters: 898 + mat - the shell matrix 899 . f - the function 900 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 901 - ctx - an optional context for the function 902 903 Output Parameter: 904 . flg - PETSC_TRUE if the multiply is likely correct 905 906 Options Database: 907 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 908 909 Level: advanced 910 911 Fortran Notes: 912 Not supported from Fortran 913 914 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult() 915 @*/ 916 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 917 { 918 PetscErrorCode ierr; 919 Vec x,y,z; 920 PetscInt m,n,M,N; 921 Mat mf,Dmf,Dmat,Ddiff; 922 PetscReal Diffnorm,Dmfnorm; 923 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 924 925 PetscFunctionBegin; 926 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 927 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr); 928 ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr); 929 ierr = VecDuplicate(y,&z);CHKERRQ(ierr); 930 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 931 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 932 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr); 933 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 934 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 935 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 936 ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr); 937 ierr = MatComputeExplicitOperatorTranspose(mat,&Dmat);CHKERRQ(ierr); 938 939 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 940 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 941 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 942 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 943 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 944 flag = PETSC_FALSE; 945 if (v) { 946 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); 947 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 948 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 949 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 950 } 951 } else if (v) { 952 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 953 } 954 if (flg) *flg = flag; 955 ierr = MatDestroy(&mf);CHKERRQ(ierr); 956 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 957 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 958 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 959 ierr = VecDestroy(&x);CHKERRQ(ierr); 960 ierr = VecDestroy(&y);CHKERRQ(ierr); 961 ierr = VecDestroy(&z);CHKERRQ(ierr); 962 PetscFunctionReturn(0); 963 } 964 965 /*@C 966 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 967 968 Logically Collective on Mat 969 970 Input Parameters: 971 + mat - the shell matrix 972 . op - the name of the operation 973 - f - the function that provides the operation. 974 975 Level: advanced 976 977 Usage: 978 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 979 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 980 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 981 982 Notes: 983 See the file include/petscmat.h for a complete list of matrix 984 operations, which all have the form MATOP_<OPERATION>, where 985 <OPERATION> is the name (in all capital letters) of the 986 user interface routine (e.g., MatMult() -> MATOP_MULT). 987 988 All user-provided functions (except for MATOP_DESTROY) should have the same calling 989 sequence as the usual matrix interface routines, since they 990 are intended to be accessed via the usual matrix interface 991 routines, e.g., 992 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 993 994 In particular each function MUST return an error code of 0 on success and 995 nonzero on failure. 996 997 Within each user-defined routine, the user should call 998 MatShellGetContext() to obtain the user-defined context that was 999 set by MatCreateShell(). 1000 1001 Fortran Notes: 1002 For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 1003 generate a matrix. See src/mat/examples/tests/ex120f.F 1004 1005 Use MatSetOperation() to set an operation for any matrix type 1006 1007 .keywords: matrix, shell, set, operation 1008 1009 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts() 1010 @*/ 1011 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 1012 { 1013 PetscBool flg; 1014 Mat_Shell *shell; 1015 PetscErrorCode ierr; 1016 1017 PetscFunctionBegin; 1018 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1019 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1020 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1021 shell = (Mat_Shell*)mat->data; 1022 1023 switch (op) { 1024 case MATOP_DESTROY: 1025 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1026 break; 1027 case MATOP_VIEW: 1028 if (!mat->ops->viewnative) { 1029 mat->ops->viewnative = mat->ops->view; 1030 } 1031 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1032 break; 1033 case MATOP_COPY: 1034 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1035 break; 1036 case MATOP_DIAGONAL_SET: 1037 case MATOP_DIAGONAL_SCALE: 1038 case MATOP_SHIFT: 1039 case MATOP_SCALE: 1040 case MATOP_AXPY: 1041 if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1042 (((void(**)(void))mat->ops)[op]) = f; 1043 break; 1044 case MATOP_GET_DIAGONAL: 1045 if (shell->managescalingshifts) { 1046 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1047 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1048 } else { 1049 shell->ops->getdiagonal = NULL; 1050 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1051 } 1052 break; 1053 case MATOP_MULT: 1054 if (shell->managescalingshifts) { 1055 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1056 mat->ops->mult = MatMult_Shell; 1057 } else { 1058 shell->ops->mult = NULL; 1059 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1060 } 1061 break; 1062 case MATOP_MULT_TRANSPOSE: 1063 if (shell->managescalingshifts) { 1064 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1065 mat->ops->multtranspose = MatMultTranspose_Shell; 1066 } else { 1067 shell->ops->multtranspose = NULL; 1068 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1069 } 1070 break; 1071 default: 1072 (((void(**)(void))mat->ops)[op]) = f; 1073 break; 1074 } 1075 PetscFunctionReturn(0); 1076 } 1077 1078 /*@C 1079 MatShellGetOperation - Gets a matrix function for a shell matrix. 1080 1081 Not Collective 1082 1083 Input Parameters: 1084 + mat - the shell matrix 1085 - op - the name of the operation 1086 1087 Output Parameter: 1088 . f - the function that provides the operation. 1089 1090 Level: advanced 1091 1092 Notes: 1093 See the file include/petscmat.h for a complete list of matrix 1094 operations, which all have the form MATOP_<OPERATION>, where 1095 <OPERATION> is the name (in all capital letters) of the 1096 user interface routine (e.g., MatMult() -> MATOP_MULT). 1097 1098 All user-provided functions have the same calling 1099 sequence as the usual matrix interface routines, since they 1100 are intended to be accessed via the usual matrix interface 1101 routines, e.g., 1102 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1103 1104 Within each user-defined routine, the user should call 1105 MatShellGetContext() to obtain the user-defined context that was 1106 set by MatCreateShell(). 1107 1108 .keywords: matrix, shell, set, operation 1109 1110 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 1111 @*/ 1112 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 1113 { 1114 PetscBool flg; 1115 Mat_Shell *shell; 1116 PetscErrorCode ierr; 1117 1118 PetscFunctionBegin; 1119 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1120 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1121 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1122 shell = (Mat_Shell*)mat->data; 1123 1124 switch (op) { 1125 case MATOP_DESTROY: 1126 *f = (void (*)(void))shell->ops->destroy; 1127 break; 1128 case MATOP_VIEW: 1129 *f = (void (*)(void))mat->ops->view; 1130 break; 1131 case MATOP_COPY: 1132 *f = (void (*)(void))shell->ops->copy; 1133 break; 1134 case MATOP_DIAGONAL_SET: 1135 case MATOP_DIAGONAL_SCALE: 1136 case MATOP_SHIFT: 1137 case MATOP_SCALE: 1138 case MATOP_AXPY: 1139 *f = (((void (**)(void))mat->ops)[op]); 1140 break; 1141 case MATOP_GET_DIAGONAL: 1142 if (shell->ops->getdiagonal) 1143 *f = (void (*)(void))shell->ops->getdiagonal; 1144 else 1145 *f = (((void (**)(void))mat->ops)[op]); 1146 break; 1147 case MATOP_MULT: 1148 if (shell->ops->mult) 1149 *f = (void (*)(void))shell->ops->mult; 1150 else 1151 *f = (((void (**)(void))mat->ops)[op]); 1152 break; 1153 case MATOP_MULT_TRANSPOSE: 1154 if (shell->ops->multtranspose) 1155 *f = (void (*)(void))shell->ops->multtranspose; 1156 else 1157 *f = (((void (**)(void))mat->ops)[op]); 1158 break; 1159 default: 1160 *f = (((void (**)(void))mat->ops)[op]); 1161 } 1162 PetscFunctionReturn(0); 1163 } 1164