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