1 #define PETSCMAT_DLL 2 3 /* 4 This provides a simple shell for Fortran (and C programmers) to 5 create a very simple matrix class for use with KSP without coding 6 much of anything. 7 */ 8 9 #include "include/private/matimpl.h" /*I "petscmat.h" I*/ 10 #include "private/vecimpl.h" 11 12 typedef struct { 13 PetscErrorCode (*destroy)(Mat); 14 PetscErrorCode (*mult)(Mat,Vec,Vec); 15 PetscErrorCode (*multtranspose)(Mat,Vec,Vec); 16 PetscErrorCode (*getdiagonal)(Mat,Vec); 17 PetscTruth scale,shift; 18 PetscScalar vscale,vshift; 19 void *ctx; 20 } Mat_Shell; 21 22 #undef __FUNCT__ 23 #define __FUNCT__ "MatShellGetContext" 24 /*@ 25 MatShellGetContext - Returns the user-provided context associated with a shell matrix. 26 27 Not Collective 28 29 Input Parameter: 30 . mat - the matrix, should have been created with MatCreateShell() 31 32 Output Parameter: 33 . ctx - the user provided context 34 35 Level: advanced 36 37 Notes: 38 This routine is intended for use within various shell matrix routines, 39 as set with MatShellSetOperation(). 40 41 .keywords: matrix, shell, get, context 42 43 .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext() 44 @*/ 45 PetscErrorCode PETSCMAT_DLLEXPORT MatShellGetContext(Mat mat,void **ctx) 46 { 47 PetscErrorCode ierr; 48 PetscTruth flg; 49 50 PetscFunctionBegin; 51 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 52 PetscValidPointer(ctx,2); 53 ierr = PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 54 if (!flg) *ctx = 0; 55 else *ctx = ((Mat_Shell*)(mat->data))->ctx; 56 PetscFunctionReturn(0); 57 } 58 59 #undef __FUNCT__ 60 #define __FUNCT__ "MatDestroy_Shell" 61 PetscErrorCode MatDestroy_Shell(Mat mat) 62 { 63 PetscErrorCode ierr; 64 Mat_Shell *shell; 65 66 PetscFunctionBegin; 67 shell = (Mat_Shell*)mat->data; 68 if (shell->destroy) {ierr = (*shell->destroy)(mat);CHKERRQ(ierr);} 69 ierr = PetscFree(shell);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 #undef __FUNCT__ 74 #define __FUNCT__ "MatMult_Shell" 75 PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y) 76 { 77 Mat_Shell *shell = (Mat_Shell*)A->data; 78 PetscErrorCode ierr; 79 80 PetscFunctionBegin; 81 ierr = (*shell->mult)(A,x,y);CHKERRQ(ierr); 82 if (shell->shift && shell->scale) { 83 ierr = VecAXPBY(y,shell->vshift,shell->vscale,x);CHKERRQ(ierr); 84 } else if (shell->scale) { 85 ierr = VecScale(y,shell->vscale);CHKERRQ(ierr); 86 } else { 87 ierr = VecAXPY(y,shell->vshift,x);CHKERRQ(ierr); 88 } 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "MatMultTranspose_Shell" 94 PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y) 95 { 96 Mat_Shell *shell = (Mat_Shell*)A->data; 97 PetscErrorCode ierr; 98 99 PetscFunctionBegin; 100 ierr = (*shell->multtranspose)(A,x,y);CHKERRQ(ierr); 101 if (shell->shift && shell->scale) { 102 ierr = VecAXPBY(y,shell->vshift,shell->vscale,x);CHKERRQ(ierr); 103 } else if (shell->scale) { 104 ierr = VecScale(y,shell->vscale);CHKERRQ(ierr); 105 } else { 106 ierr = VecAXPY(y,shell->vshift,x);CHKERRQ(ierr); 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "MatGetDiagonal_Shell" 113 PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v) 114 { 115 Mat_Shell *shell = (Mat_Shell*)A->data; 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 ierr = (*shell->getdiagonal)(A,v);CHKERRQ(ierr); 120 if (shell->scale) { 121 ierr = VecScale(v,shell->vscale);CHKERRQ(ierr); 122 } 123 if (shell->shift) { 124 ierr = VecShift(v,shell->vshift);CHKERRQ(ierr); 125 } 126 PetscFunctionReturn(0); 127 } 128 129 #undef __FUNCT__ 130 #define __FUNCT__ "MatShift_Shell" 131 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a) 132 { 133 Mat_Shell *shell = (Mat_Shell*)Y->data; 134 135 PetscFunctionBegin; 136 if (shell->scale || shell->shift) { 137 shell->vshift += a; 138 } else { 139 shell->mult = Y->ops->mult; 140 Y->ops->mult = MatMult_Shell; 141 if (Y->ops->multtranspose) { 142 shell->multtranspose = Y->ops->multtranspose; 143 Y->ops->multtranspose = MatMultTranspose_Shell; 144 } 145 if (Y->ops->getdiagonal) { 146 shell->getdiagonal = Y->ops->getdiagonal; 147 Y->ops->getdiagonal = MatGetDiagonal_Shell; 148 } 149 shell->vshift = a; 150 } 151 shell->shift = PETSC_TRUE; 152 PetscFunctionReturn(0); 153 } 154 155 #undef __FUNCT__ 156 #define __FUNCT__ "MatScale_Shell" 157 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a) 158 { 159 Mat_Shell *shell = (Mat_Shell*)Y->data; 160 161 PetscFunctionBegin; 162 if (shell->scale || shell->shift) { 163 shell->vscale *= a; 164 } else { 165 shell->mult = Y->ops->mult; 166 Y->ops->mult = MatMult_Shell; 167 if (Y->ops->multtranspose) { 168 shell->multtranspose = Y->ops->multtranspose; 169 Y->ops->multtranspose = MatMultTranspose_Shell; 170 } 171 if (Y->ops->getdiagonal) { 172 shell->getdiagonal = Y->ops->getdiagonal; 173 Y->ops->getdiagonal = MatGetDiagonal_Shell; 174 } 175 shell->vscale = a; 176 } 177 shell->scale = PETSC_TRUE; 178 PetscFunctionReturn(0); 179 } 180 181 #undef __FUNCT__ 182 #define __FUNCT__ "MatAssemblyEnd_Shell" 183 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t) 184 { 185 Mat_Shell *shell = (Mat_Shell*)Y->data; 186 187 PetscFunctionBegin; 188 if ((shell->shift || shell->scale) && t == MAT_FINAL_ASSEMBLY) { 189 shell->scale = PETSC_FALSE; 190 shell->shift = PETSC_FALSE; 191 shell->vshift = 0.0; 192 shell->vscale = 1.0; 193 Y->ops->mult = shell->mult; 194 Y->ops->multtranspose = shell->multtranspose; 195 Y->ops->getdiagonal = shell->getdiagonal; 196 } 197 PetscFunctionReturn(0); 198 } 199 200 EXTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*); 201 202 #undef __FUNCT__ 203 #define __FUNCT__ "MatSetBlockSize_Shell" 204 PetscErrorCode MatSetBlockSize_Shell(Mat A,PetscInt bs) 205 { 206 PetscFunctionBegin; 207 PetscFunctionReturn(0); 208 } 209 210 static struct _MatOps MatOps_Values = {0, 211 0, 212 0, 213 0, 214 /* 4*/ 0, 215 0, 216 0, 217 0, 218 0, 219 0, 220 /*10*/ 0, 221 0, 222 0, 223 0, 224 0, 225 /*15*/ 0, 226 0, 227 0, 228 0, 229 0, 230 /*20*/ 0, 231 MatAssemblyEnd_Shell, 232 0, 233 0, 234 0, 235 /*25*/ 0, 236 0, 237 0, 238 0, 239 0, 240 /*30*/ 0, 241 0, 242 0, 243 0, 244 0, 245 /*35*/ 0, 246 0, 247 0, 248 0, 249 0, 250 /*40*/ 0, 251 0, 252 0, 253 0, 254 0, 255 /*45*/ 0, 256 MatScale_Shell, 257 MatShift_Shell, 258 0, 259 0, 260 /*50*/ MatSetBlockSize_Shell, 261 0, 262 0, 263 0, 264 0, 265 /*55*/ 0, 266 0, 267 0, 268 0, 269 0, 270 /*60*/ 0, 271 MatDestroy_Shell, 272 0, 273 0, 274 0, 275 /*65*/ 0, 276 0, 277 0, 278 0, 279 0, 280 /*70*/ 0, 281 MatConvert_Shell, 282 0, 283 0, 284 0, 285 /*75*/ 0, 286 0, 287 0, 288 0, 289 0, 290 /*80*/ 0, 291 0, 292 0, 293 0, 294 0, 295 /*85*/ 0, 296 0, 297 0, 298 0, 299 0, 300 /*90*/ 0, 301 0, 302 0, 303 0, 304 0, 305 /*95*/ 0, 306 0, 307 0, 308 0}; 309 310 /*MC 311 MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free. 312 313 Level: advanced 314 315 .seealso: MatCreateShell 316 M*/ 317 318 EXTERN_C_BEGIN 319 #undef __FUNCT__ 320 #define __FUNCT__ "MatCreate_Shell" 321 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_Shell(Mat A) 322 { 323 Mat_Shell *b; 324 PetscErrorCode ierr; 325 326 PetscFunctionBegin; 327 ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 328 329 ierr = PetscNewLog(A,Mat_Shell,&b);CHKERRQ(ierr); 330 A->data = (void*)b; 331 332 if (A->rmap.n == PETSC_DECIDE || A->cmap.n == PETSC_DECIDE) { 333 SETERRQ(PETSC_ERR_ARG_WRONG,"Must give local row and column count for matrix"); 334 } 335 336 ierr = PetscMapSetBlockSize(&A->rmap,1);CHKERRQ(ierr); 337 ierr = PetscMapSetBlockSize(&A->cmap,1);CHKERRQ(ierr); 338 ierr = PetscMapSetUp(&A->rmap);CHKERRQ(ierr); 339 ierr = PetscMapSetUp(&A->cmap);CHKERRQ(ierr); 340 341 b->ctx = 0; 342 b->scale = PETSC_FALSE; 343 b->shift = PETSC_FALSE; 344 b->vshift = 0.0; 345 b->vscale = 1.0; 346 b->mult = 0; 347 b->multtranspose = 0; 348 b->getdiagonal = 0; 349 A->assembled = PETSC_TRUE; 350 A->preallocated = PETSC_FALSE; 351 ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr); 352 PetscFunctionReturn(0); 353 } 354 EXTERN_C_END 355 356 #undef __FUNCT__ 357 #define __FUNCT__ "MatCreateShell" 358 /*@C 359 MatCreateShell - Creates a new matrix class for use with a user-defined 360 private data storage format. 361 362 Collective on MPI_Comm 363 364 Input Parameters: 365 + comm - MPI communicator 366 . m - number of local rows (must be given) 367 . n - number of local columns (must be given) 368 . M - number of global rows (may be PETSC_DETERMINE) 369 . N - number of global columns (may be PETSC_DETERMINE) 370 - ctx - pointer to data needed by the shell matrix routines 371 372 Output Parameter: 373 . A - the matrix 374 375 Level: advanced 376 377 Usage: 378 $ extern int mult(Mat,Vec,Vec); 379 $ MatCreateShell(comm,m,n,M,N,ctx,&mat); 380 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 381 $ [ Use matrix for operations that have been set ] 382 $ MatDestroy(mat); 383 384 Notes: 385 The shell matrix type is intended to provide a simple class to use 386 with KSP (such as, for use with matrix-free methods). You should not 387 use the shell type if you plan to define a complete matrix class. 388 389 Fortran Notes: The context can only be an integer or a PetscObject 390 unfortunately it cannot be a Fortran array or derived type. 391 392 PETSc requires that matrices and vectors being used for certain 393 operations are partitioned accordingly. For example, when 394 creating a shell matrix, A, that supports parallel matrix-vector 395 products using MatMult(A,x,y) the user should set the number 396 of local matrix rows to be the number of local elements of the 397 corresponding result vector, y. Note that this is information is 398 required for use of the matrix interface routines, even though 399 the shell matrix may not actually be physically partitioned. 400 For example, 401 402 $ 403 $ Vec x, y 404 $ extern int mult(Mat,Vec,Vec); 405 $ Mat A 406 $ 407 $ VecCreateMPI(comm,PETSC_DECIDE,M,&y); 408 $ VecCreateMPI(comm,PETSC_DECIDE,N,&x); 409 $ VecGetLocalSize(y,&m); 410 $ VecGetLocalSize(x,&n); 411 $ MatCreateShell(comm,m,n,M,N,ctx,&A); 412 $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult); 413 $ MatMult(A,x,y); 414 $ MatDestroy(A); 415 $ VecDestroy(y); VecDestroy(x); 416 $ 417 418 .keywords: matrix, shell, create 419 420 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext() 421 @*/ 422 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A) 423 { 424 PetscErrorCode ierr; 425 426 PetscFunctionBegin; 427 ierr = MatCreate(comm,A);CHKERRQ(ierr); 428 ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); 429 430 ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr); 431 ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr); 432 PetscFunctionReturn(0); 433 } 434 435 #undef __FUNCT__ 436 #define __FUNCT__ "MatShellSetContext" 437 /*@ 438 MatShellSetContext - sets the context for a shell matrix 439 440 Collective on Mat 441 442 Input Parameters: 443 + mat - the shell matrix 444 - ctx - the context 445 446 Level: advanced 447 448 Fortran Notes: The context can only be an integer or a PetscObject 449 unfortunately it cannot be a Fortran array or derived type. 450 451 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation() 452 @*/ 453 PetscErrorCode PETSCMAT_DLLEXPORT MatShellSetContext(Mat mat,void *ctx) 454 { 455 Mat_Shell *shell = (Mat_Shell*)mat->data; 456 PetscErrorCode ierr; 457 PetscTruth flg; 458 459 PetscFunctionBegin; 460 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 461 ierr = PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 462 if (flg) { 463 shell->ctx = ctx; 464 } 465 PetscFunctionReturn(0); 466 } 467 468 #undef __FUNCT__ 469 #define __FUNCT__ "MatShellSetOperation" 470 /*@C 471 MatShellSetOperation - Allows user to set a matrix operation for 472 a shell matrix. 473 474 Collective on Mat 475 476 Input Parameters: 477 + mat - the shell matrix 478 . op - the name of the operation 479 - f - the function that provides the operation. 480 481 Level: advanced 482 483 Usage: 484 $ extern PetscErrorCode usermult(Mat,Vec,Vec); 485 $ ierr = MatCreateShell(comm,m,n,M,N,ctx,&A); 486 $ ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult); 487 488 Notes: 489 See the file include/petscmat.h for a complete list of matrix 490 operations, which all have the form MATOP_<OPERATION>, where 491 <OPERATION> is the name (in all capital letters) of the 492 user interface routine (e.g., MatMult() -> MATOP_MULT). 493 494 All user-provided functions should have the same calling 495 sequence as the usual matrix interface routines, since they 496 are intended to be accessed via the usual matrix interface 497 routines, e.g., 498 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 499 500 Within each user-defined routine, the user should call 501 MatShellGetContext() to obtain the user-defined context that was 502 set by MatCreateShell(). 503 504 .keywords: matrix, shell, set, operation 505 506 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext() 507 @*/ 508 PetscErrorCode PETSCMAT_DLLEXPORT MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void)) 509 { 510 PetscErrorCode ierr; 511 PetscTruth flg; 512 513 PetscFunctionBegin; 514 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 515 if (op == MATOP_DESTROY) { 516 ierr = PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 517 if (flg) { 518 Mat_Shell *shell = (Mat_Shell*)mat->data; 519 shell->destroy = (PetscErrorCode (*)(Mat)) f; 520 } else mat->ops->destroy = (PetscErrorCode (*)(Mat)) f; 521 } 522 else if (op == MATOP_VIEW) mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer)) f; 523 else (((void(**)(void))mat->ops)[op]) = f; 524 525 PetscFunctionReturn(0); 526 } 527 528 #undef __FUNCT__ 529 #define __FUNCT__ "MatShellGetOperation" 530 /*@C 531 MatShellGetOperation - Gets a matrix function for a shell matrix. 532 533 Not Collective 534 535 Input Parameters: 536 + mat - the shell matrix 537 - op - the name of the operation 538 539 Output Parameter: 540 . f - the function that provides the operation. 541 542 Level: advanced 543 544 Notes: 545 See the file include/petscmat.h for a complete list of matrix 546 operations, which all have the form MATOP_<OPERATION>, where 547 <OPERATION> is the name (in all capital letters) of the 548 user interface routine (e.g., MatMult() -> MATOP_MULT). 549 550 All user-provided functions have the same calling 551 sequence as the usual matrix interface routines, since they 552 are intended to be accessed via the usual matrix interface 553 routines, e.g., 554 $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec) 555 556 Within each user-defined routine, the user should call 557 MatShellGetContext() to obtain the user-defined context that was 558 set by MatCreateShell(). 559 560 .keywords: matrix, shell, set, operation 561 562 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext() 563 @*/ 564 PetscErrorCode PETSCMAT_DLLEXPORT MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void)) 565 { 566 PetscErrorCode ierr; 567 PetscTruth flg; 568 569 PetscFunctionBegin; 570 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 571 if (op == MATOP_DESTROY) { 572 ierr = PetscTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr); 573 if (flg) { 574 Mat_Shell *shell = (Mat_Shell*)mat->data; 575 *f = (void(*)(void))shell->destroy; 576 } else { 577 *f = (void(*)(void))mat->ops->destroy; 578 } 579 } else if (op == MATOP_VIEW) { 580 *f = (void(*)(void))mat->ops->view; 581 } else { 582 *f = (((void(**)(void))mat->ops)[op]); 583 } 584 585 PetscFunctionReturn(0); 586 } 587 588