1 2 /* 3 The basic KSP routines, Create, View etc. are here. 4 */ 5 #include <petsc-private/kspimpl.h> /*I "petscksp.h" I*/ 6 7 /* Logging support */ 8 PetscClassId KSP_CLASSID; 9 PetscClassId DMKSP_CLASSID; 10 PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve; 11 12 /* 13 Contains the list of registered KSP routines 14 */ 15 PetscFList KSPList = 0; 16 PetscBool KSPRegisterAllCalled = PETSC_FALSE; 17 18 #undef __FUNCT__ 19 #define __FUNCT__ "KSPLoad" 20 /*@C 21 KSPLoad - Loads a KSP that has been stored in binary with KSPView(). 22 23 Collective on PetscViewer 24 25 Input Parameters: 26 + newdm - the newly loaded KSP, this needs to have been created with KSPCreate() or 27 some related function before a call to KSPLoad(). 28 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() 29 30 Level: intermediate 31 32 Notes: 33 The type is determined by the data in the file, any type set into the KSP before this call is ignored. 34 35 Notes for advanced users: 36 Most users should not need to know the details of the binary storage 37 format, since KSPLoad() and KSPView() completely hide these details. 38 But for anyone who's interested, the standard binary matrix storage 39 format is 40 .vb 41 has not yet been determined 42 .ve 43 44 .seealso: PetscViewerBinaryOpen(), KSPView(), MatLoad(), VecLoad() 45 @*/ 46 PetscErrorCode KSPLoad(KSP newdm, PetscViewer viewer) 47 { 48 PetscErrorCode ierr; 49 PetscBool isbinary; 50 PetscInt classid; 51 char type[256]; 52 PC pc; 53 54 PetscFunctionBegin; 55 PetscValidHeaderSpecific(newdm,KSP_CLASSID,1); 56 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 57 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 58 if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()"); 59 60 ierr = PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);CHKERRQ(ierr); 61 if (classid != KSP_FILE_CLASSID) SETERRQ(((PetscObject)newdm)->comm,PETSC_ERR_ARG_WRONG,"Not KSP next in file"); 62 ierr = PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);CHKERRQ(ierr); 63 ierr = KSPSetType(newdm, type);CHKERRQ(ierr); 64 if (newdm->ops->load) { 65 ierr = (*newdm->ops->load)(newdm,viewer);CHKERRQ(ierr); 66 } 67 ierr = KSPGetPC(newdm,&pc);CHKERRQ(ierr); 68 ierr = PCLoad(pc,viewer);CHKERRQ(ierr); 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "KSPView" 74 /*@C 75 KSPView - Prints the KSP data structure. 76 77 Collective on KSP 78 79 Input Parameters: 80 + ksp - the Krylov space context 81 - viewer - visualization context 82 83 Options Database Keys: 84 . -ksp_view - print the ksp data structure at the end of a KSPSolve call 85 86 Note: 87 The available visualization contexts include 88 + PETSC_VIEWER_STDOUT_SELF - standard output (default) 89 - PETSC_VIEWER_STDOUT_WORLD - synchronized standard 90 output where only the first processor opens 91 the file. All other processors send their 92 data to the first processor to print. 93 94 The user can open an alternative visualization context with 95 PetscViewerASCIIOpen() - output to a specified file. 96 97 Level: beginner 98 99 .keywords: KSP, view 100 101 .seealso: PCView(), PetscViewerASCIIOpen() 102 @*/ 103 PetscErrorCode KSPView(KSP ksp,PetscViewer viewer) 104 { 105 PetscErrorCode ierr; 106 PetscBool iascii,isbinary; 107 108 PetscFunctionBegin; 109 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 110 if (!viewer) viewer = PETSC_VIEWER_STDOUT_(((PetscObject)ksp)->comm); 111 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 112 PetscCheckSameComm(ksp,1,viewer,2); 113 114 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 115 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 116 if (iascii) { 117 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)ksp,viewer,"KSP Object");CHKERRQ(ierr); 118 if (ksp->ops->view) { 119 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 120 ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr); 121 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 122 } 123 if (ksp->guess_zero) {ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D, initial guess is zero\n",ksp->max_it);CHKERRQ(ierr);} 124 else {ierr = PetscViewerASCIIPrintf(viewer," maximum iterations=%D\n", ksp->max_it);CHKERRQ(ierr);} 125 if (ksp->guess_knoll) {ierr = PetscViewerASCIIPrintf(viewer," using preconditioner applied to right hand side for initial guess\n");CHKERRQ(ierr);} 126 ierr = PetscViewerASCIIPrintf(viewer," tolerances: relative=%G, absolute=%G, divergence=%G\n",ksp->rtol,ksp->abstol,ksp->divtol);CHKERRQ(ierr); 127 if (ksp->pc_side == PC_RIGHT) {ierr = PetscViewerASCIIPrintf(viewer," right preconditioning\n");CHKERRQ(ierr);} 128 else if (ksp->pc_side == PC_SYMMETRIC) {ierr = PetscViewerASCIIPrintf(viewer," symmetric preconditioning\n");CHKERRQ(ierr);} 129 else {ierr = PetscViewerASCIIPrintf(viewer," left preconditioning\n");CHKERRQ(ierr);} 130 if (ksp->guess) {ierr = PetscViewerASCIIPrintf(viewer," using Fischers initial guess method %D with size %D\n",ksp->guess->method,ksp->guess->maxl);CHKERRQ(ierr);} 131 if (ksp->dscale) {ierr = PetscViewerASCIIPrintf(viewer," diagonally scaled system\n");CHKERRQ(ierr);} 132 if (ksp->nullsp) {ierr = PetscViewerASCIIPrintf(viewer," has attached null space\n");CHKERRQ(ierr);} 133 if (!ksp->guess_zero) {ierr = PetscViewerASCIIPrintf(viewer," using nonzero initial guess\n");CHKERRQ(ierr);} 134 ierr = PetscViewerASCIIPrintf(viewer," using %s norm type for convergence test\n",KSPNormTypes[ksp->normtype]);CHKERRQ(ierr); 135 } else if (isbinary) { 136 PetscInt classid = KSP_FILE_CLASSID; 137 MPI_Comm comm; 138 PetscMPIInt rank; 139 char type[256]; 140 141 ierr = PetscObjectGetComm((PetscObject)ksp,&comm);CHKERRQ(ierr); 142 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 143 if (!rank) { 144 ierr = PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 145 ierr = PetscStrncpy(type,((PetscObject)ksp)->type_name,256);CHKERRQ(ierr); 146 ierr = PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);CHKERRQ(ierr); 147 } 148 if (ksp->ops->view) { 149 ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr); 150 } 151 } else { 152 if (ksp->ops->view) { 153 ierr = (*ksp->ops->view)(ksp,viewer);CHKERRQ(ierr); 154 } 155 } 156 if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);} 157 ierr = PCView(ksp->pc,viewer);CHKERRQ(ierr); 158 PetscFunctionReturn(0); 159 } 160 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "KSPSetNormType" 164 /*@ 165 KSPSetNormType - Sets the norm that is used for convergence testing. 166 167 Logically Collective on KSP 168 169 Input Parameter: 170 + ksp - Krylov solver context 171 - normtype - one of 172 $ KSP_NORM_NONE - skips computing the norm, this should only be used if you are using 173 $ the Krylov method as a smoother with a fixed small number of iterations. 174 $ Implicitly sets KSPSkipConverged as KSP convergence test. 175 $ Supported only by CG, Richardson, Bi-CG-stab, CR, and CGS methods. 176 $ KSP_NORM_PRECONDITIONED - the default for left preconditioned solves, uses the l2 norm 177 $ of the preconditioned residual 178 $ KSP_NORM_UNPRECONDITIONED - uses the l2 norm of the true b - Ax residual, supported only by 179 $ CG, CHEBYSHEV, and RICHARDSON, automatically true for right (see KSPSetPCSide()) 180 $ preconditioning.. 181 $ KSP_NORM_NATURAL - supported by KSPCG, KSPCR, KSPCGNE, KSPCGS 182 183 184 Options Database Key: 185 . -ksp_norm_type <none,preconditioned,unpreconditioned,natural> 186 187 Notes: 188 Currently only works with the CG, Richardson, Bi-CG-stab, CR, and CGS methods. 189 190 Level: advanced 191 192 .keywords: KSP, create, context, norms 193 194 .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetCheckNormIteration() 195 @*/ 196 PetscErrorCode KSPSetNormType(KSP ksp,KSPNormType normtype) 197 { 198 PetscErrorCode ierr; 199 200 PetscFunctionBegin; 201 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 202 PetscValidLogicalCollectiveEnum(ksp,normtype,2); 203 ksp->normtype = normtype; 204 if (normtype == KSP_NORM_NONE) { 205 ierr = KSPSetConvergenceTest(ksp,KSPSkipConverged,0,0);CHKERRQ(ierr); 206 ierr = PetscInfo(ksp,"Warning: setting KSPNormType to skip computing the norm\n\ 207 KSP convergence test is implicitly set to KSPSkipConverged\n");CHKERRQ(ierr); 208 } 209 PetscFunctionReturn(0); 210 } 211 212 #undef __FUNCT__ 213 #define __FUNCT__ "KSPSetCheckNormIteration" 214 /*@ 215 KSPSetCheckNormIteration - Sets the first iteration at which the norm of the residual will be 216 computed and used in the convergence test. 217 218 Logically Collective on KSP 219 220 Input Parameter: 221 + ksp - Krylov solver context 222 - it - use -1 to check at all iterations 223 224 Notes: 225 Currently only works with KSPCG, KSPBCGS and KSPIBCGS 226 227 Use KSPSetNormType(ksp,KSP_NORM_NONE) to never check the norm 228 229 On steps where the norm is not computed, the previous norm is still in the variable, so if you run with, for example, 230 -ksp_monitor the residual norm will appear to be unchanged for several iterations (though it is not really unchanged). 231 Level: advanced 232 233 .keywords: KSP, create, context, norms 234 235 .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetNormType() 236 @*/ 237 PetscErrorCode KSPSetCheckNormIteration(KSP ksp,PetscInt it) 238 { 239 PetscFunctionBegin; 240 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 241 PetscValidLogicalCollectiveInt(ksp,it,2); 242 ksp->chknorm = it; 243 PetscFunctionReturn(0); 244 } 245 246 #undef __FUNCT__ 247 #define __FUNCT__ "KSPSetLagNorm" 248 /*@ 249 KSPSetLagNorm - Lags the residual norm calculation so that it is computed as part of the MPI_Allreduce() for 250 computing the inner products for the next iteration. This can reduce communication costs at the expense of doing 251 one additional iteration. 252 253 254 Logically Collective on KSP 255 256 Input Parameter: 257 + ksp - Krylov solver context 258 - flg - PETSC_TRUE or PETSC_FALSE 259 260 Options Database Keys: 261 . -ksp_lag_norm - lag the calculated residual norm 262 263 Notes: 264 Currently only works with KSPIBCGS. 265 266 Use KSPSetNormType(ksp,KSP_NORM_NONE) to never check the norm 267 268 If you lag the norm and run with, for example, -ksp_monitor, the residual norm reported will be the lagged one. 269 Level: advanced 270 271 .keywords: KSP, create, context, norms 272 273 .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetNormType(), KSPSetCheckNormIteration() 274 @*/ 275 PetscErrorCode KSPSetLagNorm(KSP ksp,PetscBool flg) 276 { 277 PetscFunctionBegin; 278 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 279 PetscValidLogicalCollectiveBool(ksp,flg,2); 280 ksp->lagnorm = flg; 281 PetscFunctionReturn(0); 282 } 283 284 #undef __FUNCT__ 285 #define __FUNCT__ "KSPSetSupportedNorm" 286 /*@ 287 KSPSetSupportedNorm - Sets a norm and preconditioner side supported by a KSP 288 289 Logically Collective 290 291 Input Arguments: 292 + ksp - Krylov method 293 . normtype - supported norm type 294 . pcside - preconditioner side that can be used with this norm 295 - preference - integer preference for this combination, larger values have higher priority 296 297 Level: developer 298 299 Notes: 300 This function should be called from the implementation files KSPCreate_XXX() to declare 301 which norms and preconditioner sides are supported. Users should not need to call this 302 function. 303 304 KSP_NORM_NONE is supported by default with all KSP methods and any PC side. If a KSP explicitly does not support 305 KSP_NORM_NONE, it should set this by setting priority=0. 306 307 .seealso: KSPSetNormType(), KSPSetPCSide() 308 @*/ 309 PetscErrorCode KSPSetSupportedNorm(KSP ksp,KSPNormType normtype,PCSide pcside,PetscInt priority) 310 { 311 312 PetscFunctionBegin; 313 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 314 ksp->normsupporttable[normtype][pcside] = priority; 315 PetscFunctionReturn(0); 316 } 317 318 #undef __FUNCT__ 319 #define __FUNCT__ "KSPNormSupportTableReset_Private" 320 PetscErrorCode KSPNormSupportTableReset_Private(KSP ksp) 321 { 322 PetscErrorCode ierr; 323 324 PetscFunctionBegin; 325 ierr = PetscMemzero(ksp->normsupporttable,sizeof(ksp->normsupporttable));CHKERRQ(ierr); 326 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);CHKERRQ(ierr); 327 ierr = KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330 331 #undef __FUNCT__ 332 #define __FUNCT__ "KSPSetUpNorms_Private" 333 PetscErrorCode KSPSetUpNorms_Private(KSP ksp,KSPNormType *normtype,PCSide *pcside) 334 { 335 PetscInt i,j,best,ibest = 0,jbest = 0; 336 337 PetscFunctionBegin; 338 best = 0; 339 for (i=0; i<KSP_NORM_MAX; i++) { 340 for (j=0; j<PC_SIDE_MAX; j++) { 341 if ((ksp->normtype == KSP_NORM_DEFAULT || ksp->normtype == i) 342 && (ksp->pc_side == PC_SIDE_DEFAULT || ksp->pc_side == j) 343 && (ksp->normsupporttable[i][j] > best)) { 344 if (ksp->normtype == KSP_NORM_DEFAULT && i == KSP_NORM_NONE && ksp->normsupporttable[i][j] <= 1) 345 continue; /* Skip because we don't want to default to no norms unless set by the KSP (preonly). */ 346 best = ksp->normsupporttable[i][j]; 347 ibest = i; 348 jbest = j; 349 } 350 } 351 } 352 if (best < 1) { 353 if (ksp->normtype == KSP_NORM_DEFAULT && ksp->pc_side == PC_SIDE_DEFAULT) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_PLIB,"The %s KSP implementation did not call KSPSetSupportedNorm()",((PetscObject)ksp)->type_name); 354 if (ksp->normtype == KSP_NORM_DEFAULT) SETERRQ2(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"KSP %s does not support %s",((PetscObject)ksp)->type_name,PCSides[ksp->pc_side]); 355 if (ksp->pc_side == PC_SIDE_DEFAULT) SETERRQ2(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"KSP %s does not support %s",((PetscObject)ksp)->type_name,KSPNormTypes[ksp->normtype]); 356 SETERRQ3(((PetscObject)ksp)->comm,PETSC_ERR_SUP,"KSP %s does not support %s with %s",((PetscObject)ksp)->type_name,KSPNormTypes[ksp->normtype],PCSides[ksp->pc_side]); 357 } 358 *normtype = (KSPNormType)ibest; 359 *pcside = (PCSide)jbest; 360 PetscFunctionReturn(0); 361 } 362 363 #undef __FUNCT__ 364 #define __FUNCT__ "KSPGetNormType" 365 /*@ 366 KSPGetNormType - Gets the norm that is used for convergence testing. 367 368 Not Collective 369 370 Input Parameter: 371 . ksp - Krylov solver context 372 373 Output Parameter: 374 . normtype - norm that is used for convergence testing 375 376 Level: advanced 377 378 .keywords: KSP, create, context, norms 379 380 .seealso: KSPNormType, KSPSetNormType(), KSPSkipConverged() 381 @*/ 382 PetscErrorCode KSPGetNormType(KSP ksp, KSPNormType *normtype) 383 { 384 PetscErrorCode ierr; 385 386 PetscFunctionBegin; 387 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 388 PetscValidPointer(normtype,2); 389 ierr = KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);CHKERRQ(ierr); 390 *normtype = ksp->normtype; 391 PetscFunctionReturn(0); 392 } 393 394 #if defined(PETSC_HAVE_AMS) 395 #undef __FUNCT__ 396 #define __FUNCT__ "KSPPublish_Petsc" 397 static PetscErrorCode KSPPublish_Petsc(PetscObject obj) 398 { 399 KSP ksp = (KSP) obj; 400 PetscErrorCode ierr; 401 402 PetscFunctionBegin; 403 ierr = AMS_Memory_add_field(obj->amem,"its",&ksp->its,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF);CHKERRQ(ierr); 404 PetscFunctionReturn(0); 405 } 406 #endif 407 408 #undef __FUNCT__ 409 #define __FUNCT__ "KSPSetOperators" 410 /*@ 411 KSPSetOperators - Sets the matrix associated with the linear system 412 and a (possibly) different one associated with the preconditioner. 413 414 Collective on KSP and Mat 415 416 Input Parameters: 417 + ksp - the KSP context 418 . Amat - the matrix associated with the linear system 419 . Pmat - the matrix to be used in constructing the preconditioner, usually the 420 same as Amat. 421 - flag - flag indicating information about the preconditioner matrix structure 422 during successive linear solves. This flag is ignored the first time a 423 linear system is solved, and thus is irrelevant when solving just one linear 424 system. 425 426 Notes: 427 The flag can be used to eliminate unnecessary work in the preconditioner 428 during the repeated solution of linear systems of the same size. The 429 available options are 430 $ SAME_PRECONDITIONER - 431 $ Pmat is identical during successive linear solves. 432 $ This option is intended for folks who are using 433 $ different Amat and Pmat matrices and want to reuse the 434 $ same preconditioner matrix. For example, this option 435 $ saves work by not recomputing incomplete factorization 436 $ for ILU/ICC preconditioners. 437 $ SAME_NONZERO_PATTERN - 438 $ Pmat has the same nonzero structure during 439 $ successive linear solves. 440 $ DIFFERENT_NONZERO_PATTERN - 441 $ Pmat does not have the same nonzero structure. 442 443 All future calls to KSPSetOperators() must use the same size matrices! 444 445 Passing a PETSC_NULL for Amat or Pmat removes the matrix that is currently used. 446 447 If you wish to replace either Amat or Pmat but leave the other one untouched then 448 first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference() 449 on it and then pass it back in in your call to KSPSetOperators(). 450 451 Caution: 452 If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion 453 and does not check the structure of the matrix. If you erroneously 454 claim that the structure is the same when it actually is not, the new 455 preconditioner will not function correctly. Thus, use this optimization 456 feature carefully! 457 458 If in doubt about whether your preconditioner matrix has changed 459 structure or not, use the flag DIFFERENT_NONZERO_PATTERN. 460 461 Level: beginner 462 463 Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators 464 are created in PC and returned to the user. In this case, if both operators 465 mat and pmat are requested, two DIFFERENT operators will be returned. If 466 only one is requested both operators in the PC will be the same (i.e. as 467 if one had called KSP/PCSetOperators() with the same argument for both Mats). 468 The user must set the sizes of the returned matrices and their type etc just 469 as if the user created them with MatCreate(). For example, 470 471 $ KSP/PCGetOperators(ksp/pc,&mat,PETSC_NULL,PETSC_NULL); is equivalent to 472 $ set size, type, etc of mat 473 474 $ MatCreate(comm,&mat); 475 $ KSP/PCSetOperators(ksp/pc,mat,mat,SAME_NONZERO_PATTERN); 476 $ PetscObjectDereference((PetscObject)mat); 477 $ set size, type, etc of mat 478 479 and 480 481 $ KSP/PCGetOperators(ksp/pc,&mat,&pmat,PETSC_NULL); is equivalent to 482 $ set size, type, etc of mat and pmat 483 484 $ MatCreate(comm,&mat); 485 $ MatCreate(comm,&pmat); 486 $ KSP/PCSetOperators(ksp/pc,mat,pmat,SAME_NONZERO_PATTERN); 487 $ PetscObjectDereference((PetscObject)mat); 488 $ PetscObjectDereference((PetscObject)pmat); 489 $ set size, type, etc of mat and pmat 490 491 The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy 492 of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely 493 managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look 494 at this is when you create a SNES you do not NEED to create a KSP and attach it to 495 the SNES object (the SNES object manages it for you). Similarly when you create a KSP 496 you do not need to attach a PC to it (the KSP object manages the PC object for you). 497 Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when 498 it can be created for you? 499 500 .keywords: KSP, set, operators, matrix, preconditioner, linear system 501 502 .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPGetOperators() 503 @*/ 504 PetscErrorCode KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat,MatStructure flag) 505 { 506 MatNullSpace nullsp; 507 PetscErrorCode ierr; 508 509 PetscFunctionBegin; 510 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 511 if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); 512 if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); 513 if (Amat) PetscCheckSameComm(ksp,1,Amat,2); 514 if (Pmat) PetscCheckSameComm(ksp,1,Pmat,3); 515 if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);} 516 ierr = PCSetOperators(ksp->pc,Amat,Pmat,flag);CHKERRQ(ierr); 517 if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX; /* so that next solve call will call PCSetUp() on new matrix */ 518 if (ksp->guess) { 519 ierr = KSPFischerGuessReset(ksp->guess);CHKERRQ(ierr); 520 } 521 if (Pmat) { 522 ierr = MatGetNullSpace(Pmat, &nullsp);CHKERRQ(ierr); 523 if (nullsp) { 524 ierr = KSPSetNullSpace(ksp, nullsp);CHKERRQ(ierr); 525 } 526 } 527 PetscFunctionReturn(0); 528 } 529 530 #undef __FUNCT__ 531 #define __FUNCT__ "KSPGetOperators" 532 /*@ 533 KSPGetOperators - Gets the matrix associated with the linear system 534 and a (possibly) different one associated with the preconditioner. 535 536 Collective on KSP and Mat 537 538 Input Parameter: 539 . ksp - the KSP context 540 541 Output Parameters: 542 + Amat - the matrix associated with the linear system 543 . Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat. 544 - flag - flag indicating information about the preconditioner matrix structure 545 during successive linear solves. This flag is ignored the first time a 546 linear system is solved, and thus is irrelevant when solving just one linear 547 system. 548 549 Level: intermediate 550 551 Notes: DOES NOT increase the reference counts of the matrix, so you should NOT destroy them. 552 553 .keywords: KSP, set, get, operators, matrix, preconditioner, linear system 554 555 .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPSetOperators(), KSPGetOperatorsSet() 556 @*/ 557 PetscErrorCode KSPGetOperators(KSP ksp,Mat *Amat,Mat *Pmat,MatStructure *flag) 558 { 559 PetscErrorCode ierr; 560 561 PetscFunctionBegin; 562 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 563 if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);} 564 ierr = PCGetOperators(ksp->pc,Amat,Pmat,flag);CHKERRQ(ierr); 565 PetscFunctionReturn(0); 566 } 567 568 #undef __FUNCT__ 569 #define __FUNCT__ "KSPGetOperatorsSet" 570 /*@C 571 KSPGetOperatorsSet - Determines if the matrix associated with the linear system and 572 possibly a different one associated with the preconditioner have been set in the KSP. 573 574 Not collective, though the results on all processes should be the same 575 576 Input Parameter: 577 . pc - the KSP context 578 579 Output Parameters: 580 + mat - the matrix associated with the linear system was set 581 - pmat - matrix associated with the preconditioner was set, usually the same 582 583 Level: intermediate 584 585 .keywords: KSP, get, operators, matrix, linear system 586 587 .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators(), PCGetOperatorsSet() 588 @*/ 589 PetscErrorCode KSPGetOperatorsSet(KSP ksp,PetscBool *mat,PetscBool *pmat) 590 { 591 PetscErrorCode ierr; 592 593 PetscFunctionBegin; 594 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 595 if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);} 596 ierr = PCGetOperatorsSet(ksp->pc,mat,pmat);CHKERRQ(ierr); 597 PetscFunctionReturn(0); 598 } 599 600 #undef __FUNCT__ 601 #define __FUNCT__ "KSPCreate" 602 /*@ 603 KSPCreate - Creates the default KSP context. 604 605 Collective on MPI_Comm 606 607 Input Parameter: 608 . comm - MPI communicator 609 610 Output Parameter: 611 . ksp - location to put the KSP context 612 613 Notes: 614 The default KSP type is GMRES with a restart of 30, using modified Gram-Schmidt 615 orthogonalization. 616 617 Level: beginner 618 619 .keywords: KSP, create, context 620 621 .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP 622 @*/ 623 PetscErrorCode KSPCreate(MPI_Comm comm,KSP *inksp) 624 { 625 KSP ksp; 626 PetscErrorCode ierr; 627 void *ctx; 628 629 PetscFunctionBegin; 630 PetscValidPointer(inksp,2); 631 *inksp = 0; 632 #ifndef PETSC_USE_DYNAMIC_LIBRARIES 633 ierr = KSPInitializePackage(PETSC_NULL);CHKERRQ(ierr); 634 #endif 635 636 ierr = PetscHeaderCreate(ksp,_p_KSP,struct _KSPOps,KSP_CLASSID,-1,"KSP","Krylov Method","KSP",comm,KSPDestroy,KSPView);CHKERRQ(ierr); 637 638 ksp->max_it = 10000; 639 ksp->pc_side = PC_SIDE_DEFAULT; 640 ksp->rtol = 1.e-5; 641 ksp->abstol = 1.e-50; 642 ksp->divtol = 1.e4; 643 644 ksp->chknorm = -1; 645 ksp->normtype = KSP_NORM_DEFAULT; 646 ksp->rnorm = 0.0; 647 ksp->its = 0; 648 ksp->guess_zero = PETSC_TRUE; 649 ksp->calc_sings = PETSC_FALSE; 650 ksp->res_hist = PETSC_NULL; 651 ksp->res_hist_alloc = PETSC_NULL; 652 ksp->res_hist_len = 0; 653 ksp->res_hist_max = 0; 654 ksp->res_hist_reset = PETSC_TRUE; 655 ksp->numbermonitors = 0; 656 657 ierr = KSPDefaultConvergedCreate(&ctx);CHKERRQ(ierr); 658 ierr = KSPSetConvergenceTest(ksp,KSPDefaultConverged,ctx,KSPDefaultConvergedDestroy);CHKERRQ(ierr); 659 ksp->ops->buildsolution = KSPDefaultBuildSolution; 660 ksp->ops->buildresidual = KSPDefaultBuildResidual; 661 #if defined(PETSC_HAVE_AMS) 662 ((PetscObject)ksp)->bops->publish = KSPPublish_Petsc; 663 #endif 664 665 ksp->vec_sol = 0; 666 ksp->vec_rhs = 0; 667 ksp->pc = 0; 668 ksp->data = 0; 669 ksp->nwork = 0; 670 ksp->work = 0; 671 ksp->reason = KSP_CONVERGED_ITERATING; 672 ksp->setupstage = KSP_SETUP_NEW; 673 674 ierr = KSPNormSupportTableReset_Private(ksp);CHKERRQ(ierr); 675 676 *inksp = ksp; 677 PetscFunctionReturn(0); 678 } 679 680 #undef __FUNCT__ 681 #define __FUNCT__ "KSPSetType" 682 /*@C 683 KSPSetType - Builds KSP for a particular solver. 684 685 Logically Collective on KSP 686 687 Input Parameters: 688 + ksp - the Krylov space context 689 - type - a known method 690 691 Options Database Key: 692 . -ksp_type <method> - Sets the method; use -help for a list 693 of available methods (for instance, cg or gmres) 694 695 Notes: 696 See "petsc/include/petscksp.h" for available methods (for instance, 697 KSPCG or KSPGMRES). 698 699 Normally, it is best to use the KSPSetFromOptions() command and 700 then set the KSP type from the options database rather than by using 701 this routine. Using the options database provides the user with 702 maximum flexibility in evaluating the many different Krylov methods. 703 The KSPSetType() routine is provided for those situations where it 704 is necessary to set the iterative solver independently of the command 705 line or options database. This might be the case, for example, when 706 the choice of iterative solver changes during the execution of the 707 program, and the user's application is taking responsibility for 708 choosing the appropriate method. In other words, this routine is 709 not for beginners. 710 711 Level: intermediate 712 713 .keywords: KSP, set, method 714 715 .seealso: PCSetType(), KSPType 716 717 @*/ 718 PetscErrorCode KSPSetType(KSP ksp, KSPType type) 719 { 720 PetscErrorCode ierr,(*r)(KSP); 721 PetscBool match; 722 723 PetscFunctionBegin; 724 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 725 PetscValidCharPointer(type,2); 726 727 ierr = PetscObjectTypeCompare((PetscObject)ksp,type,&match);CHKERRQ(ierr); 728 if (match) PetscFunctionReturn(0); 729 730 ierr = PetscFListFind(KSPList,((PetscObject)ksp)->comm,type,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 731 if (!r) SETERRQ1(((PetscObject)ksp)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested KSP type %s",type); 732 /* Destroy the previous private KSP context */ 733 if (ksp->ops->destroy) { 734 ierr = (*ksp->ops->destroy)(ksp);CHKERRQ(ierr); 735 ksp->ops->destroy = PETSC_NULL; 736 } 737 /* Reinitialize function pointers in KSPOps structure */ 738 ierr = PetscMemzero(ksp->ops,sizeof(struct _KSPOps));CHKERRQ(ierr); 739 ksp->ops->buildsolution = KSPDefaultBuildSolution; 740 ksp->ops->buildresidual = KSPDefaultBuildResidual; 741 ierr = KSPNormSupportTableReset_Private(ksp);CHKERRQ(ierr); 742 /* Call the KSPCreate_XXX routine for this particular Krylov solver */ 743 ksp->setupstage = KSP_SETUP_NEW; 744 ierr = PetscObjectChangeTypeName((PetscObject)ksp,type);CHKERRQ(ierr); 745 ierr = (*r)(ksp);CHKERRQ(ierr); 746 #if defined(PETSC_HAVE_AMS) 747 if (PetscAMSPublishAll) { 748 ierr = PetscObjectAMSPublish((PetscObject)ksp);CHKERRQ(ierr); 749 } 750 #endif 751 PetscFunctionReturn(0); 752 } 753 754 #undef __FUNCT__ 755 #define __FUNCT__ "KSPRegisterDestroy" 756 /*@ 757 KSPRegisterDestroy - Frees the list of KSP methods that were 758 registered by KSPRegisterDynamic(). 759 760 Not Collective 761 762 Level: advanced 763 764 .keywords: KSP, register, destroy 765 766 .seealso: KSPRegisterDynamic(), KSPRegisterAll() 767 @*/ 768 PetscErrorCode KSPRegisterDestroy(void) 769 { 770 PetscErrorCode ierr; 771 772 PetscFunctionBegin; 773 ierr = PetscFListDestroy(&KSPList);CHKERRQ(ierr); 774 KSPRegisterAllCalled = PETSC_FALSE; 775 PetscFunctionReturn(0); 776 } 777 778 #undef __FUNCT__ 779 #define __FUNCT__ "KSPGetType" 780 /*@C 781 KSPGetType - Gets the KSP type as a string from the KSP object. 782 783 Not Collective 784 785 Input Parameter: 786 . ksp - Krylov context 787 788 Output Parameter: 789 . name - name of KSP method 790 791 Level: intermediate 792 793 .keywords: KSP, get, method, name 794 795 .seealso: KSPSetType() 796 @*/ 797 PetscErrorCode KSPGetType(KSP ksp,KSPType *type) 798 { 799 PetscFunctionBegin; 800 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 801 PetscValidPointer(type,2); 802 *type = ((PetscObject)ksp)->type_name; 803 PetscFunctionReturn(0); 804 } 805 806 #undef __FUNCT__ 807 #define __FUNCT__ "KSPRegister" 808 /*@C 809 KSPRegister - See KSPRegisterDynamic() 810 811 Level: advanced 812 @*/ 813 PetscErrorCode KSPRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(KSP)) 814 { 815 PetscErrorCode ierr; 816 char fullname[PETSC_MAX_PATH_LEN]; 817 818 PetscFunctionBegin; 819 ierr = PetscFListConcat(path,name,fullname);CHKERRQ(ierr); 820 ierr = PetscFListAdd(&KSPList,sname,fullname,(void (*)(void))function);CHKERRQ(ierr); 821 PetscFunctionReturn(0); 822 } 823 824 #undef __FUNCT__ 825 #define __FUNCT__ "KSPSetNullSpace" 826 /*@ 827 KSPSetNullSpace - Sets the null space of the operator 828 829 Logically Collective on KSP 830 831 Input Parameters: 832 + ksp - the Krylov space object 833 - nullsp - the null space of the operator 834 835 Level: advanced 836 837 .seealso: KSPSetOperators(), MatNullSpaceCreate(), KSPGetNullSpace(), MatSetNullSpace() 838 @*/ 839 PetscErrorCode KSPSetNullSpace(KSP ksp,MatNullSpace nullsp) 840 { 841 PetscErrorCode ierr; 842 843 PetscFunctionBegin; 844 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 845 PetscValidHeaderSpecific(nullsp,MAT_NULLSPACE_CLASSID,2); 846 ierr = PetscObjectReference((PetscObject)nullsp);CHKERRQ(ierr); 847 if (ksp->nullsp) { ierr = MatNullSpaceDestroy(&ksp->nullsp);CHKERRQ(ierr); } 848 ksp->nullsp = nullsp; 849 PetscFunctionReturn(0); 850 } 851 852 #undef __FUNCT__ 853 #define __FUNCT__ "KSPGetNullSpace" 854 /*@ 855 KSPGetNullSpace - Gets the null space of the operator 856 857 Not Collective 858 859 Input Parameters: 860 + ksp - the Krylov space object 861 - nullsp - the null space of the operator 862 863 Level: advanced 864 865 .seealso: KSPSetOperators(), MatNullSpaceCreate(), KSPSetNullSpace() 866 @*/ 867 PetscErrorCode KSPGetNullSpace(KSP ksp,MatNullSpace *nullsp) 868 { 869 PetscFunctionBegin; 870 PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); 871 PetscValidPointer(nullsp,2); 872 *nullsp = ksp->nullsp; 873 PetscFunctionReturn(0); 874 } 875 876