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 .keywords: matrix, shell, create 763 764 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts() 765 @*/ 766 PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 767 { 768 PetscErrorCode ierr; 769 770 PetscFunctionBegin; 771 ierr = MatCreate(comm,A);CHKERRQ(ierr); 772 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 773 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 774 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 775 ierr = MatSetUp(*A);CHKERRQ(ierr); 776 PetscFunctionReturn(0); 777 } 778 779 /*@ 780 MatShellSetContext - sets the context for a shell matrix 781 782 Logically Collective on Mat 783 784 Input Parameters: 785 + mat - the shell matrix 786 - ctx - the context 787 788 Level: advanced 789 790 Fortran Notes: 791 To use this from Fortran you must write a Fortran interface definition for this 792 function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument. 793 794 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 795 @*/ 796 PetscErrorCode MatShellSetContext(Mat mat,void *ctx) 797 { 798 Mat_Shell *shell = (Mat_Shell*)mat->data; 799 PetscErrorCode ierr; 800 PetscBool flg; 801 802 PetscFunctionBegin; 803 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 804 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 805 if (flg) { 806 shell->ctx = ctx; 807 } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix"); 808 PetscFunctionReturn(0); 809 } 810 811 /*@ 812 MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately 813 after MatCreateShell() 814 815 Logically Collective on Mat 816 817 Input Parameter: 818 . mat - the shell matrix 819 820 Level: advanced 821 822 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation() 823 @*/ 824 PetscErrorCode MatShellSetManageScalingShifts(Mat A) 825 { 826 PetscErrorCode ierr; 827 Mat_Shell *shell; 828 PetscBool flg; 829 830 PetscFunctionBegin; 831 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 832 ierr = PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);CHKERRQ(ierr); 833 if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 834 shell = (Mat_Shell*)A->data; 835 shell->managescalingshifts = PETSC_FALSE; 836 A->ops->diagonalset = NULL; 837 A->ops->diagonalscale = NULL; 838 A->ops->scale = NULL; 839 A->ops->shift = NULL; 840 A->ops->axpy = NULL; 841 PetscFunctionReturn(0); 842 } 843 844 /*@C 845 MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function. 846 847 Logically Collective on Mat 848 849 Input Parameters: 850 + mat - the shell matrix 851 . f - the function 852 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 853 - ctx - an optional context for the function 854 855 Output Parameter: 856 . flg - PETSC_TRUE if the multiply is likely correct 857 858 Options Database: 859 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 860 861 Level: advanced 862 863 Fortran Notes: 864 Not supported from Fortran 865 866 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose() 867 @*/ 868 PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 869 { 870 PetscErrorCode ierr; 871 PetscInt m,n; 872 Mat mf,Dmf,Dmat,Ddiff; 873 PetscReal Diffnorm,Dmfnorm; 874 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 875 876 PetscFunctionBegin; 877 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 878 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr); 879 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 880 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr); 881 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 882 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 883 884 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 885 ierr = MatComputeExplicitOperator(mat,&Dmat);CHKERRQ(ierr); 886 887 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 888 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 889 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 890 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 891 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 892 flag = PETSC_FALSE; 893 if (v) { 894 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); 895 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 896 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 897 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr); 898 } 899 } else if (v) { 900 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 901 } 902 if (flg) *flg = flag; 903 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 904 ierr = MatDestroy(&mf);CHKERRQ(ierr); 905 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 906 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 907 PetscFunctionReturn(0); 908 } 909 910 /*@C 911 MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function. 912 913 Logically Collective on Mat 914 915 Input Parameters: 916 + mat - the shell matrix 917 . f - the function 918 . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated 919 - ctx - an optional context for the function 920 921 Output Parameter: 922 . flg - PETSC_TRUE if the multiply is likely correct 923 924 Options Database: 925 . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference 926 927 Level: advanced 928 929 Fortran Notes: 930 Not supported from Fortran 931 932 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult() 933 @*/ 934 PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg) 935 { 936 PetscErrorCode ierr; 937 Vec x,y,z; 938 PetscInt m,n,M,N; 939 Mat mf,Dmf,Dmat,Ddiff; 940 PetscReal Diffnorm,Dmfnorm; 941 PetscBool v = PETSC_FALSE, flag = PETSC_TRUE; 942 943 PetscFunctionBegin; 944 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 945 ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr); 946 ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr); 947 ierr = VecDuplicate(y,&z);CHKERRQ(ierr); 948 ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr); 949 ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr); 950 ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr); 951 ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr); 952 ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr); 953 ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr); 954 ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr); 955 ierr = MatComputeExplicitOperatorTranspose(mat,&Dmat);CHKERRQ(ierr); 956 957 ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr); 958 ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 959 ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr); 960 ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr); 961 if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) { 962 flag = PETSC_FALSE; 963 if (v) { 964 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); 965 ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 966 ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 967 ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr); 968 } 969 } else if (v) { 970 ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");CHKERRQ(ierr); 971 } 972 if (flg) *flg = flag; 973 ierr = MatDestroy(&mf);CHKERRQ(ierr); 974 ierr = MatDestroy(&Dmat);CHKERRQ(ierr); 975 ierr = MatDestroy(&Ddiff);CHKERRQ(ierr); 976 ierr = MatDestroy(&Dmf);CHKERRQ(ierr); 977 ierr = VecDestroy(&x);CHKERRQ(ierr); 978 ierr = VecDestroy(&y);CHKERRQ(ierr); 979 ierr = VecDestroy(&z);CHKERRQ(ierr); 980 PetscFunctionReturn(0); 981 } 982 983 /*@C 984 MatShellSetOperation - Allows user to set a matrix operation for a shell matrix. 985 986 Logically Collective on Mat 987 988 Input Parameters: 989 + mat - the shell matrix 990 . op - the name of the operation 991 - f - the function that provides the operation. 992 993 Level: advanced 994 995 Usage: 996 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 997 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 998 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 999 1000 Notes: 1001 See the file include/petscmat.h for a complete list of matrix 1002 operations, which all have the form MATOP_<OPERATION>, where 1003 <OPERATION> is the name (in all capital letters) of the 1004 user interface routine (e.g., MatMult() -> MATOP_MULT). 1005 1006 All user-provided functions (except for MATOP_DESTROY) should have the same calling 1007 sequence as the usual matrix interface routines, since they 1008 are intended to be accessed via the usual matrix interface 1009 routines, e.g., 1010 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1011 1012 In particular each function MUST return an error code of 0 on success and 1013 nonzero on failure. 1014 1015 Within each user-defined routine, the user should call 1016 MatShellGetContext() to obtain the user-defined context that was 1017 set by MatCreateShell(). 1018 1019 Fortran Notes: 1020 For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not 1021 generate a matrix. See src/mat/examples/tests/ex120f.F 1022 1023 Use MatSetOperation() to set an operation for any matrix type 1024 1025 .keywords: matrix, shell, set, operation 1026 1027 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts() 1028 @*/ 1029 PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 1030 { 1031 PetscBool flg; 1032 Mat_Shell *shell; 1033 PetscErrorCode ierr; 1034 1035 PetscFunctionBegin; 1036 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1037 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1038 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1039 shell = (Mat_Shell*)mat->data; 1040 1041 switch (op) { 1042 case MATOP_DESTROY: 1043 shell->ops->destroy = (PetscErrorCode (*)(Mat))f; 1044 break; 1045 case MATOP_VIEW: 1046 if (!mat->ops->viewnative) { 1047 mat->ops->viewnative = mat->ops->view; 1048 } 1049 mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f; 1050 break; 1051 case MATOP_COPY: 1052 shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f; 1053 break; 1054 case MATOP_DIAGONAL_SET: 1055 case MATOP_DIAGONAL_SCALE: 1056 case MATOP_SHIFT: 1057 case MATOP_SCALE: 1058 case MATOP_AXPY: 1059 if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()"); 1060 (((void(**)(void))mat->ops)[op]) = f; 1061 break; 1062 case MATOP_GET_DIAGONAL: 1063 if (shell->managescalingshifts) { 1064 shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1065 mat->ops->getdiagonal = MatGetDiagonal_Shell; 1066 } else { 1067 shell->ops->getdiagonal = NULL; 1068 mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f; 1069 } 1070 break; 1071 case MATOP_MULT: 1072 if (shell->managescalingshifts) { 1073 shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1074 mat->ops->mult = MatMult_Shell; 1075 } else { 1076 shell->ops->mult = NULL; 1077 mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1078 } 1079 break; 1080 case MATOP_MULT_TRANSPOSE: 1081 if (shell->managescalingshifts) { 1082 shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1083 mat->ops->multtranspose = MatMultTranspose_Shell; 1084 } else { 1085 shell->ops->multtranspose = NULL; 1086 mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f; 1087 } 1088 break; 1089 default: 1090 (((void(**)(void))mat->ops)[op]) = f; 1091 break; 1092 } 1093 PetscFunctionReturn(0); 1094 } 1095 1096 /*@C 1097 MatShellGetOperation - Gets a matrix function for a shell matrix. 1098 1099 Not Collective 1100 1101 Input Parameters: 1102 + mat - the shell matrix 1103 - op - the name of the operation 1104 1105 Output Parameter: 1106 . f - the function that provides the operation. 1107 1108 Level: advanced 1109 1110 Notes: 1111 See the file include/petscmat.h for a complete list of matrix 1112 operations, which all have the form MATOP_<OPERATION>, where 1113 <OPERATION> is the name (in all capital letters) of the 1114 user interface routine (e.g., MatMult() -> MATOP_MULT). 1115 1116 All user-provided functions have the same calling 1117 sequence as the usual matrix interface routines, since they 1118 are intended to be accessed via the usual matrix interface 1119 routines, e.g., 1120 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 1121 1122 Within each user-defined routine, the user should call 1123 MatShellGetContext() to obtain the user-defined context that was 1124 set by MatCreateShell(). 1125 1126 .keywords: matrix, shell, set, operation 1127 1128 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 1129 @*/ 1130 PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 1131 { 1132 PetscBool flg; 1133 Mat_Shell *shell; 1134 PetscErrorCode ierr; 1135 1136 PetscFunctionBegin; 1137 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 1138 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 1139 if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices"); 1140 shell = (Mat_Shell*)mat->data; 1141 1142 switch (op) { 1143 case MATOP_DESTROY: 1144 *f = (void (*)(void))shell->ops->destroy; 1145 break; 1146 case MATOP_VIEW: 1147 *f = (void (*)(void))mat->ops->view; 1148 break; 1149 case MATOP_COPY: 1150 *f = (void (*)(void))shell->ops->copy; 1151 break; 1152 case MATOP_DIAGONAL_SET: 1153 case MATOP_DIAGONAL_SCALE: 1154 case MATOP_SHIFT: 1155 case MATOP_SCALE: 1156 case MATOP_AXPY: 1157 *f = (((void (**)(void))mat->ops)[op]); 1158 break; 1159 case MATOP_GET_DIAGONAL: 1160 if (shell->ops->getdiagonal) 1161 *f = (void (*)(void))shell->ops->getdiagonal; 1162 else 1163 *f = (((void (**)(void))mat->ops)[op]); 1164 break; 1165 case MATOP_MULT: 1166 if (shell->ops->mult) 1167 *f = (void (*)(void))shell->ops->mult; 1168 else 1169 *f = (((void (**)(void))mat->ops)[op]); 1170 break; 1171 case MATOP_MULT_TRANSPOSE: 1172 if (shell->ops->multtranspose) 1173 *f = (void (*)(void))shell->ops->multtranspose; 1174 else 1175 *f = (((void (**)(void))mat->ops)[op]); 1176 break; 1177 default: 1178 *f = (((void (**)(void))mat->ops)[op]); 1179 } 1180 PetscFunctionReturn(0); 1181 } 1182