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