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 static struct _MatOps MatOps_Values = {0, 421 0, 422 0, 423 0, 424 /* 4*/ 0, 425 0, 426 0, 427 0, 428 0, 429 0, 430 /*10*/ 0, 431 0, 432 0, 433 0, 434 0, 435 /*15*/ 0, 436 0, 437 0, 438 MatDiagonalScale_Shell, 439 0, 440 /*20*/ 0, 441 MatAssemblyEnd_Shell, 442 0, 443 0, 444 /*24*/ 0, 445 0, 446 0, 447 0, 448 0, 449 /*29*/ 0, 450 0, 451 0, 452 0, 453 0, 454 /*34*/ 0, 455 0, 456 0, 457 0, 458 0, 459 /*39*/ 0, 460 0, 461 0, 462 0, 463 0, 464 /*44*/ 0, 465 MatScale_Shell, 466 MatShift_Shell, 467 0, 468 0, 469 /*49*/ 0, 470 0, 471 0, 472 0, 473 0, 474 /*54*/ 0, 475 0, 476 0, 477 0, 478 0, 479 /*59*/ 0, 480 MatDestroy_Shell, 481 0, 482 0, 483 0, 484 /*64*/ 0, 485 0, 486 0, 487 0, 488 0, 489 /*69*/ 0, 490 0, 491 MatConvert_Shell, 492 0, 493 0, 494 /*74*/ 0, 495 0, 496 0, 497 0, 498 0, 499 /*79*/ 0, 500 0, 501 0, 502 0, 503 0, 504 /*84*/ 0, 505 0, 506 0, 507 0, 508 0, 509 /*89*/ 0, 510 0, 511 0, 512 0, 513 0, 514 /*94*/ 0, 515 0, 516 0, 517 0, 518 0, 519 /*99*/ 0, 520 0, 521 0, 522 0, 523 0, 524 /*104*/ 0, 525 0, 526 0, 527 0, 528 0, 529 /*109*/ 0, 530 0, 531 0, 532 0, 533 0, 534 /*114*/ 0, 535 0, 536 0, 537 0, 538 0, 539 /*119*/ 0, 540 0, 541 0, 542 0, 543 0, 544 /*124*/ 0, 545 0, 546 0, 547 0, 548 0, 549 /*129*/ 0, 550 0, 551 0, 552 0, 553 0, 554 /*134*/ 0, 555 0, 556 0, 557 0, 558 0, 559 /*139*/ 0, 560 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,&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: To use this from Fortran with a ctx you must write an interface definition for this 635 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 636 in as the ctx argument. 637 638 PETSc requires that matrices and vectors being used for certain 639 operations are partitioned accordingly. For example, when 640 creating a shell matrix, A, that supports parallel matrix-vector 641 products using MatMult(A,x,y) the user should set the number 642 of local matrix rows to be the number of local elements of the 643 corresponding result vector, y. Note that this is information is 644 required for use of the matrix interface routines, even though 645 the shell matrix may not actually be physically partitioned. 646 For example, 647 648 $ 649 $ Vec x, y 650 $ extern int mult(Mat,Vec,Vec); 651 $ Mat A 652 $ 653 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 654 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 655 $ VecGetLocalSize(y,&m); 656 $ VecGetLocalSize(x,&n); 657 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 658 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 659 $ MatMult(A,x,y); 660 $ MatDestroy(A); 661 $ VecDestroy(y); VecDestroy(x); 662 $ 663 664 .keywords: matrix, shell, create 665 666 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext() 667 @*/ 668 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 669 { 670 PetscErrorCode ierr; 671 672 PetscFunctionBegin; 673 ierr = MatCreate(comm,A);CHKERRQ(ierr); 674 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 675 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 676 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 677 ierr = MatSetUp(*A);CHKERRQ(ierr); 678 PetscFunctionReturn(0); 679 } 680 681 #undef __FUNCT__ 682 #define __FUNCT__ "MatShellSetContext" 683 /*@ 684 MatShellSetContext - sets the context for a shell matrix 685 686 Logically Collective on Mat 687 688 Input Parameters: 689 + mat - the shell matrix 690 - ctx - the context 691 692 Level: advanced 693 694 Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this 695 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 696 697 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 698 @*/ 699 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 700 { 701 Mat_Shell *shell = (Mat_Shell*)mat->data; 702 PetscErrorCode ierr; 703 PetscBool flg; 704 705 PetscFunctionBegin; 706 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 707 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 708 if (flg) { 709 shell->ctx = ctx; 710 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix"); 711 PetscFunctionReturn(0); 712 } 713 714 #undef __FUNCT__ 715 #define __FUNCT__ "MatShellSetOperation" 716 /*@C 717 MatShellSetOperation - Allows user to set a matrix operation for 718 a shell matrix. 719 720 Logically Collective on Mat 721 722 Input Parameters: 723 + mat - the shell matrix 724 . op - the name of the operation 725 - f - the function that provides the operation. 726 727 Level: advanced 728 729 Usage: 730 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 731 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 732 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 733 734 Notes: 735 See the file include/petscmat.h for a complete list of matrix 736 operations, which all have the form MATOP_<OPERATION>, where 737 <OPERATION> is the name (in all capital letters) of the 738 user interface routine (e.g., MatMult() -> MATOP_MULT). 739 740 All user-provided functions (execept for MATOP_DESTROY) should have the same calling 741 sequence as the usual matrix interface routines, since they 742 are intended to be accessed via the usual matrix interface 743 routines, e.g., 744 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 745 746 In particular each function MUST return an error code of 0 on success and 747 nonzero on failure. 748 749 Within each user-defined routine, the user should call 750 MatShellGetContext() to obtain the user-defined context that was 751 set by MatCreateShell(). 752 753 Fortran Notes: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 754 generate a matrix. See src/mat/examples/tests/ex120f.F 755 756 .keywords: matrix, shell, set, operation 757 758 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext() 759 @*/ 760 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 761 { 762 PetscErrorCode ierr; 763 PetscBool flg; 764 765 PetscFunctionBegin; 766 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 767 switch (op) { 768 case MATOP_DESTROY: 769 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 770 if (flg) { 771 Mat_Shell *shell = (Mat_Shell*)mat->data; 772 shell->destroy = (PetscErrorCode (*)(Mat))f; 773 } else mat->ops->destroy = (PetscErrorCode (*)(Mat))f; 774 break; 775 case MATOP_VIEW: 776 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 777 break; 778 case MATOP_MULT: 779 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 780 if (!mat->ops->multadd) mat->ops->multadd = MatMultAdd_Shell; 781 break; 782 case MATOP_MULT_TRANSPOSE: 783 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 784 if (!mat->ops->multtransposeadd) mat->ops->multtransposeadd = MatMultTransposeAdd_Shell; 785 break; 786 default: 787 (((void(**)(void))mat->ops)[op]) = f; 788 } 789 PetscFunctionReturn(0); 790 } 791 792 #undef __FUNCT__ 793 #define __FUNCT__ "MatShellGetOperation" 794 /*@C 795 MatShellGetOperation - Gets a matrix function for a shell matrix. 796 797 Not Collective 798 799 Input Parameters: 800 + mat - the shell matrix 801 - op - the name of the operation 802 803 Output Parameter: 804 . f - the function that provides the operation. 805 806 Level: advanced 807 808 Notes: 809 See the file include/petscmat.h for a complete list of matrix 810 operations, which all have the form MATOP_<OPERATION>, where 811 <OPERATION> is the name (in all capital letters) of the 812 user interface routine (e.g., MatMult() -> MATOP_MULT). 813 814 All user-provided functions have the same calling 815 sequence as the usual matrix interface routines, since they 816 are intended to be accessed via the usual matrix interface 817 routines, e.g., 818 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 819 820 Within each user-defined routine, the user should call 821 MatShellGetContext() to obtain the user-defined context that was 822 set by MatCreateShell(). 823 824 .keywords: matrix, shell, set, operation 825 826 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 827 @*/ 828 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 829 { 830 PetscErrorCode ierr; 831 PetscBool flg; 832 833 PetscFunctionBegin; 834 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 835 if (op == MATOP_DESTROY) { 836 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 837 if (flg) { 838 Mat_Shell *shell = (Mat_Shell*)mat->data; 839 *f = (void (*)(void))shell->destroy; 840 } else { 841 *f = (void (*)(void))mat->ops->destroy; 842 } 843 } else if (op == MATOP_VIEW) { 844 *f = (void (*)(void))mat->ops->view; 845 } else { 846 *f = (((void (**)(void))mat->ops)[op]); 847 } 848 PetscFunctionReturn(0); 849 } 850 851