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