1 2 /* 3 This provides a simple shell for Fortran (and C programmers) to 4 create a very simple matrix class for use with KSP without coding 5 much of anything. 6 */ 7 8 #include <petsc-private/matimpl.h> /*I "petscmat.h" I*/ 9 #include <petsc-private/vecimpl.h> 10 11 typedef struct { 12 PetscErrorCode (*destroy)(Mat); 13 PetscErrorCode (*mult)(Mat,Vec,Vec); 14 PetscErrorCode (*multtranspose)(Mat,Vec,Vec); 15 PetscErrorCode (*getdiagonal)(Mat,Vec); 16 17 PetscScalar vscale,vshift; 18 Vec dshift; 19 Vec left,right; 20 Vec dshift_owned,left_owned,right_owned; 21 Vec left_work,right_work; 22 Vec left_add_work,right_add_work; 23 PetscBool usingscaled; 24 void *ctx; 25 } Mat_Shell; 26 /* 27 The most general expression for the matrix is 28 29 S = L (a A + B) R 30 31 where 32 A is the matrix defined by the user's function 33 a is a scalar multiple 34 L is left scaling 35 R is right scaling 36 B is a diagonal shift defined by 37 diag(dshift) if the vector dshift is non-NULL 38 vscale*identity otherwise 39 40 The following identities apply: 41 42 Scale by c: 43 c [L (a A + B) R] = L [(a c) A + c B] R 44 45 Shift by c: 46 [L (a A + B) R] + c = L [a A + (B + c Linv Rinv)] R 47 48 Diagonal scaling is achieved by simply multiplying with existing L and R vectors 49 50 In the data structure: 51 52 vscale=1.0 means no special scaling will be applied 53 dshift=NULL means a constant diagonal shift (fall back to vshift) 54 vshift=0.0 means no constant diagonal shift, note that vshift is only used if dshift is NULL 55 */ 56 57 static PetscErrorCode MatMult_Shell(Mat,Vec,Vec); 58 static PetscErrorCode MatMultTranspose_Shell(Mat,Vec,Vec); 59 static PetscErrorCode MatGetDiagonal_Shell(Mat,Vec); 60 61 #undef __FUNCT__ 62 #define __FUNCT__ "MatShellUseScaledMethods" 63 static PetscErrorCode MatShellUseScaledMethods(Mat Y) 64 { 65 Mat_Shell *shell = (Mat_Shell*)Y->data; 66 67 PetscFunctionBegin; 68 if (shell->usingscaled) PetscFunctionReturn(0); 69 shell->mult = Y->ops->mult; 70 Y->ops->mult = MatMult_Shell; 71 if (Y->ops->multtranspose) { 72 shell->multtranspose = Y->ops->multtranspose; 73 Y->ops->multtranspose = MatMultTranspose_Shell; 74 } 75 if (Y->ops->getdiagonal) { 76 shell->getdiagonal = Y->ops->getdiagonal; 77 Y->ops->getdiagonal = MatGetDiagonal_Shell; 78 } 79 shell->usingscaled = PETSC_TRUE; 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "MatShellPreScaleLeft" 85 static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx) 86 { 87 Mat_Shell *shell = (Mat_Shell*)A->data; 88 PetscErrorCode ierr; 89 90 PetscFunctionBegin; 91 *xx = NULL; 92 if (!shell->left) { 93 *xx = x; 94 } else { 95 if (!shell->left_work) {ierr = VecDuplicate(shell->left,&shell->left_work);CHKERRQ(ierr);} 96 ierr = VecPointwiseMult(shell->left_work,x,shell->left);CHKERRQ(ierr); 97 *xx = shell->left_work; 98 } 99 PetscFunctionReturn(0); 100 } 101 102 #undef __FUNCT__ 103 #define __FUNCT__ "MatShellPreScaleRight" 104 static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx) 105 { 106 Mat_Shell *shell = (Mat_Shell*)A->data; 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 *xx = NULL; 111 if (!shell->right) { 112 *xx = x; 113 } else { 114 if (!shell->right_work) {ierr = VecDuplicate(shell->right,&shell->right_work);CHKERRQ(ierr);} 115 ierr = VecPointwiseMult(shell->right_work,x,shell->right);CHKERRQ(ierr); 116 *xx = shell->right_work; 117 } 118 PetscFunctionReturn(0); 119 } 120 121 #undef __FUNCT__ 122 #define __FUNCT__ "MatShellPostScaleLeft" 123 static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x) 124 { 125 Mat_Shell *shell = (Mat_Shell*)A->data; 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 if (shell->left) {ierr = VecPointwiseMult(x,x,shell->left);CHKERRQ(ierr);} 130 PetscFunctionReturn(0); 131 } 132 133 #undef __FUNCT__ 134 #define __FUNCT__ "MatShellPostScaleRight" 135 static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x) 136 { 137 Mat_Shell *shell = (Mat_Shell*)A->data; 138 PetscErrorCode ierr; 139 140 PetscFunctionBegin; 141 if (shell->right) {ierr = VecPointwiseMult(x,x,shell->right);CHKERRQ(ierr);} 142 PetscFunctionReturn(0); 143 } 144 145 #undef __FUNCT__ 146 #define __FUNCT__ "MatShellShiftAndScale" 147 static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y) 148 { 149 Mat_Shell *shell = (Mat_Shell*)A->data; 150 PetscErrorCode ierr; 151 152 PetscFunctionBegin; 153 if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */ 154 PetscInt i,m; 155 const PetscScalar *x,*d; 156 PetscScalar *y; 157 ierr = VecGetLocalSize(X,&m);CHKERRQ(ierr); 158 ierr = VecGetArrayRead(shell->dshift,&d);CHKERRQ(ierr); 159 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 160 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 161 for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i]; 162 ierr = VecRestoreArrayRead(shell->dshift,&d);CHKERRQ(ierr); 163 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 164 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 165 } else if (PetscAbsScalar(shell->vshift) != 0) { 166 ierr = VecAXPBY(Y,shell->vshift,shell->vscale,X);CHKERRQ(ierr); 167 } else if (shell->vscale != 1.0) { 168 ierr = VecScale(Y,shell->vscale);CHKERRQ(ierr); 169 } 170 PetscFunctionReturn(0); 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "MatShellGetContext" 175 /*@ 176 MatShellGetContext - Returns the user-provided context associated with a shell matrix. 177 178 Not Collective 179 180 Input Parameter: 181 . mat - the matrix, should have been created with MatCreateShell() 182 183 Output Parameter: 184 . ctx - the user provided context 185 186 Level: advanced 187 188 Notes: 189 This routine is intended for use within various shell matrix routines, 190 as set with MatShellSetOperation(). 191 192 .keywords: matrix, shell, get, context 193 194 .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext() 195 @*/ 196 PetscErrorCode MatShellGetContext(Mat mat,void *ctx) 197 { 198 PetscErrorCode ierr; 199 PetscBool flg; 200 201 PetscFunctionBegin; 202 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 203 PetscValidPointer(ctx,2); 204 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 205 if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx; 206 else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix"); 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "MatDestroy_Shell" 212 PetscErrorCode MatDestroy_Shell(Mat mat) 213 { 214 PetscErrorCode ierr; 215 Mat_Shell *shell = (Mat_Shell*)mat->data; 216 217 PetscFunctionBegin; 218 if (shell->destroy) { 219 ierr = (*shell->destroy)(mat);CHKERRQ(ierr); 220 } 221 ierr = VecDestroy(&shell->left_owned);CHKERRQ(ierr); 222 ierr = VecDestroy(&shell->right_owned);CHKERRQ(ierr); 223 ierr = VecDestroy(&shell->dshift_owned);CHKERRQ(ierr); 224 ierr = VecDestroy(&shell->left_work);CHKERRQ(ierr); 225 ierr = VecDestroy(&shell->right_work);CHKERRQ(ierr); 226 ierr = VecDestroy(&shell->left_add_work);CHKERRQ(ierr); 227 ierr = VecDestroy(&shell->right_add_work);CHKERRQ(ierr); 228 ierr = PetscFree(mat->data);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatMult_Shell" 234 PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y) 235 { 236 Mat_Shell *shell = (Mat_Shell*)A->data; 237 PetscErrorCode ierr; 238 Vec xx; 239 240 PetscFunctionBegin; 241 ierr = MatShellPreScaleRight(A,x,&xx);CHKERRQ(ierr); 242 ierr = (*shell->mult)(A,xx,y);CHKERRQ(ierr); 243 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 244 ierr = MatShellPostScaleLeft(A,y);CHKERRQ(ierr); 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "MatMultAdd_Shell" 250 PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z) 251 { 252 Mat_Shell *shell = (Mat_Shell*)A->data; 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 if (y == z) { 257 if (!shell->right_add_work) {ierr = VecDuplicate(z,&shell->right_add_work);CHKERRQ(ierr);} 258 ierr = MatMult(A,x,shell->right_add_work);CHKERRQ(ierr); 259 ierr = VecWAXPY(z,1.0,shell->right_add_work,y);CHKERRQ(ierr); 260 } else { 261 ierr = MatMult(A,x,z);CHKERRQ(ierr); 262 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 263 } 264 PetscFunctionReturn(0); 265 } 266 267 #undef __FUNCT__ 268 #define __FUNCT__ "MatMultTranspose_Shell" 269 PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y) 270 { 271 Mat_Shell *shell = (Mat_Shell*)A->data; 272 PetscErrorCode ierr; 273 Vec xx; 274 275 PetscFunctionBegin; 276 ierr = MatShellPreScaleLeft(A,x,&xx);CHKERRQ(ierr); 277 ierr = (*shell->multtranspose)(A,xx,y);CHKERRQ(ierr); 278 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 279 ierr = MatShellPostScaleRight(A,y);CHKERRQ(ierr); 280 PetscFunctionReturn(0); 281 } 282 283 #undef __FUNCT__ 284 #define __FUNCT__ "MatMultTransposeAdd_Shell" 285 PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z) 286 { 287 Mat_Shell *shell = (Mat_Shell*)A->data; 288 PetscErrorCode ierr; 289 290 PetscFunctionBegin; 291 if (y == z) { 292 if (!shell->left_add_work) {ierr = VecDuplicate(z,&shell->left_add_work);CHKERRQ(ierr);} 293 ierr = MatMultTranspose(A,x,shell->left_add_work);CHKERRQ(ierr); 294 ierr = VecWAXPY(z,1.0,shell->left_add_work,y);CHKERRQ(ierr); 295 } else { 296 ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr); 297 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 298 } 299 PetscFunctionReturn(0); 300 } 301 302 #undef __FUNCT__ 303 #define __FUNCT__ "MatGetDiagonal_Shell" 304 PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v) 305 { 306 Mat_Shell *shell = (Mat_Shell*)A->data; 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 ierr = (*shell->getdiagonal)(A,v);CHKERRQ(ierr); 311 ierr = VecScale(v,shell->vscale);CHKERRQ(ierr); 312 if (shell->dshift) { 313 ierr = VecPointwiseMult(v,v,shell->dshift);CHKERRQ(ierr); 314 } else { 315 ierr = VecShift(v,shell->vshift);CHKERRQ(ierr); 316 } 317 if (shell->left) {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);} 318 if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);} 319 PetscFunctionReturn(0); 320 } 321 322 #undef __FUNCT__ 323 #define __FUNCT__ "MatShift_Shell" 324 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 325 { 326 Mat_Shell *shell = (Mat_Shell*)Y->data; 327 PetscErrorCode ierr; 328 329 PetscFunctionBegin; 330 if (shell->left || shell->right || shell->dshift) { 331 if (!shell->dshift) { 332 if (!shell->dshift_owned) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift_owned);CHKERRQ(ierr);} 333 shell->dshift = shell->dshift_owned; 334 ierr = VecSet(shell->dshift,shell->vshift+a);CHKERRQ(ierr); 335 } else {ierr = VecScale(shell->dshift,a);CHKERRQ(ierr);} 336 if (shell->left) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 337 if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 338 } else shell->vshift += a; 339 ierr = MatShellUseScaledMethods(Y);CHKERRQ(ierr); 340 PetscFunctionReturn(0); 341 } 342 343 #undef __FUNCT__ 344 #define __FUNCT__ "MatScale_Shell" 345 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 346 { 347 Mat_Shell *shell = (Mat_Shell*)Y->data; 348 PetscErrorCode ierr; 349 350 PetscFunctionBegin; 351 shell->vscale *= a; 352 if (shell->dshift) { 353 ierr = VecScale(shell->dshift,a);CHKERRQ(ierr); 354 } else shell->vshift *= a; 355 ierr = MatShellUseScaledMethods(Y);CHKERRQ(ierr); 356 PetscFunctionReturn(0); 357 } 358 359 #undef __FUNCT__ 360 #define __FUNCT__ "MatDiagonalScale_Shell" 361 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right) 362 { 363 Mat_Shell *shell = (Mat_Shell*)Y->data; 364 PetscErrorCode ierr; 365 366 PetscFunctionBegin; 367 if (left) { 368 if (!shell->left_owned) {ierr = VecDuplicate(left,&shell->left_owned);CHKERRQ(ierr);} 369 if (shell->left) { 370 ierr = VecPointwiseMult(shell->left,shell->left,left);CHKERRQ(ierr); 371 } else { 372 shell->left = shell->left_owned; 373 ierr = VecCopy(left,shell->left);CHKERRQ(ierr); 374 } 375 } 376 if (right) { 377 if (!shell->right_owned) {ierr = VecDuplicate(right,&shell->right_owned);CHKERRQ(ierr);} 378 if (shell->right) { 379 ierr = VecPointwiseMult(shell->right,shell->right,right);CHKERRQ(ierr); 380 } else { 381 shell->right = shell->right_owned; 382 ierr = VecCopy(right,shell->right);CHKERRQ(ierr); 383 } 384 } 385 ierr = MatShellUseScaledMethods(Y);CHKERRQ(ierr); 386 PetscFunctionReturn(0); 387 } 388 389 #undef __FUNCT__ 390 #define __FUNCT__ "MatAssemblyEnd_Shell" 391 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 392 { 393 Mat_Shell *shell = (Mat_Shell*)Y->data; 394 395 PetscFunctionBegin; 396 if (t == MAT_FINAL_ASSEMBLY) { 397 shell->vshift = 0.0; 398 shell->vscale = 1.0; 399 shell->dshift = NULL; 400 shell->left = NULL; 401 shell->right = NULL; 402 if (shell->mult) { 403 Y->ops->mult = shell->mult; 404 shell->mult = NULL; 405 } 406 if (shell->multtranspose) { 407 Y->ops->multtranspose = shell->multtranspose; 408 shell->multtranspose = NULL; 409 } 410 if (shell->getdiagonal) { 411 Y->ops->getdiagonal = shell->getdiagonal; 412 shell->getdiagonal = NULL; 413 } 414 shell->usingscaled = PETSC_FALSE; 415 } 416 PetscFunctionReturn(0); 417 } 418 419 extern PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*); 420 421 static struct _MatOps MatOps_Values = {0, 422 0, 423 0, 424 0, 425 /* 4*/ 0, 426 0, 427 0, 428 0, 429 0, 430 0, 431 /*10*/ 0, 432 0, 433 0, 434 0, 435 0, 436 /*15*/ 0, 437 0, 438 0, 439 MatDiagonalScale_Shell, 440 0, 441 /*20*/ 0, 442 MatAssemblyEnd_Shell, 443 0, 444 0, 445 /*24*/ 0, 446 0, 447 0, 448 0, 449 0, 450 /*29*/ 0, 451 0, 452 0, 453 0, 454 0, 455 /*34*/ 0, 456 0, 457 0, 458 0, 459 0, 460 /*39*/ 0, 461 0, 462 0, 463 0, 464 0, 465 /*44*/ 0, 466 MatScale_Shell, 467 MatShift_Shell, 468 0, 469 0, 470 /*49*/ 0, 471 0, 472 0, 473 0, 474 0, 475 /*54*/ 0, 476 0, 477 0, 478 0, 479 0, 480 /*59*/ 0, 481 MatDestroy_Shell, 482 0, 483 0, 484 0, 485 /*64*/ 0, 486 0, 487 0, 488 0, 489 0, 490 /*69*/ 0, 491 0, 492 MatConvert_Shell, 493 0, 494 0, 495 /*74*/ 0, 496 0, 497 0, 498 0, 499 0, 500 /*79*/ 0, 501 0, 502 0, 503 0, 504 0, 505 /*84*/ 0, 506 0, 507 0, 508 0, 509 0, 510 /*89*/ 0, 511 0, 512 0, 513 0, 514 0, 515 /*94*/ 0, 516 0, 517 0, 518 0, 519 0, 520 /*99*/ 0, 521 0, 522 0, 523 0, 524 0, 525 /*104*/ 0, 526 0, 527 0, 528 0, 529 0, 530 /*109*/ 0, 531 0, 532 0, 533 0, 534 0, 535 /*114*/ 0, 536 0, 537 0, 538 0, 539 0, 540 /*119*/ 0, 541 0, 542 0, 543 0, 544 0, 545 /*124*/ 0, 546 0, 547 0, 548 0, 549 0, 550 /*129*/ 0, 551 0, 552 0, 553 0, 554 0, 555 /*134*/ 0, 556 0, 557 0, 558 0, 559 0, 560 /*139*/ 0, 561 0 562 }; 563 564 /*MC 565 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 566 567 Level: advanced 568 569 .seealso: MatCreateShell 570 M*/ 571 572 #undef __FUNCT__ 573 #define __FUNCT__ "MatCreate_Shell" 574 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 575 { 576 Mat_Shell *b; 577 PetscErrorCode ierr; 578 579 PetscFunctionBegin; 580 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 581 582 ierr = PetscNewLog(A,Mat_Shell,&b);CHKERRQ(ierr); 583 A->data = (void*)b; 584 585 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 586 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 587 588 b->ctx = 0; 589 b->vshift = 0.0; 590 b->vscale = 1.0; 591 b->mult = 0; 592 b->multtranspose = 0; 593 b->getdiagonal = 0; 594 A->assembled = PETSC_TRUE; 595 A->preallocated = PETSC_FALSE; 596 597 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 598 PetscFunctionReturn(0); 599 } 600 601 #undef __FUNCT__ 602 #define __FUNCT__ "MatCreateShell" 603 /*@C 604 MatCreateShell - Creates a new matrix class for use with a user-defined 605 private data storage format. 606 607 Collective on MPI_Comm 608 609 Input Parameters: 610 + comm - MPI communicator 611 . m - number of local rows (must be given) 612 . n - number of local columns (must be given) 613 . M - number of global rows (may be PETSC_DETERMINE) 614 . N - number of global columns (may be PETSC_DETERMINE) 615 - ctx - pointer to data needed by the shell matrix routines 616 617 Output Parameter: 618 . A - the matrix 619 620 Level: advanced 621 622 Usage: 623 $ extern int mult(Mat,Vec,Vec); 624 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 625 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 626 $ [ Use matrix for operations that have been set ] 627 $ MatDestroy(mat); 628 629 Notes: 630 The shell matrix type is intended to provide a simple class to use 631 with KSP (such as, for use with matrix-free methods). You should not 632 use the shell type if you plan to define a complete matrix class. 633 634 Fortran Notes: The context can only be an integer or a PetscObject 635 unfortunately it cannot be a Fortran array or derived type. 636 637 PETSc requires that matrices and vectors being used for certain 638 operations are partitioned accordingly. For example, when 639 creating a shell matrix, A, that supports parallel matrix-vector 640 products using MatMult(A,x,y) the user should set the number 641 of local matrix rows to be the number of local elements of the 642 corresponding result vector, y. Note that this is information is 643 required for use of the matrix interface routines, even though 644 the shell matrix may not actually be physically partitioned. 645 For example, 646 647 $ 648 $ Vec x, y 649 $ extern int mult(Mat,Vec,Vec); 650 $ Mat A 651 $ 652 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 653 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 654 $ VecGetLocalSize(y,&m); 655 $ VecGetLocalSize(x,&n); 656 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 657 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 658 $ MatMult(A,x,y); 659 $ MatDestroy(A); 660 $ VecDestroy(y); VecDestroy(x); 661 $ 662 663 .keywords: matrix, shell, create 664 665 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext() 666 @*/ 667 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 668 { 669 PetscErrorCode ierr; 670 671 PetscFunctionBegin; 672 ierr = MatCreate(comm,A);CHKERRQ(ierr); 673 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 674 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 675 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 676 ierr = MatSetUp(*A);CHKERRQ(ierr); 677 PetscFunctionReturn(0); 678 } 679 680 #undef __FUNCT__ 681 #define __FUNCT__ "MatShellSetContext" 682 /*@ 683 MatShellSetContext - sets the context for a shell matrix 684 685 Logically Collective on Mat 686 687 Input Parameters: 688 + mat - the shell matrix 689 - ctx - the context 690 691 Level: advanced 692 693 Fortran Notes: The context can only be an integer or a PetscObject 694 unfortunately it cannot be a Fortran array or derived type. 695 696 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 697 @*/ 698 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 699 { 700 Mat_Shell *shell = (Mat_Shell*)mat->data; 701 PetscErrorCode ierr; 702 PetscBool flg; 703 704 PetscFunctionBegin; 705 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 706 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 707 if (flg) { 708 shell->ctx = ctx; 709 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix"); 710 PetscFunctionReturn(0); 711 } 712 713 #undef __FUNCT__ 714 #define __FUNCT__ "MatShellSetOperation" 715 /*@C 716 MatShellSetOperation - Allows user to set a matrix operation for 717 a shell matrix. 718 719 Logically Collective on Mat 720 721 Input Parameters: 722 + mat - the shell matrix 723 . op - the name of the operation 724 - f - the function that provides the operation. 725 726 Level: advanced 727 728 Usage: 729 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 730 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 731 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 732 733 Notes: 734 See the file include/petscmat.h for a complete list of matrix 735 operations, which all have the form MATOP_<OPERATION>, where 736 <OPERATION> is the name (in all capital letters) of the 737 user interface routine (e.g., MatMult() -> MATOP_MULT). 738 739 All user-provided functions (execept for MATOP_DESTROY) should have the same calling 740 sequence as the usual matrix interface routines, since they 741 are intended to be accessed via the usual matrix interface 742 routines, e.g., 743 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 744 745 In particular each function MUST return an error code of 0 on success and 746 nonzero on failure. 747 748 Within each user-defined routine, the user should call 749 MatShellGetContext() to obtain the user-defined context that was 750 set by MatCreateShell(). 751 752 Fortran Notes: For MatGetVecs() the user code should check if the input left or right matrix is -1 and in that case not 753 generate a matrix. See src/mat/examples/tests/ex120f.F 754 755 .keywords: matrix, shell, set, operation 756 757 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext() 758 @*/ 759 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 760 { 761 PetscErrorCode ierr; 762 PetscBool flg; 763 764 PetscFunctionBegin; 765 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 766 switch (op) { 767 case MATOP_DESTROY: 768 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 769 if (flg) { 770 Mat_Shell *shell = (Mat_Shell*)mat->data; 771 shell->destroy = (PetscErrorCode (*)(Mat))f; 772 } else mat->ops->destroy = (PetscErrorCode (*)(Mat))f; 773 break; 774 case MATOP_VIEW: 775 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 776 break; 777 case MATOP_MULT: 778 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 779 if (!mat->ops->multadd) mat->ops->multadd = MatMultAdd_Shell; 780 break; 781 case MATOP_MULT_TRANSPOSE: 782 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 783 if (!mat->ops->multtransposeadd) mat->ops->multtransposeadd = MatMultTransposeAdd_Shell; 784 break; 785 default: 786 (((void(**)(void))mat->ops)[op]) = f; 787 } 788 PetscFunctionReturn(0); 789 } 790 791 #undef __FUNCT__ 792 #define __FUNCT__ "MatShellGetOperation" 793 /*@C 794 MatShellGetOperation - Gets a matrix function for a shell matrix. 795 796 Not Collective 797 798 Input Parameters: 799 + mat - the shell matrix 800 - op - the name of the operation 801 802 Output Parameter: 803 . f - the function that provides the operation. 804 805 Level: advanced 806 807 Notes: 808 See the file include/petscmat.h for a complete list of matrix 809 operations, which all have the form MATOP_<OPERATION>, where 810 <OPERATION> is the name (in all capital letters) of the 811 user interface routine (e.g., MatMult() -> MATOP_MULT). 812 813 All user-provided functions have the same calling 814 sequence as the usual matrix interface routines, since they 815 are intended to be accessed via the usual matrix interface 816 routines, e.g., 817 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 818 819 Within each user-defined routine, the user should call 820 MatShellGetContext() to obtain the user-defined context that was 821 set by MatCreateShell(). 822 823 .keywords: matrix, shell, set, operation 824 825 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 826 @*/ 827 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 828 { 829 PetscErrorCode ierr; 830 PetscBool flg; 831 832 PetscFunctionBegin; 833 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 834 if (op == MATOP_DESTROY) { 835 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 836 if (flg) { 837 Mat_Shell *shell = (Mat_Shell*)mat->data; 838 *f = (void (*)(void))shell->destroy; 839 } else { 840 *f = (void (*)(void))mat->ops->destroy; 841 } 842 } else if (op == MATOP_VIEW) { 843 *f = (void (*)(void))mat->ops->view; 844 } else { 845 *f = (((void (**)(void))mat->ops)[op]); 846 } 847 PetscFunctionReturn(0); 848 } 849 850