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