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