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 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)"); 342 ierr = VecScale(v,shell->vscale);CHKERRQ(ierr); 343 if (shell->dshift) { 344 ierr = VecAXPY(v,1.0,shell->dshift);CHKERRQ(ierr); 345 } 346 ierr = VecShift(v,shell->vshift);CHKERRQ(ierr); 347 if (shell->left) {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);} 348 if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);} 349 if (shell->axpy) { 350 if (!shell->left_work) {ierr = VecDuplicate(v,&shell->left_work);CHKERRQ(ierr);} 351 ierr = MatGetDiagonal(shell->axpy,shell->left_work);CHKERRQ(ierr); 352 ierr = VecAXPY(v,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr); 353 } 354 PetscFunctionReturn(0); 355 } 356 357 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 358 { 359 Mat_Shell *shell = (Mat_Shell*)Y->data; 360 PetscErrorCode ierr; 361 362 PetscFunctionBegin; 363 if (shell->left || shell->right) { 364 if (!shell->dshift) { 365 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 366 ierr = VecSet(shell->dshift,a);CHKERRQ(ierr); 367 } else { 368 if (shell->left) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 369 if (shell->right) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 370 ierr = VecShift(shell->dshift,a);CHKERRQ(ierr); 371 } 372 if (shell->left) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);} 373 if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);} 374 } else shell->vshift += a; 375 PetscFunctionReturn(0); 376 } 377 378 PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins) 379 { 380 Mat_Shell *shell = (Mat_Shell*)A->data; 381 PetscErrorCode ierr; 382 383 PetscFunctionBegin; 384 if (ins == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported with INSERT_VALUES"); 385 if (!shell->dshift) {ierr = VecDuplicate(D,&shell->dshift);CHKERRQ(ierr);} 386 if (shell->left || shell->right) { 387 if (!shell->right_work) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);CHKERRQ(ierr);} 388 if (shell->left && shell->right) { 389 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 390 ierr = VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);CHKERRQ(ierr); 391 } else if (shell->left) { 392 ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr); 393 } else { 394 ierr = VecPointwiseDivide(shell->right_work,D,shell->right);CHKERRQ(ierr); 395 } 396 if (!shell->dshift) { 397 ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr); 398 ierr = VecCopy(shell->dshift,shell->right_work);CHKERRQ(ierr); 399 } else { 400 ierr = VecAXPY(shell->dshift,1.0,shell->right_work);CHKERRQ(ierr); 401 } 402 } else { 403 ierr = VecAXPY(shell->dshift,1.0,D);CHKERRQ(ierr); 404 } 405 PetscFunctionReturn(0); 406 } 407 408 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 409 { 410 Mat_Shell *shell = (Mat_Shell*)Y->data; 411 PetscErrorCode ierr; 412 413 PetscFunctionBegin; 414 shell->vscale *= a; 415 shell->vshift *= a; 416 if (shell->dshift) { 417 ierr = VecScale(shell->dshift,a);CHKERRQ(ierr); 418 } 419 shell->axpy_vscale *= a; 420 PetscFunctionReturn(0); 421 } 422 423 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right) 424 { 425 Mat_Shell *shell = (Mat_Shell*)Y->data; 426 PetscErrorCode ierr; 427 428 PetscFunctionBegin; 429 if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation"); 430 if (left) { 431 if (!shell->left) { 432 ierr = VecDuplicate(left,&shell->left);CHKERRQ(ierr); 433 ierr = VecCopy(left,shell->left);CHKERRQ(ierr); 434 } else { 435 ierr = VecPointwiseMult(shell->left,shell->left,left);CHKERRQ(ierr); 436 } 437 } 438 if (right) { 439 if (!shell->right) { 440 ierr = VecDuplicate(right,&shell->right);CHKERRQ(ierr); 441 ierr = VecCopy(right,shell->right);CHKERRQ(ierr); 442 } else { 443 ierr = VecPointwiseMult(shell->right,shell->right,right);CHKERRQ(ierr); 444 } 445 } 446 PetscFunctionReturn(0); 447 } 448 449 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 450 { 451 Mat_Shell *shell = (Mat_Shell*)Y->data; 452 PetscErrorCode ierr; 453 454 PetscFunctionBegin; 455 if (t == MAT_FINAL_ASSEMBLY) { 456 shell->vshift = 0.0; 457 shell->vscale = 1.0; 458 ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr); 459 ierr = VecDestroy(&shell->left);CHKERRQ(ierr); 460 ierr = VecDestroy(&shell->right);CHKERRQ(ierr); 461 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 462 } 463 PetscFunctionReturn(0); 464 } 465 466 PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*); 467 468 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d) 469 { 470 PetscFunctionBegin; 471 *missing = PETSC_FALSE; 472 PetscFunctionReturn(0); 473 } 474 475 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str) 476 { 477 Mat_Shell *shell = (Mat_Shell*)Y->data; 478 PetscErrorCode ierr; 479 480 PetscFunctionBegin; 481 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 482 ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr); 483 shell->axpy = X; 484 shell->axpy_vscale = a; 485 PetscFunctionReturn(0); 486 } 487 488 static struct _MatOps MatOps_Values = {0, 489 0, 490 0, 491 0, 492 /* 4*/ MatMultAdd_Shell, 493 0, 494 MatMultTransposeAdd_Shell, 495 0, 496 0, 497 0, 498 /*10*/ 0, 499 0, 500 0, 501 0, 502 0, 503 /*15*/ 0, 504 0, 505 MatGetDiagonal_Shell, 506 MatDiagonalScale_Shell, 507 0, 508 /*20*/ 0, 509 MatAssemblyEnd_Shell, 510 0, 511 0, 512 /*24*/ 0, 513 0, 514 0, 515 0, 516 0, 517 /*29*/ 0, 518 0, 519 0, 520 0, 521 0, 522 /*34*/ MatDuplicate_Shell, 523 0, 524 0, 525 0, 526 0, 527 /*39*/ MatAXPY_Shell, 528 0, 529 0, 530 0, 531 MatCopy_Shell, 532 /*44*/ 0, 533 MatScale_Shell, 534 MatShift_Shell, 535 MatDiagonalSet_Shell, 536 0, 537 /*49*/ 0, 538 0, 539 0, 540 0, 541 0, 542 /*54*/ 0, 543 0, 544 0, 545 0, 546 0, 547 /*59*/ 0, 548 MatDestroy_Shell, 549 0, 550 0, 551 0, 552 /*64*/ 0, 553 0, 554 0, 555 0, 556 0, 557 /*69*/ 0, 558 0, 559 MatConvert_Shell, 560 0, 561 0, 562 /*74*/ 0, 563 0, 564 0, 565 0, 566 0, 567 /*79*/ 0, 568 0, 569 0, 570 0, 571 0, 572 /*84*/ 0, 573 0, 574 0, 575 0, 576 0, 577 /*89*/ 0, 578 0, 579 0, 580 0, 581 0, 582 /*94*/ 0, 583 0, 584 0, 585 0, 586 0, 587 /*99*/ 0, 588 0, 589 0, 590 0, 591 0, 592 /*104*/ 0, 593 0, 594 0, 595 0, 596 0, 597 /*109*/ 0, 598 0, 599 0, 600 0, 601 MatMissingDiagonal_Shell, 602 /*114*/ 0, 603 0, 604 0, 605 0, 606 0, 607 /*119*/ 0, 608 0, 609 0, 610 0, 611 0, 612 /*124*/ 0, 613 0, 614 0, 615 0, 616 0, 617 /*129*/ 0, 618 0, 619 0, 620 0, 621 0, 622 /*134*/ 0, 623 0, 624 0, 625 0, 626 0, 627 /*139*/ 0, 628 0, 629 0 630 }; 631 632 /*MC 633 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 634 635 Level: advanced 636 637 .seealso: MatCreateShell() 638 M*/ 639 640 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A) 641 { 642 Mat_Shell *b; 643 PetscErrorCode ierr; 644 645 PetscFunctionBegin; 646 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 647 648 ierr = PetscNewLog(A,&b);CHKERRQ(ierr); 649 A->data = (void*)b; 650 651 ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr); 652 ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr); 653 654 b->ctx = 0; 655 b->vshift = 0.0; 656 b->vscale = 1.0; 657 b->managescalingshifts = PETSC_TRUE; 658 A->assembled = PETSC_TRUE; 659 A->preallocated = PETSC_FALSE; 660 661 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 662 PetscFunctionReturn(0); 663 } 664 665 /*@C 666 MatCreateShell - Creates a new matrix class for use with a user-defined 667 private data storage format. 668 669 Collective on MPI_Comm 670 671 Input Parameters: 672 + comm - MPI communicator 673 . m - number of local rows (must be given) 674 . n - number of local columns (must be given) 675 . M - number of global rows (may be PETSC_DETERMINE) 676 . N - number of global columns (may be PETSC_DETERMINE) 677 - ctx - pointer to data needed by the shell matrix routines 678 679 Output Parameter: 680 . A - the matrix 681 682 Level: advanced 683 684 Usage: 685 $ extern int mult(Mat,Vec,Vec); 686 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 687 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 688 $ [ Use matrix for operations that have been set ] 689 $ MatDestroy(mat); 690 691 Notes: 692 The shell matrix type is intended to provide a simple class to use 693 with KSP (such as, for use with matrix-free methods). You should not 694 use the shell type if you plan to define a complete matrix class. 695 696 Fortran Notes: To use this from Fortran with a ctx you must write an interface definition for this 697 function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing 698 in as the ctx argument. 699 700 PETSc requires that matrices and vectors being used for certain 701 operations are partitioned accordingly. For example, when 702 creating a shell matrix, A, that supports parallel matrix-vector 703 products using MatMult(A,x,y) the user should set the number 704 of local matrix rows to be the number of local elements of the 705 corresponding result vector, y. Note that this is information is 706 required for use of the matrix interface routines, even though 707 the shell matrix may not actually be physically partitioned. 708 For example, 709 710 $ 711 $ Vec x, y 712 $ extern int mult(Mat,Vec,Vec); 713 $ Mat A 714 $ 715 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 716 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 717 $ VecGetLocalSize(y,&m); 718 $ VecGetLocalSize(x,&n); 719 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 720 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 721 $ MatMult(A,x,y); 722 $ MatDestroy(A); 723 $ VecDestroy(y); VecDestroy(x); 724 $ 725 726 727 MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), and MatScale() internally so these 728 operations cannot be overwritten unless MatShellSetManageScalingShifts() is called. 729 730 731 For rectangular matrices do all the scalings and shifts make sense? 732 733 Developers Notes: Regarding shifting and scaling. The general form is 734 735 diag(left)(vscale*A + diag(dshift) + vshift I)diag(right) 736 737 The order you apply the operations is important. For example if you have a dshift then 738 apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift 739 you get s*vscale*A + diag(shift) 740 741 A is the user provided function. 742 743 .keywords: matrix, shell, create 744 745 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts() 746 @*/ 747 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 748 { 749 PetscErrorCode ierr; 750 751 PetscFunctionBegin; 752 ierr = MatCreate(comm,A);CHKERRQ(ierr); 753 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 754 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 755 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 756 ierr = MatSetUp(*A);CHKERRQ(ierr); 757 PetscFunctionReturn(0); 758 } 759 760 /*@ 761 MatShellSetContext - sets the context for a shell matrix 762 763 Logically Collective on Mat 764 765 Input Parameters: 766 + mat - the shell matrix 767 - ctx - the context 768 769 Level: advanced 770 771 Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this 772 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 773 774 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 775 @*/ 776 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 777 { 778 Mat_Shell *shell = (Mat_Shell*)mat->data; 779 PetscErrorCode ierr; 780 PetscBool flg; 781 782 PetscFunctionBegin; 783 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 784 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 785 if (flg) { 786 shell->ctx = ctx; 787 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix"); 788 PetscFunctionReturn(0); 789 } 790 791 /*@ 792 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 793 after MatCreateShell() 794 795 Logically Collective on Mat 796 797 Input Parameter: 798 . mat - the shell matrix 799 800 Level: advanced 801 802 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation() 803 @*/ 804 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 805 { 806 PetscErrorCode ierr; 807 Mat_Shell *shell; 808 PetscBool flg; 809 810 PetscFunctionBegin; 811 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 812 ierr = PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);CHKERRQ(ierr); 813 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 814 shell = (Mat_Shell*)A->data; 815 shell->managescalingshifts = PETSC_FALSE; 816 A->ops->diagonalset = NULL; 817 A->ops->diagonalscale = NULL; 818 A->ops->scale = NULL; 819 A->ops->shift = NULL; 820 A->ops->axpy = NULL; 821 PetscFunctionReturn(0); 822 } 823 824 /*@C 825 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 826 827 Logically Collective on Mat 828 829 Input Parameters: 830 + mat - the shell matrix 831 . f - the function 832 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 833 - ctx - an optional context for the function 834 835 Output Parameter: 836 . flg - PETSC_TRUE if the multiply is likely correct 837 838 Options Database: 839 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 840 841 Level: advanced 842 843 Fortran Notes: Not supported from Fortran 844 845 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose() 846 @*/ 847 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 848 { 849 PetscErrorCode ierr; 850 PetscInt m,n; 851 Mat mf,Dmf,Dmat,Ddiff; 852 PetscReal Diffnorm,Dmfnorm; 853 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 854 855 PetscFunctionBegin; 856 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 857 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr); 858 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 859 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr); 860 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 861 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 862 863 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 864 ierr = MatComputeExplicitOperator(mat,&Dmat);CHKERRQ(ierr); 865 866 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 867 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 868 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 869 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 870 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 871 flag = PETSC_FALSE; 872 if (v) { 873 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); 874 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 875 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 876 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 877 } 878 } else if (v) { 879 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 880 } 881 if (flg) *flg = flag; 882 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 883 ierr = MatDestroy(&mf);CHKERRQ(ierr); 884 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 885 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 886 PetscFunctionReturn(0); 887 } 888 889 /*@C 890 MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 891 892 Logically Collective on Mat 893 894 Input Parameters: 895 + mat - the shell matrix 896 . f - the function 897 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 898 - ctx - an optional context for the function 899 900 Output Parameter: 901 . flg - PETSC_TRUE if the multiply is likely correct 902 903 Options Database: 904 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 905 906 Level: advanced 907 908 Fortran Notes: Not supported from Fortran 909 910 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult() 911 @*/ 912 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 913 { 914 PetscErrorCode ierr; 915 Vec x,y,z; 916 PetscInt m,n,M,N; 917 Mat mf,Dmf,Dmat,Ddiff; 918 PetscReal Diffnorm,Dmfnorm; 919 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 920 921 PetscFunctionBegin; 922 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 923 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr); 924 ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr); 925 ierr = VecDuplicate(y,&z);CHKERRQ(ierr); 926 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 927 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 928 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr); 929 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 930 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 931 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 932 ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr); 933 ierr = MatComputeExplicitOperatorTranspose(mat,&Dmat);CHKERRQ(ierr); 934 935 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 936 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 937 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 938 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 939 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 940 flag = PETSC_FALSE; 941 if (v) { 942 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); 943 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 944 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 945 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 946 } 947 } else if (v) { 948 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 949 } 950 if (flg) *flg = flag; 951 ierr = MatDestroy(&mf);CHKERRQ(ierr); 952 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 953 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 954 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 955 ierr = VecDestroy(&x);CHKERRQ(ierr); 956 ierr = VecDestroy(&y);CHKERRQ(ierr); 957 ierr = VecDestroy(&z);CHKERRQ(ierr); 958 PetscFunctionReturn(0); 959 } 960 961 /*@C 962 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 963 964 Logically Collective on Mat 965 966 Input Parameters: 967 + mat - the shell matrix 968 . op - the name of the operation 969 - f - the function that provides the operation. 970 971 Level: advanced 972 973 Usage: 974 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 975 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 976 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 977 978 Notes: 979 See the file include/petscmat.h for a complete list of matrix 980 operations, which all have the form MATOP_<OPERATION>, where 981 <OPERATION> is the name (in all capital letters) of the 982 user interface routine (e.g., MatMult() -> MATOP_MULT). 983 984 All user-provided functions (except for MATOP_DESTROY) should have the same calling 985 sequence as the usual matrix interface routines, since they 986 are intended to be accessed via the usual matrix interface 987 routines, e.g., 988 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 989 990 In particular each function MUST return an error code of 0 on success and 991 nonzero on failure. 992 993 Within each user-defined routine, the user should call 994 MatShellGetContext() to obtain the user-defined context that was 995 set by MatCreateShell(). 996 997 Fortran Notes: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 998 generate a matrix. See src/mat/examples/tests/ex120f.F 999 1000 Use MatSetOperation() to set an operation for any matrix type 1001 1002 .keywords: matrix, shell, set, operation 1003 1004 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts() 1005 @*/ 1006 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 1007 { 1008 PetscBool flg; 1009 Mat_Shell *shell; 1010 PetscErrorCode ierr; 1011 1012 PetscFunctionBegin; 1013 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1014 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1015 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1016 shell = (Mat_Shell*)mat->data; 1017 1018 switch (op) { 1019 case MATOP_DESTROY: 1020 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1021 break; 1022 case MATOP_VIEW: 1023 if (!mat->ops->viewnative) { 1024 mat->ops->viewnative = mat->ops->view; 1025 } 1026 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1027 break; 1028 case MATOP_COPY: 1029 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1030 break; 1031 case MATOP_DIAGONAL_SET: 1032 case MATOP_DIAGONAL_SCALE: 1033 case MATOP_SHIFT: 1034 case MATOP_SCALE: 1035 case MATOP_AXPY: 1036 if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1037 (((void(**)(void))mat->ops)[op]) = f; 1038 break; 1039 case MATOP_GET_DIAGONAL: 1040 if (shell->managescalingshifts) { 1041 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1042 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1043 } else { 1044 shell->ops->getdiagonal = NULL; 1045 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1046 } 1047 break; 1048 case MATOP_MULT: 1049 if (shell->managescalingshifts) { 1050 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1051 mat->ops->mult = MatMult_Shell; 1052 } else { 1053 shell->ops->mult = NULL; 1054 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1055 } 1056 break; 1057 case MATOP_MULT_TRANSPOSE: 1058 if (shell->managescalingshifts) { 1059 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1060 mat->ops->multtranspose = MatMultTranspose_Shell; 1061 } else { 1062 shell->ops->multtranspose = NULL; 1063 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1064 } 1065 break; 1066 default: 1067 (((void(**)(void))mat->ops)[op]) = f; 1068 break; 1069 } 1070 PetscFunctionReturn(0); 1071 } 1072 1073 /*@C 1074 MatShellGetOperation - Gets a matrix function for a shell matrix. 1075 1076 Not Collective 1077 1078 Input Parameters: 1079 + mat - the shell matrix 1080 - op - the name of the operation 1081 1082 Output Parameter: 1083 . f - the function that provides the operation. 1084 1085 Level: advanced 1086 1087 Notes: 1088 See the file include/petscmat.h for a complete list of matrix 1089 operations, which all have the form MATOP_<OPERATION>, where 1090 <OPERATION> is the name (in all capital letters) of the 1091 user interface routine (e.g., MatMult() -> MATOP_MULT). 1092 1093 All user-provided functions have the same calling 1094 sequence as the usual matrix interface routines, since they 1095 are intended to be accessed via the usual matrix interface 1096 routines, e.g., 1097 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1098 1099 Within each user-defined routine, the user should call 1100 MatShellGetContext() to obtain the user-defined context that was 1101 set by MatCreateShell(). 1102 1103 .keywords: matrix, shell, set, operation 1104 1105 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 1106 @*/ 1107 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 1108 { 1109 PetscBool flg; 1110 Mat_Shell *shell; 1111 PetscErrorCode ierr; 1112 1113 PetscFunctionBegin; 1114 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1115 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1116 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1117 shell = (Mat_Shell*)mat->data; 1118 1119 switch (op) { 1120 case MATOP_DESTROY: 1121 *f = (void (*)(void))shell->ops->destroy; 1122 break; 1123 case MATOP_VIEW: 1124 *f = (void (*)(void))mat->ops->view; 1125 break; 1126 case MATOP_COPY: 1127 *f = (void (*)(void))shell->ops->copy; 1128 break; 1129 case MATOP_DIAGONAL_SET: 1130 case MATOP_DIAGONAL_SCALE: 1131 case MATOP_SHIFT: 1132 case MATOP_SCALE: 1133 case MATOP_AXPY: 1134 *f = (((void (**)(void))mat->ops)[op]); 1135 break; 1136 case MATOP_GET_DIAGONAL: 1137 if (shell->ops->getdiagonal) 1138 *f = (void (*)(void))shell->ops->getdiagonal; 1139 else 1140 *f = (((void (**)(void))mat->ops)[op]); 1141 break; 1142 case MATOP_MULT: 1143 if (shell->ops->mult) 1144 *f = (void (*)(void))shell->ops->mult; 1145 else 1146 *f = (((void (**)(void))mat->ops)[op]); 1147 break; 1148 case MATOP_MULT_TRANSPOSE: 1149 if (shell->ops->multtranspose) 1150 *f = (void (*)(void))shell->ops->multtranspose; 1151 else 1152 *f = (((void (**)(void))mat->ops)[op]); 1153 break; 1154 default: 1155 *f = (((void (**)(void))mat->ops)[op]); 1156 } 1157 PetscFunctionReturn(0); 1158 } 1159