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