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