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