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