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 0 563 }; 564 565 /*MC 566 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 567 568 Level: advanced 569 570 .seealso: MatCreateShell 571 M*/ 572 573 #undef __FUNCT__ 574 #define __FUNCT__ "MatCreate_Shell" 575 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 576 { 577 Mat_Shell *b; 578 PetscErrorCode ierr; 579 580 PetscFunctionBegin; 581 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 582 583 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 584 A->data = (void*)b; 585 586 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 587 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 588 589 b->ctx = 0; 590 b->vshift = 0.0; 591 b->vscale = 1.0; 592 b->mult = 0; 593 b->multtranspose = 0; 594 b->getdiagonal = 0; 595 A->assembled = PETSC_TRUE; 596 A->preallocated = PETSC_FALSE; 597 598 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 599 PetscFunctionReturn(0); 600 } 601 602 #undef __FUNCT__ 603 #define __FUNCT__ "MatCreateShell" 604 /*@C 605 MatCreateShell - Creates a new matrix class for use with a user-defined 606 private data storage format. 607 608 Collective on MPI_Comm 609 610 Input Parameters: 611 + comm - MPI communicator 612 . m - number of local rows (must be given) 613 . n - number of local columns (must be given) 614 . M - number of global rows (may be PETSC_DETERMINE) 615 . N - number of global columns (may be PETSC_DETERMINE) 616 - ctx - pointer to data needed by the shell matrix routines 617 618 Output Parameter: 619 . A - the matrix 620 621 Level: advanced 622 623 Usage: 624 $ extern int mult(Mat,Vec,Vec); 625 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 626 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 627 $ [ Use matrix for operations that have been set ] 628 $ MatDestroy(mat); 629 630 Notes: 631 The shell matrix type is intended to provide a simple class to use 632 with KSP (such as, for use with matrix-free methods). You should not 633 use the shell type if you plan to define a complete matrix class. 634 635 Fortran Notes: The context can only be an integer or a PetscObject 636 unfortunately it cannot be a Fortran array or derived type. 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: The context can only be an integer or a PetscObject 695 unfortunately it cannot be a Fortran array or derived type. 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