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