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