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