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 /* 3 */ PetscErrorCode (*mult)(Mat,Vec,Vec); 12 /* 5 */ PetscErrorCode (*multtranspose)(Mat,Vec,Vec); 13 /* 17 */ PetscErrorCode (*getdiagonal)(Mat,Vec); 14 /* 43 */ PetscErrorCode (*copy)(Mat,Mat,MatStructure); 15 /* 60 */ PetscErrorCode (*destroy)(Mat); 16 }; 17 18 typedef struct { 19 struct _MatShellOps ops[1]; 20 21 PetscScalar vscale,vshift; 22 Vec dshift; 23 Vec left,right; 24 Vec left_work,right_work; 25 Vec left_add_work,right_add_work; 26 Mat axpy; 27 PetscScalar axpy_vscale; 28 PetscBool managescalingshifts; /* The user will manage the scaling and shifts for the MATSHELL, not the default */ 29 void *ctx; 30 } Mat_Shell; 31 32 /* 33 xx = diag(left)*x 34 */ 35 static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx) 36 { 37 Mat_Shell *shell = (Mat_Shell*)A->data; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 *xx = NULL; 42 if (!shell->left) { 43 *xx = x; 44 } else { 45 if (!shell->left_work) {ierr = VecDuplicate(shell->left,&shell->left_work);CHKERRQ(ierr);} 46 ierr = VecPointwiseMult(shell->left_work,x,shell->left);CHKERRQ(ierr); 47 *xx = shell->left_work; 48 } 49 PetscFunctionReturn(0); 50 } 51 52 /* 53 xx = diag(right)*x 54 */ 55 static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx) 56 { 57 Mat_Shell *shell = (Mat_Shell*)A->data; 58 PetscErrorCode ierr; 59 60 PetscFunctionBegin; 61 *xx = NULL; 62 if (!shell->right) { 63 *xx = x; 64 } else { 65 if (!shell->right_work) {ierr = VecDuplicate(shell->right,&shell->right_work);CHKERRQ(ierr);} 66 ierr = VecPointwiseMult(shell->right_work,x,shell->right);CHKERRQ(ierr); 67 *xx = shell->right_work; 68 } 69 PetscFunctionReturn(0); 70 } 71 72 /* 73 x = diag(left)*x 74 */ 75 static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x) 76 { 77 Mat_Shell *shell = (Mat_Shell*)A->data; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 if (shell->left) {ierr = VecPointwiseMult(x,x,shell->left);CHKERRQ(ierr);} 82 PetscFunctionReturn(0); 83 } 84 85 /* 86 x = diag(right)*x 87 */ 88 static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x) 89 { 90 Mat_Shell *shell = (Mat_Shell*)A->data; 91 PetscErrorCode ierr; 92 93 PetscFunctionBegin; 94 if (shell->right) {ierr = VecPointwiseMult(x,x,shell->right);CHKERRQ(ierr);} 95 PetscFunctionReturn(0); 96 } 97 98 /* 99 Y = vscale*Y + diag(dshift)*X + vshift*X 100 101 On input Y already contains A*x 102 */ 103 static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y) 104 { 105 Mat_Shell *shell = (Mat_Shell*)A->data; 106 PetscErrorCode ierr; 107 108 PetscFunctionBegin; 109 if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */ 110 PetscInt i,m; 111 const PetscScalar *x,*d; 112 PetscScalar *y; 113 ierr = VecGetLocalSize(X,&m);CHKERRQ(ierr); 114 ierr = VecGetArrayRead(shell->dshift,&d);CHKERRQ(ierr); 115 ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); 116 ierr = VecGetArray(Y,&y);CHKERRQ(ierr); 117 for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i]; 118 ierr = VecRestoreArrayRead(shell->dshift,&d);CHKERRQ(ierr); 119 ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); 120 ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr); 121 } else { 122 ierr = VecScale(Y,shell->vscale);CHKERRQ(ierr); 123 } 124 if (shell->vshift != 0.0) {ierr = VecAXPY(Y,shell->vshift,X);CHKERRQ(ierr);} /* if test is for non-square matrices */ 125 PetscFunctionReturn(0); 126 } 127 128 /*@ 129 MatShellGetContext - Returns the user-provided context associated with a shell matrix. 130 131 Not Collective 132 133 Input Parameter: 134 . mat - the matrix, should have been created with MatCreateShell() 135 136 Output Parameter: 137 . ctx - the user provided context 138 139 Level: advanced 140 141 Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this 142 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 143 144 .keywords: matrix, shell, get, context 145 146 .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext() 147 @*/ 148 PetscErrorCode MatShellGetContext(Mat mat,void *ctx) 149 { 150 PetscErrorCode ierr; 151 PetscBool flg; 152 153 PetscFunctionBegin; 154 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 155 PetscValidPointer(ctx,2); 156 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 157 if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx; 158 else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix"); 159 PetscFunctionReturn(0); 160 } 161 162 PetscErrorCode MatDestroy_Shell(Mat mat) 163 { 164 PetscErrorCode ierr; 165 Mat_Shell *shell = (Mat_Shell*)mat->data; 166 167 PetscFunctionBegin; 168 if (shell->ops->destroy) { 169 ierr = (*shell->ops->destroy)(mat);CHKERRQ(ierr); 170 } 171 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 172 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 173 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 174 ierr = VecDestroy(&shell->left_work);CHKERRQ(ierr); 175 ierr = VecDestroy(&shell->right_work);CHKERRQ(ierr); 176 ierr = VecDestroy(&shell->left_add_work);CHKERRQ(ierr); 177 ierr = VecDestroy(&shell->right_add_work);CHKERRQ(ierr); 178 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 179 ierr = PetscFree(mat->data);CHKERRQ(ierr); 180 PetscFunctionReturn(0); 181 } 182 183 PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str) 184 { 185 Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data; 186 PetscErrorCode ierr; 187 PetscBool matflg; 188 189 PetscFunctionBegin; 190 ierr = PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);CHKERRQ(ierr); 191 if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL"); 192 193 ierr = PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 194 ierr = PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));CHKERRQ(ierr); 195 196 if (shellA->ops->copy) { 197 ierr = (*shellA->ops->copy)(A,B,str);CHKERRQ(ierr); 198 } 199 shellB->vscale = shellA->vscale; 200 shellB->vshift = shellA->vshift; 201 if (shellA->dshift) { 202 if (!shellB->dshift) { 203 ierr = VecDuplicate(shellA->dshift,&shellB->dshift);CHKERRQ(ierr); 204 } 205 ierr = VecCopy(shellA->dshift,shellB->dshift);CHKERRQ(ierr); 206 } else { 207 ierr = VecDestroy(&shellB->dshift);CHKERRQ(ierr); 208 } 209 if (shellA->left) { 210 if (!shellB->left) { 211 ierr = VecDuplicate(shellA->left,&shellB->left);CHKERRQ(ierr); 212 } 213 ierr = VecCopy(shellA->left,shellB->left);CHKERRQ(ierr); 214 } else { 215 ierr = VecDestroy(&shellB->left);CHKERRQ(ierr); 216 } 217 if (shellA->right) { 218 if (!shellB->right) { 219 ierr = VecDuplicate(shellA->right,&shellB->right);CHKERRQ(ierr); 220 } 221 ierr = VecCopy(shellA->right,shellB->right);CHKERRQ(ierr); 222 } else { 223 ierr = VecDestroy(&shellB->right);CHKERRQ(ierr); 224 } 225 ierr = MatDestroy(&shellB->axpy);CHKERRQ(ierr); 226 if (shellA->axpy) { 227 ierr = PetscObjectReference((PetscObject)shellA->axpy);CHKERRQ(ierr); 228 shellB->axpy = shellA->axpy; 229 shellB->axpy_vscale = shellA->axpy_vscale; 230 } 231 PetscFunctionReturn(0); 232 } 233 234 PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M) 235 { 236 PetscErrorCode ierr; 237 void *ctx; 238 239 PetscFunctionBegin; 240 ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr); 241 ierr = MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);CHKERRQ(ierr); 242 ierr = MatCopy(mat,*M,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 243 PetscFunctionReturn(0); 244 } 245 246 PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y) 247 { 248 Mat_Shell *shell = (Mat_Shell*)A->data; 249 PetscErrorCode ierr; 250 Vec xx; 251 PetscObjectState instate,outstate; 252 253 PetscFunctionBegin; 254 ierr = MatShellPreScaleRight(A,x,&xx);CHKERRQ(ierr); 255 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 256 if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL"); 257 ierr = (*shell->ops->mult)(A,xx,y);CHKERRQ(ierr); 258 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 259 if (instate == outstate) { 260 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 261 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 262 } 263 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 264 ierr = MatShellPostScaleLeft(A,y);CHKERRQ(ierr); 265 266 if (shell->axpy) { 267 if (!shell->left_work) {ierr = MatCreateVecs(A,&shell->left_work,NULL);CHKERRQ(ierr);} 268 ierr = MatMult(shell->axpy,x,shell->left_work);CHKERRQ(ierr); 269 ierr = VecAXPY(y,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr); 270 } 271 PetscFunctionReturn(0); 272 } 273 274 PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z) 275 { 276 Mat_Shell *shell = (Mat_Shell*)A->data; 277 PetscErrorCode ierr; 278 279 PetscFunctionBegin; 280 if (y == z) { 281 if (!shell->right_add_work) {ierr = VecDuplicate(z,&shell->right_add_work);CHKERRQ(ierr);} 282 ierr = MatMult(A,x,shell->right_add_work);CHKERRQ(ierr); 283 ierr = VecAXPY(z,1.0,shell->right_add_work);CHKERRQ(ierr); 284 } else { 285 ierr = MatMult(A,x,z);CHKERRQ(ierr); 286 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 287 } 288 PetscFunctionReturn(0); 289 } 290 291 PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y) 292 { 293 Mat_Shell *shell = (Mat_Shell*)A->data; 294 PetscErrorCode ierr; 295 Vec xx; 296 PetscObjectState instate,outstate; 297 298 PetscFunctionBegin; 299 ierr = MatShellPreScaleLeft(A,x,&xx);CHKERRQ(ierr); 300 ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr); 301 if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL"); 302 ierr = (*shell->ops->multtranspose)(A,xx,y);CHKERRQ(ierr); 303 ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr); 304 if (instate == outstate) { 305 /* increase the state of the output vector since the user did not update its state themself as should have been done */ 306 ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr); 307 } 308 ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr); 309 ierr = MatShellPostScaleRight(A,y);CHKERRQ(ierr); 310 PetscFunctionReturn(0); 311 } 312 313 PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z) 314 { 315 Mat_Shell *shell = (Mat_Shell*)A->data; 316 PetscErrorCode ierr; 317 318 PetscFunctionBegin; 319 if (y == z) { 320 if (!shell->left_add_work) {ierr = VecDuplicate(z,&shell->left_add_work);CHKERRQ(ierr);} 321 ierr = MatMultTranspose(A,x,shell->left_add_work);CHKERRQ(ierr); 322 ierr = VecWAXPY(z,1.0,shell->left_add_work,y);CHKERRQ(ierr); 323 } else { 324 ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr); 325 ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr); 326 } 327 PetscFunctionReturn(0); 328 } 329 330 /* 331 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 332 */ 333 PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v) 334 { 335 Mat_Shell *shell = (Mat_Shell*)A->data; 336 PetscErrorCode ierr; 337 338 PetscFunctionBegin; 339 if (shell->ops->getdiagonal) { 340 ierr = (*shell->ops->getdiagonal)(A,v);CHKERRQ(ierr); 341 } else { 342 ierr = VecSet(v,0.0);CHKERRQ(ierr); 343 } 344 ierr = VecScale(v,shell->vscale);CHKERRQ(ierr); 345 if (shell->dshift) { 346 ierr = VecAXPY(v,1.0,shell->dshift);CHKERRQ(ierr); 347 } 348 ierr = VecShift(v,shell->vshift);CHKERRQ(ierr); 349 if (shell->left) {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);} 350 if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);} 351 if (shell->axpy) { 352 if (!shell->left_work) {ierr = VecDuplicate(v,&shell->left_work);CHKERRQ(ierr);} 353 ierr = MatGetDiagonal(shell->axpy,shell->left_work);CHKERRQ(ierr); 354 ierr = VecAXPY(v,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr); 355 } 356 PetscFunctionReturn(0); 357 } 358 359 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 360 { 361 Mat_Shell *shell = (Mat_Shell*)Y->data; 362 PetscErrorCode ierr; 363 364 PetscFunctionBegin; 365 if (shell->left || shell->right) { 366 if (!shell->dshift) { 367 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 368 ierr = VecSet(shell->dshift,a);CHKERRQ(ierr); 369 } else { 370 if (shell->left) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 371 if (shell->right) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 372 ierr = VecShift(shell->dshift,a);CHKERRQ(ierr); 373 } 374 if (shell->left) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 375 if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 376 } else shell->vshift += a; 377 PetscFunctionReturn(0); 378 } 379 380 PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins) 381 { 382 Mat_Shell *shell = (Mat_Shell*)A->data; 383 PetscErrorCode ierr; 384 385 PetscFunctionBegin; 386 if (ins == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported with INSERT_VALUES"); 387 if (!shell->dshift) {ierr = VecDuplicate(D,&shell->dshift);CHKERRQ(ierr);} 388 if (shell->left || shell->right) { 389 if (!shell->right_work) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);CHKERRQ(ierr);} 390 if (shell->left && shell->right) { 391 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 392 ierr = VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);CHKERRQ(ierr); 393 } else if (shell->left) { 394 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 395 } else { 396 ierr = VecPointwiseDivide(shell->right_work,D,shell->right);CHKERRQ(ierr); 397 } 398 if (!shell->dshift) { 399 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 400 ierr = VecCopy(shell->dshift,shell->right_work);CHKERRQ(ierr); 401 } else { 402 ierr = VecAXPY(shell->dshift,1.0,shell->right_work);CHKERRQ(ierr); 403 } 404 } else { 405 ierr = VecAXPY(shell->dshift,1.0,D);CHKERRQ(ierr); 406 } 407 PetscFunctionReturn(0); 408 } 409 410 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 411 { 412 Mat_Shell *shell = (Mat_Shell*)Y->data; 413 PetscErrorCode ierr; 414 415 PetscFunctionBegin; 416 shell->vscale *= a; 417 shell->vshift *= a; 418 if (shell->dshift) { 419 ierr = VecScale(shell->dshift,a);CHKERRQ(ierr); 420 } 421 shell->axpy_vscale *= a; 422 PetscFunctionReturn(0); 423 } 424 425 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right) 426 { 427 Mat_Shell *shell = (Mat_Shell*)Y->data; 428 PetscErrorCode ierr; 429 430 PetscFunctionBegin; 431 if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation"); 432 if (left) { 433 if (!shell->left) { 434 ierr = VecDuplicate(left,&shell->left);CHKERRQ(ierr); 435 ierr = VecCopy(left,shell->left);CHKERRQ(ierr); 436 } else { 437 ierr = VecPointwiseMult(shell->left,shell->left,left);CHKERRQ(ierr); 438 } 439 } 440 if (right) { 441 if (!shell->right) { 442 ierr = VecDuplicate(right,&shell->right);CHKERRQ(ierr); 443 ierr = VecCopy(right,shell->right);CHKERRQ(ierr); 444 } else { 445 ierr = VecPointwiseMult(shell->right,shell->right,right);CHKERRQ(ierr); 446 } 447 } 448 PetscFunctionReturn(0); 449 } 450 451 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 452 { 453 Mat_Shell *shell = (Mat_Shell*)Y->data; 454 PetscErrorCode ierr; 455 456 PetscFunctionBegin; 457 if (t == MAT_FINAL_ASSEMBLY) { 458 shell->vshift = 0.0; 459 shell->vscale = 1.0; 460 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 461 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 462 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 463 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 464 } 465 PetscFunctionReturn(0); 466 } 467 468 PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*); 469 470 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d) 471 { 472 PetscFunctionBegin; 473 *missing = PETSC_FALSE; 474 PetscFunctionReturn(0); 475 } 476 477 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str) 478 { 479 Mat_Shell *shell = (Mat_Shell*)Y->data; 480 PetscErrorCode ierr; 481 482 PetscFunctionBegin; 483 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 484 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 485 shell->axpy = X; 486 shell->axpy_vscale = a; 487 PetscFunctionReturn(0); 488 } 489 490 static struct _MatOps MatOps_Values = {0, 491 0, 492 0, 493 0, 494 /* 4*/ MatMultAdd_Shell, 495 0, 496 MatMultTransposeAdd_Shell, 497 0, 498 0, 499 0, 500 /*10*/ 0, 501 0, 502 0, 503 0, 504 0, 505 /*15*/ 0, 506 0, 507 MatGetDiagonal_Shell, 508 MatDiagonalScale_Shell, 509 0, 510 /*20*/ 0, 511 MatAssemblyEnd_Shell, 512 0, 513 0, 514 /*24*/ 0, 515 0, 516 0, 517 0, 518 0, 519 /*29*/ 0, 520 0, 521 0, 522 0, 523 0, 524 /*34*/ MatDuplicate_Shell, 525 0, 526 0, 527 0, 528 0, 529 /*39*/ MatAXPY_Shell, 530 0, 531 0, 532 0, 533 MatCopy_Shell, 534 /*44*/ 0, 535 MatScale_Shell, 536 MatShift_Shell, 537 MatDiagonalSet_Shell, 538 0, 539 /*49*/ 0, 540 0, 541 0, 542 0, 543 0, 544 /*54*/ 0, 545 0, 546 0, 547 0, 548 0, 549 /*59*/ 0, 550 MatDestroy_Shell, 551 0, 552 0, 553 0, 554 /*64*/ 0, 555 0, 556 0, 557 0, 558 0, 559 /*69*/ 0, 560 0, 561 MatConvert_Shell, 562 0, 563 0, 564 /*74*/ 0, 565 0, 566 0, 567 0, 568 0, 569 /*79*/ 0, 570 0, 571 0, 572 0, 573 0, 574 /*84*/ 0, 575 0, 576 0, 577 0, 578 0, 579 /*89*/ 0, 580 0, 581 0, 582 0, 583 0, 584 /*94*/ 0, 585 0, 586 0, 587 0, 588 0, 589 /*99*/ 0, 590 0, 591 0, 592 0, 593 0, 594 /*104*/ 0, 595 0, 596 0, 597 0, 598 0, 599 /*109*/ 0, 600 0, 601 0, 602 0, 603 MatMissingDiagonal_Shell, 604 /*114*/ 0, 605 0, 606 0, 607 0, 608 0, 609 /*119*/ 0, 610 0, 611 0, 612 0, 613 0, 614 /*124*/ 0, 615 0, 616 0, 617 0, 618 0, 619 /*129*/ 0, 620 0, 621 0, 622 0, 623 0, 624 /*134*/ 0, 625 0, 626 0, 627 0, 628 0, 629 /*139*/ 0, 630 0, 631 0 632 }; 633 634 /*MC 635 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 636 637 Level: advanced 638 639 .seealso: MatCreateShell() 640 M*/ 641 642 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 643 { 644 Mat_Shell *b; 645 PetscErrorCode ierr; 646 647 PetscFunctionBegin; 648 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 649 650 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 651 A->data = (void*)b; 652 653 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 654 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 655 656 b->ctx = 0; 657 b->vshift = 0.0; 658 b->vscale = 1.0; 659 b->managescalingshifts = PETSC_TRUE; 660 A->assembled = PETSC_TRUE; 661 A->preallocated = PETSC_FALSE; 662 663 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 664 PetscFunctionReturn(0); 665 } 666 667 /*@C 668 MatCreateShell - Creates a new matrix class for use with a user-defined 669 private data storage format. 670 671 Collective on MPI_Comm 672 673 Input Parameters: 674 + comm - MPI communicator 675 . m - number of local rows (must be given) 676 . n - number of local columns (must be given) 677 . M - number of global rows (may be PETSC_DETERMINE) 678 . N - number of global columns (may be PETSC_DETERMINE) 679 - ctx - pointer to data needed by the shell matrix routines 680 681 Output Parameter: 682 . A - the matrix 683 684 Level: advanced 685 686 Usage: 687 $ extern int mult(Mat,Vec,Vec); 688 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 689 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 690 $ [ Use matrix for operations that have been set ] 691 $ MatDestroy(mat); 692 693 Notes: 694 The shell matrix type is intended to provide a simple class to use 695 with KSP (such as, for use with matrix-free methods). You should not 696 use the shell type if you plan to define a complete matrix class. 697 698 Fortran Notes: To use this from Fortran with a ctx you must write an interface definition for this 699 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 700 in as the ctx argument. 701 702 PETSc requires that matrices and vectors being used for certain 703 operations are partitioned accordingly. For example, when 704 creating a shell matrix, A, that supports parallel matrix-vector 705 products using MatMult(A,x,y) the user should set the number 706 of local matrix rows to be the number of local elements of the 707 corresponding result vector, y. Note that this is information is 708 required for use of the matrix interface routines, even though 709 the shell matrix may not actually be physically partitioned. 710 For example, 711 712 $ 713 $ Vec x, y 714 $ extern int mult(Mat,Vec,Vec); 715 $ Mat A 716 $ 717 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 718 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 719 $ VecGetLocalSize(y,&m); 720 $ VecGetLocalSize(x,&n); 721 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 722 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 723 $ MatMult(A,x,y); 724 $ MatDestroy(A); 725 $ VecDestroy(y); VecDestroy(x); 726 $ 727 728 729 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), and MatScale() internally so these 730 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 731 732 733 For rectangular matrices do all the scalings and shifts make sense? 734 735 Developers Notes: Regarding shifting and scaling. The general form is 736 737 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 738 739 The order you apply the operations is important. For example if you have a dshift then 740 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 741 you get s*vscale*A + diag(shift) 742 743 A is the user provided function. 744 745 .keywords: matrix, shell, create 746 747 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts() 748 @*/ 749 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 750 { 751 PetscErrorCode ierr; 752 753 PetscFunctionBegin; 754 ierr = MatCreate(comm,A);CHKERRQ(ierr); 755 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 756 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 757 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 758 ierr = MatSetUp(*A);CHKERRQ(ierr); 759 PetscFunctionReturn(0); 760 } 761 762 /*@ 763 MatShellSetContext - sets the context for a shell matrix 764 765 Logically Collective on Mat 766 767 Input Parameters: 768 + mat - the shell matrix 769 - ctx - the context 770 771 Level: advanced 772 773 Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this 774 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 775 776 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 777 @*/ 778 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 779 { 780 Mat_Shell *shell = (Mat_Shell*)mat->data; 781 PetscErrorCode ierr; 782 PetscBool flg; 783 784 PetscFunctionBegin; 785 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 786 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 787 if (flg) { 788 shell->ctx = ctx; 789 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix"); 790 PetscFunctionReturn(0); 791 } 792 793 /*@ 794 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 795 after MatCreateShell() 796 797 Logically Collective on Mat 798 799 Input Parameter: 800 . mat - the shell matrix 801 802 Level: advanced 803 804 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation() 805 @*/ 806 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 807 { 808 PetscErrorCode ierr; 809 Mat_Shell *shell; 810 PetscBool flg; 811 812 PetscFunctionBegin; 813 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 814 ierr = PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);CHKERRQ(ierr); 815 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 816 shell = (Mat_Shell*)A->data; 817 shell->managescalingshifts = PETSC_FALSE; 818 A->ops->diagonalset = NULL; 819 A->ops->diagonalscale = NULL; 820 A->ops->scale = NULL; 821 A->ops->shift = NULL; 822 A->ops->axpy = NULL; 823 PetscFunctionReturn(0); 824 } 825 826 /*@C 827 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 828 829 Logically Collective on Mat 830 831 Input Parameters: 832 + mat - the shell matrix 833 . f - the function 834 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 835 - ctx - an optional context for the function 836 837 Output Parameter: 838 . flg - PETSC_TRUE if the multiply is likely correct 839 840 Options Database: 841 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 842 843 Level: advanced 844 845 Fortran Notes: Not supported from Fortran 846 847 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose() 848 @*/ 849 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 850 { 851 PetscErrorCode ierr; 852 PetscInt m,n; 853 Mat mf,Dmf,Dmat,Ddiff; 854 PetscReal Diffnorm,Dmfnorm; 855 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 856 857 PetscFunctionBegin; 858 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 859 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr); 860 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 861 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr); 862 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 863 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 864 865 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 866 ierr = MatComputeExplicitOperator(mat,&Dmat);CHKERRQ(ierr); 867 868 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 869 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 870 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 871 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 872 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 873 flag = PETSC_FALSE; 874 if (v) { 875 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n",(double)(Diffnorm/Dmfnorm));CHKERRQ(ierr); 876 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 877 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 878 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 879 } 880 } else if (v) { 881 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 882 } 883 if (flg) *flg = flag; 884 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 885 ierr = MatDestroy(&mf);CHKERRQ(ierr); 886 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 887 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 888 PetscFunctionReturn(0); 889 } 890 891 /*@C 892 MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 893 894 Logically Collective on Mat 895 896 Input Parameters: 897 + mat - the shell matrix 898 . f - the function 899 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 900 - ctx - an optional context for the function 901 902 Output Parameter: 903 . flg - PETSC_TRUE if the multiply is likely correct 904 905 Options Database: 906 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 907 908 Level: advanced 909 910 Fortran Notes: Not supported from Fortran 911 912 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult() 913 @*/ 914 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 915 { 916 PetscErrorCode ierr; 917 Vec x,y,z; 918 PetscInt m,n,M,N; 919 Mat mf,Dmf,Dmat,Ddiff; 920 PetscReal Diffnorm,Dmfnorm; 921 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 922 923 PetscFunctionBegin; 924 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 925 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr); 926 ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr); 927 ierr = VecDuplicate(y,&z);CHKERRQ(ierr); 928 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 929 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 930 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr); 931 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 932 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 933 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 934 ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr); 935 ierr = MatComputeExplicitOperatorTranspose(mat,&Dmat);CHKERRQ(ierr); 936 937 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 938 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 939 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 940 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 941 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 942 flag = PETSC_FALSE; 943 if (v) { 944 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n",(double)(Diffnorm/Dmfnorm));CHKERRQ(ierr); 945 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 946 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 947 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 948 } 949 } else if (v) { 950 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 951 } 952 if (flg) *flg = flag; 953 ierr = MatDestroy(&mf);CHKERRQ(ierr); 954 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 955 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 956 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 957 ierr = VecDestroy(&x);CHKERRQ(ierr); 958 ierr = VecDestroy(&y);CHKERRQ(ierr); 959 ierr = VecDestroy(&z);CHKERRQ(ierr); 960 PetscFunctionReturn(0); 961 } 962 963 /*@C 964 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 965 966 Logically Collective on Mat 967 968 Input Parameters: 969 + mat - the shell matrix 970 . op - the name of the operation 971 - f - the function that provides the operation. 972 973 Level: advanced 974 975 Usage: 976 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 977 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 978 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 979 980 Notes: 981 See the file include/petscmat.h for a complete list of matrix 982 operations, which all have the form MATOP_<OPERATION>, where 983 <OPERATION> is the name (in all capital letters) of the 984 user interface routine (e.g., MatMult() -> MATOP_MULT). 985 986 All user-provided functions (except for MATOP_DESTROY) should have the same calling 987 sequence as the usual matrix interface routines, since they 988 are intended to be accessed via the usual matrix interface 989 routines, e.g., 990 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 991 992 In particular each function MUST return an error code of 0 on success and 993 nonzero on failure. 994 995 Within each user-defined routine, the user should call 996 MatShellGetContext() to obtain the user-defined context that was 997 set by MatCreateShell(). 998 999 Fortran Notes: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 1000 generate a matrix. See src/mat/examples/tests/ex120f.F 1001 1002 Use MatSetOperation() to set an operation for any matrix type 1003 1004 .keywords: matrix, shell, set, operation 1005 1006 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts() 1007 @*/ 1008 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 1009 { 1010 PetscBool flg; 1011 Mat_Shell *shell; 1012 PetscErrorCode ierr; 1013 1014 PetscFunctionBegin; 1015 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1016 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1017 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1018 shell = (Mat_Shell*)mat->data; 1019 1020 switch (op) { 1021 case MATOP_DESTROY: 1022 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1023 break; 1024 case MATOP_VIEW: 1025 if (!mat->ops->viewnative) { 1026 mat->ops->viewnative = mat->ops->view; 1027 } 1028 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1029 break; 1030 case MATOP_COPY: 1031 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1032 break; 1033 case MATOP_DIAGONAL_SET: 1034 case MATOP_DIAGONAL_SCALE: 1035 case MATOP_SHIFT: 1036 case MATOP_SCALE: 1037 case MATOP_AXPY: 1038 if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1039 (((void(**)(void))mat->ops)[op]) = f; 1040 break; 1041 case MATOP_GET_DIAGONAL: 1042 if (shell->managescalingshifts) { 1043 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1044 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1045 } else { 1046 shell->ops->getdiagonal = NULL; 1047 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1048 } 1049 break; 1050 case MATOP_MULT: 1051 if (shell->managescalingshifts) { 1052 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1053 mat->ops->mult = MatMult_Shell; 1054 } else { 1055 shell->ops->mult = NULL; 1056 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1057 } 1058 break; 1059 case MATOP_MULT_TRANSPOSE: 1060 if (shell->managescalingshifts) { 1061 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1062 mat->ops->multtranspose = MatMultTranspose_Shell; 1063 } else { 1064 shell->ops->multtranspose = NULL; 1065 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1066 } 1067 break; 1068 default: 1069 (((void(**)(void))mat->ops)[op]) = f; 1070 break; 1071 } 1072 PetscFunctionReturn(0); 1073 } 1074 1075 /*@C 1076 MatShellGetOperation - Gets a matrix function for a shell matrix. 1077 1078 Not Collective 1079 1080 Input Parameters: 1081 + mat - the shell matrix 1082 - op - the name of the operation 1083 1084 Output Parameter: 1085 . f - the function that provides the operation. 1086 1087 Level: advanced 1088 1089 Notes: 1090 See the file include/petscmat.h for a complete list of matrix 1091 operations, which all have the form MATOP_<OPERATION>, where 1092 <OPERATION> is the name (in all capital letters) of the 1093 user interface routine (e.g., MatMult() -> MATOP_MULT). 1094 1095 All user-provided functions have the same calling 1096 sequence as the usual matrix interface routines, since they 1097 are intended to be accessed via the usual matrix interface 1098 routines, e.g., 1099 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1100 1101 Within each user-defined routine, the user should call 1102 MatShellGetContext() to obtain the user-defined context that was 1103 set by MatCreateShell(). 1104 1105 .keywords: matrix, shell, set, operation 1106 1107 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 1108 @*/ 1109 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 1110 { 1111 PetscBool flg; 1112 Mat_Shell *shell; 1113 PetscErrorCode ierr; 1114 1115 PetscFunctionBegin; 1116 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1117 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1118 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1119 shell = (Mat_Shell*)mat->data; 1120 1121 switch (op) { 1122 case MATOP_DESTROY: 1123 *f = (void (*)(void))shell->ops->destroy; 1124 break; 1125 case MATOP_VIEW: 1126 *f = (void (*)(void))mat->ops->view; 1127 break; 1128 case MATOP_COPY: 1129 *f = (void (*)(void))shell->ops->copy; 1130 break; 1131 case MATOP_DIAGONAL_SET: 1132 case MATOP_DIAGONAL_SCALE: 1133 case MATOP_SHIFT: 1134 case MATOP_SCALE: 1135 case MATOP_AXPY: 1136 *f = (((void (**)(void))mat->ops)[op]); 1137 break; 1138 case MATOP_GET_DIAGONAL: 1139 if (shell->ops->getdiagonal) 1140 *f = (void (*)(void))shell->ops->getdiagonal; 1141 else 1142 *f = (((void (**)(void))mat->ops)[op]); 1143 break; 1144 case MATOP_MULT: 1145 if (shell->ops->mult) 1146 *f = (void (*)(void))shell->ops->mult; 1147 else 1148 *f = (((void (**)(void))mat->ops)[op]); 1149 break; 1150 case MATOP_MULT_TRANSPOSE: 1151 if (shell->ops->multtranspose) 1152 *f = (void (*)(void))shell->ops->multtranspose; 1153 else 1154 *f = (((void (**)(void))mat->ops)[op]); 1155 break; 1156 default: 1157 *f = (((void (**)(void))mat->ops)[op]); 1158 } 1159 PetscFunctionReturn(0); 1160 } 1161