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