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