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