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