Actual source code: itcreate.c

petsc-3.4.5 2014-06-29
  2: /*
  3:      The basic KSP routines, Create, View etc. are here.
  4: */
  5: #include <petsc-private/kspimpl.h>      /*I "petscksp.h" I*/

  7: /* Logging support */
  8: PetscClassId  KSP_CLASSID;
  9: PetscClassId  DMKSP_CLASSID;
 10: PetscLogEvent KSP_GMRESOrthogonalization, KSP_SetUp, KSP_Solve;

 12: /*
 13:    Contains the list of registered KSP routines
 14: */
 15: PetscFunctionList KSPList              = 0;
 16: PetscBool         KSPRegisterAllCalled = PETSC_FALSE;

 20: /*@C
 21:   KSPLoad - Loads a KSP that has been stored in binary  with KSPView().

 23:   Collective on PetscViewer

 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()

 30:    Level: intermediate

 32:   Notes:
 33:    The type is determined by the data in the file, any type set into the KSP before this call is ignored.

 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

 44: .seealso: PetscViewerBinaryOpen(), KSPView(), MatLoad(), VecLoad()
 45: @*/
 46: PetscErrorCode  KSPLoad(KSP newdm, PetscViewer viewer)
 47: {
 49:   PetscBool      isbinary;
 50:   PetscInt       classid;
 51:   char           type[256];
 52:   PC             pc;

 57:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
 58:   if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

 60:   PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
 61:   if (classid != KSP_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not KSP next in file");
 62:   PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
 63:   KSPSetType(newdm, type);
 64:   if (newdm->ops->load) {
 65:     (*newdm->ops->load)(newdm,viewer);
 66:   }
 67:   KSPGetPC(newdm,&pc);
 68:   PCLoad(pc,viewer);
 69:   return(0);
 70: }

 72: #include <petscdraw.h>
 73: #if defined(PETSC_HAVE_AMS)
 74: #include <petscviewerams.h>
 75: #endif
 78: /*@C
 79:    KSPView - Prints the KSP data structure.

 81:    Collective on KSP

 83:    Input Parameters:
 84: +  ksp - the Krylov space context
 85: -  viewer - visualization context

 87:    Options Database Keys:
 88: .  -ksp_view - print the ksp data structure at the end of a KSPSolve call

 90:    Note:
 91:    The available visualization contexts include
 92: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 93: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 94:          output where only the first processor opens
 95:          the file.  All other processors send their
 96:          data to the first processor to print.

 98:    The user can open an alternative visualization context with
 99:    PetscViewerASCIIOpen() - output to a specified file.

101:    Level: beginner

103: .keywords: KSP, view

105: .seealso: PCView(), PetscViewerASCIIOpen()
106: @*/
107: PetscErrorCode  KSPView(KSP ksp,PetscViewer viewer)
108: {
110:   PetscBool      iascii,isbinary,isdraw;
111: #if defined(PETSC_HAVE_AMS)
112:   PetscBool      isams;
113: #endif

117:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)ksp));

121:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
122:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
123:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
124: #if defined(PETSC_HAVE_AMS)
125:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERAMS,&isams);
126: #endif
127:   if (iascii) {
128:     PetscObjectPrintClassNamePrefixType((PetscObject)ksp,viewer,"KSP Object");
129:     if (ksp->ops->view) {
130:       PetscViewerASCIIPushTab(viewer);
131:       (*ksp->ops->view)(ksp,viewer);
132:       PetscViewerASCIIPopTab(viewer);
133:     }
134:     if (ksp->guess_zero) {
135:       PetscViewerASCIIPrintf(viewer,"  maximum iterations=%D, initial guess is zero\n",ksp->max_it);
136:     } else {
137:       PetscViewerASCIIPrintf(viewer,"  maximum iterations=%D\n", ksp->max_it);
138:     }
139:     if (ksp->guess_knoll) {PetscViewerASCIIPrintf(viewer,"  using preconditioner applied to right hand side for initial guess\n");}
140:     PetscViewerASCIIPrintf(viewer,"  tolerances:  relative=%G, absolute=%G, divergence=%G\n",ksp->rtol,ksp->abstol,ksp->divtol);
141:     if (ksp->pc_side == PC_RIGHT) {
142:       PetscViewerASCIIPrintf(viewer,"  right preconditioning\n");
143:     } else if (ksp->pc_side == PC_SYMMETRIC) {
144:       PetscViewerASCIIPrintf(viewer,"  symmetric preconditioning\n");
145:     } else {
146:       PetscViewerASCIIPrintf(viewer,"  left preconditioning\n");
147:     }
148:     if (ksp->guess) {PetscViewerASCIIPrintf(viewer,"  using Fischers initial guess method %D with size %D\n",ksp->guess->method,ksp->guess->maxl);}
149:     if (ksp->dscale) {PetscViewerASCIIPrintf(viewer,"  diagonally scaled system\n");}
150:     if (ksp->nullsp) {PetscViewerASCIIPrintf(viewer,"  has attached null space\n");}
151:     if (!ksp->guess_zero) {PetscViewerASCIIPrintf(viewer,"  using nonzero initial guess\n");}
152:     PetscViewerASCIIPrintf(viewer,"  using %s norm type for convergence test\n",KSPNormTypes[ksp->normtype]);
153:   } else if (isbinary) {
154:     PetscInt    classid = KSP_FILE_CLASSID;
155:     MPI_Comm    comm;
156:     PetscMPIInt rank;
157:     char        type[256];

159:     PetscObjectGetComm((PetscObject)ksp,&comm);
160:     MPI_Comm_rank(comm,&rank);
161:     if (!rank) {
162:       PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
163:       PetscStrncpy(type,((PetscObject)ksp)->type_name,256);
164:       PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
165:     }
166:     if (ksp->ops->view) {
167:       (*ksp->ops->view)(ksp,viewer);
168:     }
169:   } else if (isdraw) {
170:     PetscDraw draw;
171:     char      str[36];
172:     PetscReal x,y,bottom,h;
173:     PetscBool flg;

175:     PetscViewerDrawGetDraw(viewer,0,&draw);
176:     PetscDrawGetCurrentPoint(draw,&x,&y);
177:     PetscObjectTypeCompare((PetscObject)ksp,KSPPREONLY,&flg);
178:     if (!flg) {
179:       PetscStrcpy(str,"KSP: ");
180:       PetscStrcat(str,((PetscObject)ksp)->type_name);
181:       PetscDrawBoxedString(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
182:       bottom = y - h;
183:     } else {
184:       bottom = y;
185:     }
186:     PetscDrawPushCurrentPoint(draw,x,bottom);
187: #if defined(PETSC_HAVE_AMS)
188:   } else if (isams) {
189:     if (((PetscObject)ksp)->amsmem == -1) {
190:       PetscObjectViewAMS((PetscObject)ksp,viewer);
191:       PetscStackCallAMS(AMS_Memory_take_access,(((PetscObject)ksp)->amsmem));
192:       PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ksp)->amsmem,"its",&ksp->its,1,AMS_INT,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
193:       if (!ksp->res_hist) {
194:         KSPSetResidualHistory(ksp,NULL,PETSC_DECIDE,PETSC_FALSE);
195:       }
196:       PetscStackCallAMS(AMS_Memory_add_field,(((PetscObject)ksp)->amsmem,"res_hist",ksp->res_hist,10,AMS_DOUBLE,AMS_READ,AMS_COMMON,AMS_REDUCT_UNDEF));
197:       PetscStackCallAMS(AMS_Memory_grant_access,(((PetscObject)ksp)->amsmem));
198:     }
199: #endif
200:   } else if (ksp->ops->view) {
201:     (*ksp->ops->view)(ksp,viewer);
202:   }
203:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
204:   PCView(ksp->pc,viewer);
205:   if (isdraw) {
206:     PetscDraw draw;
207:     PetscViewerDrawGetDraw(viewer,0,&draw);
208:     PetscDrawPopCurrentPoint(draw);
209:   }
210:   return(0);
211: }


216: /*@
217:    KSPSetNormType - Sets the norm that is used for convergence testing.

219:    Logically Collective on KSP

221:    Input Parameter:
222: +  ksp - Krylov solver context
223: -  normtype - one of
224: $   KSP_NORM_NONE - skips computing the norm, this should only be used if you are using
225: $                 the Krylov method as a smoother with a fixed small number of iterations.
226: $                 Implicitly sets KSPSkipConverged as KSP convergence test.
227: $                 Supported only by CG, Richardson, Bi-CG-stab, CR, and CGS methods.
228: $   KSP_NORM_PRECONDITIONED - the default for left preconditioned solves, uses the l2 norm
229: $                 of the preconditioned residual
230: $   KSP_NORM_UNPRECONDITIONED - uses the l2 norm of the true b - Ax residual, supported only by
231: $                 CG, CHEBYSHEV, and RICHARDSON, automatically true for right (see KSPSetPCSide())
232: $                 preconditioning..
233: $   KSP_NORM_NATURAL - supported  by KSPCG, KSPCR, KSPCGNE, KSPCGS


236:    Options Database Key:
237: .   -ksp_norm_type <none,preconditioned,unpreconditioned,natural>

239:    Notes:
240:    Currently only works with the CG, Richardson, Bi-CG-stab, CR, and CGS methods.

242:    Level: advanced

244: .keywords: KSP, create, context, norms

246: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetCheckNormIteration()
247: @*/
248: PetscErrorCode  KSPSetNormType(KSP ksp,KSPNormType normtype)
249: {

255:   ksp->normtype = normtype;
256:   if (normtype == KSP_NORM_NONE) {
257:     KSPSetConvergenceTest(ksp,KSPSkipConverged,0,0);
258:     PetscInfo(ksp,"Warning: setting KSPNormType to skip computing the norm\n\
259:  KSP convergence test is implicitly set to KSPSkipConverged\n");
260:   }
261:   return(0);
262: }

266: /*@
267:    KSPSetCheckNormIteration - Sets the first iteration at which the norm of the residual will be
268:      computed and used in the convergence test.

270:    Logically Collective on KSP

272:    Input Parameter:
273: +  ksp - Krylov solver context
274: -  it  - use -1 to check at all iterations

276:    Notes:
277:    Currently only works with KSPCG, KSPBCGS and KSPIBCGS

279:    Use KSPSetNormType(ksp,KSP_NORM_NONE) to never check the norm

281:    On steps where the norm is not computed, the previous norm is still in the variable, so if you run with, for example,
282:     -ksp_monitor the residual norm will appear to be unchanged for several iterations (though it is not really unchanged).
283:    Level: advanced

285: .keywords: KSP, create, context, norms

287: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetNormType()
288: @*/
289: PetscErrorCode  KSPSetCheckNormIteration(KSP ksp,PetscInt it)
290: {
294:   ksp->chknorm = it;
295:   return(0);
296: }

300: /*@
301:    KSPSetLagNorm - Lags the residual norm calculation so that it is computed as part of the MPI_Allreduce() for
302:    computing the inner products for the next iteration.  This can reduce communication costs at the expense of doing
303:    one additional iteration.


306:    Logically Collective on KSP

308:    Input Parameter:
309: +  ksp - Krylov solver context
310: -  flg - PETSC_TRUE or PETSC_FALSE

312:    Options Database Keys:
313: .  -ksp_lag_norm - lag the calculated residual norm

315:    Notes:
316:    Currently only works with KSPIBCGS.

318:    Use KSPSetNormType(ksp,KSP_NORM_NONE) to never check the norm

320:    If you lag the norm and run with, for example, -ksp_monitor, the residual norm reported will be the lagged one.
321:    Level: advanced

323: .keywords: KSP, create, context, norms

325: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged(), KSPSetNormType(), KSPSetCheckNormIteration()
326: @*/
327: PetscErrorCode  KSPSetLagNorm(KSP ksp,PetscBool flg)
328: {
332:   ksp->lagnorm = flg;
333:   return(0);
334: }

338: /*@
339:    KSPSetSupportedNorm - Sets a norm and preconditioner side supported by a KSP

341:    Logically Collective

343:    Input Arguments:
344: +  ksp - Krylov method
345: .  normtype - supported norm type
346: .  pcside - preconditioner side that can be used with this norm
347: -  preference - integer preference for this combination, larger values have higher priority

349:    Level: developer

351:    Notes:
352:    This function should be called from the implementation files KSPCreate_XXX() to declare
353:    which norms and preconditioner sides are supported. Users should not need to call this
354:    function.

356:    KSP_NORM_NONE is supported by default with all KSP methods and any PC side. If a KSP explicitly does not support
357:    KSP_NORM_NONE, it should set this by setting priority=0.

359: .seealso: KSPSetNormType(), KSPSetPCSide()
360: @*/
361: PetscErrorCode KSPSetSupportedNorm(KSP ksp,KSPNormType normtype,PCSide pcside,PetscInt priority)
362: {

366:   ksp->normsupporttable[normtype][pcside] = priority;
367:   return(0);
368: }

372: PetscErrorCode KSPNormSupportTableReset_Private(KSP ksp)
373: {

377:   PetscMemzero(ksp->normsupporttable,sizeof(ksp->normsupporttable));
378:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_LEFT,1);
379:   KSPSetSupportedNorm(ksp,KSP_NORM_NONE,PC_RIGHT,1);
380:   return(0);
381: }

385: PetscErrorCode KSPSetUpNorms_Private(KSP ksp,KSPNormType *normtype,PCSide *pcside)
386: {
387:   PetscInt i,j,best,ibest = 0,jbest = 0;

390:   best = 0;
391:   for (i=0; i<KSP_NORM_MAX; i++) {
392:     for (j=0; j<PC_SIDE_MAX; j++) {
393:       if ((ksp->normtype == KSP_NORM_DEFAULT || ksp->normtype == i)
394:           && (ksp->pc_side == PC_SIDE_DEFAULT || ksp->pc_side == j)
395:           && (ksp->normsupporttable[i][j] > best)) {
396:         if (ksp->normtype == KSP_NORM_DEFAULT && i == KSP_NORM_NONE && ksp->normsupporttable[i][j] <= 1) {
397:           continue; /* Skip because we don't want to default to no norms unless set by the KSP (preonly). */
398:         }
399:         best  = ksp->normsupporttable[i][j];
400:         ibest = i;
401:         jbest = j;
402:       }
403:     }
404:   }
405:   if (best < 1) {
406:     if (ksp->normtype == KSP_NORM_DEFAULT && ksp->pc_side == PC_SIDE_DEFAULT) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"The %s KSP implementation did not call KSPSetSupportedNorm()",((PetscObject)ksp)->type_name);
407:     if (ksp->normtype == KSP_NORM_DEFAULT) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP %s does not support %s",((PetscObject)ksp)->type_name,PCSides[ksp->pc_side]);
408:     if (ksp->pc_side == PC_SIDE_DEFAULT) SETERRQ2(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP %s does not support %s",((PetscObject)ksp)->type_name,KSPNormTypes[ksp->normtype]);
409:     SETERRQ3(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"KSP %s does not support %s with %s",((PetscObject)ksp)->type_name,KSPNormTypes[ksp->normtype],PCSides[ksp->pc_side]);
410:   }
411:   *normtype = (KSPNormType)ibest;
412:   *pcside   = (PCSide)jbest;
413:   return(0);
414: }

418: /*@
419:    KSPGetNormType - Gets the norm that is used for convergence testing.

421:    Not Collective

423:    Input Parameter:
424: .  ksp - Krylov solver context

426:    Output Parameter:
427: .  normtype - norm that is used for convergence testing

429:    Level: advanced

431: .keywords: KSP, create, context, norms

433: .seealso: KSPNormType, KSPSetNormType(), KSPSkipConverged()
434: @*/
435: PetscErrorCode  KSPGetNormType(KSP ksp, KSPNormType *normtype)
436: {

442:   KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);
443:   *normtype = ksp->normtype;
444:   return(0);
445: }

447: #if defined(PETSC_HAVE_AMS)
448: #include <petscviewerams.h>
449: #endif

453: /*@
454:    KSPSetOperators - Sets the matrix associated with the linear system
455:    and a (possibly) different one associated with the preconditioner.

457:    Collective on KSP and Mat

459:    Input Parameters:
460: +  ksp - the KSP context
461: .  Amat - the matrix that defines the linear system
462: .  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
463: -  flag - flag indicating information about the preconditioner matrix structure
464:    during successive linear solves.  This flag is ignored the first time a
465:    linear system is solved, and thus is irrelevant when solving just one linear
466:    system.

468:    Notes:
469:    The flag can be used to eliminate unnecessary work in the preconditioner
470:    during the repeated solution of linear systems of the same size.  The
471:    available options are
472: $    SAME_PRECONDITIONER -
473: $      Pmat is identical during successive linear solves.
474: $      This option is intended for folks who are using
475: $      different Amat and Pmat matrices and want to reuse the
476: $      same preconditioner matrix.  For example, this option
477: $      saves work by not recomputing incomplete factorization
478: $      for ILU/ICC preconditioners.
479: $    SAME_NONZERO_PATTERN -
480: $      Pmat has the same nonzero structure during
481: $      successive linear solves.
482: $    DIFFERENT_NONZERO_PATTERN -
483: $      Pmat does not have the same nonzero structure.

485:     All future calls to KSPSetOperators() must use the same size matrices!

487:     Passing a NULL for Amat or Pmat removes the matrix that is currently used.

489:     If you wish to replace either Amat or Pmat but leave the other one untouched then
490:     first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
491:     on it and then pass it back in in your call to KSPSetOperators().

493:     Caution:
494:     If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion
495:     and does not check the structure of the matrix.  If you erroneously
496:     claim that the structure is the same when it actually is not, the new
497:     preconditioner will not function correctly.  Thus, use this optimization
498:     feature carefully!

500:     If in doubt about whether your preconditioner matrix has changed
501:     structure or not, use the flag DIFFERENT_NONZERO_PATTERN.

503:     Level: beginner

505:    Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
506:       are created in PC and returned to the user. In this case, if both operators
507:       mat and pmat are requested, two DIFFERENT operators will be returned. If
508:       only one is requested both operators in the PC will be the same (i.e. as
509:       if one had called KSP/PCSetOperators() with the same argument for both Mats).
510:       The user must set the sizes of the returned matrices and their type etc just
511:       as if the user created them with MatCreate(). For example,

513: $         KSP/PCGetOperators(ksp/pc,&mat,NULL,NULL); is equivalent to
514: $           set size, type, etc of mat

516: $         MatCreate(comm,&mat);
517: $         KSP/PCSetOperators(ksp/pc,mat,mat,SAME_NONZERO_PATTERN);
518: $         PetscObjectDereference((PetscObject)mat);
519: $           set size, type, etc of mat

521:      and

523: $         KSP/PCGetOperators(ksp/pc,&mat,&pmat,NULL); is equivalent to
524: $           set size, type, etc of mat and pmat

526: $         MatCreate(comm,&mat);
527: $         MatCreate(comm,&pmat);
528: $         KSP/PCSetOperators(ksp/pc,mat,pmat,SAME_NONZERO_PATTERN);
529: $         PetscObjectDereference((PetscObject)mat);
530: $         PetscObjectDereference((PetscObject)pmat);
531: $           set size, type, etc of mat and pmat

533:     The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
534:     of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
535:     managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
536:     at this is when you create a SNES you do not NEED to create a KSP and attach it to
537:     the SNES object (the SNES object manages it for you). Similarly when you create a KSP
538:     you do not need to attach a PC to it (the KSP object manages the PC object for you).
539:     Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
540:     it can be created for you?

542: .keywords: KSP, set, operators, matrix, preconditioner, linear system

544: .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPGetOperators()
545: @*/
546: PetscErrorCode  KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat,MatStructure flag)
547: {
548:   MatNullSpace   nullsp;

557:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
558:   PCSetOperators(ksp->pc,Amat,Pmat,flag);
559:   if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;  /* so that next solve call will call PCSetUp() on new matrix */
560:   if (ksp->guess) {
561:     KSPFischerGuessReset(ksp->guess);
562:   }
563:   if (Pmat) {
564:     MatGetNullSpace(Pmat, &nullsp);
565:     if (nullsp) {
566:       KSPSetNullSpace(ksp, nullsp);
567:     }
568:   }
569:   return(0);
570: }

574: /*@
575:    KSPGetOperators - Gets the matrix associated with the linear system
576:    and a (possibly) different one associated with the preconditioner.

578:    Collective on KSP and Mat

580:    Input Parameter:
581: .  ksp - the KSP context

583:    Output Parameters:
584: +  Amat - the matrix that defines the linear system
585: .  Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
586: -  flag - flag indicating information about the preconditioner matrix structure
587:    during successive linear solves.  This flag is ignored the first time a
588:    linear system is solved, and thus is irrelevant when solving just one linear
589:    system.

591:     Level: intermediate

593:    Notes: DOES NOT increase the reference counts of the matrix, so you should NOT destroy them.

595: .keywords: KSP, set, get, operators, matrix, preconditioner, linear system

597: .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPSetOperators(), KSPGetOperatorsSet()
598: @*/
599: PetscErrorCode  KSPGetOperators(KSP ksp,Mat *Amat,Mat *Pmat,MatStructure *flag)
600: {

605:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
606:   PCGetOperators(ksp->pc,Amat,Pmat,flag);
607:   return(0);
608: }

612: /*@C
613:    KSPGetOperatorsSet - Determines if the matrix associated with the linear system and
614:    possibly a different one associated with the preconditioner have been set in the KSP.

616:    Not collective, though the results on all processes should be the same

618:    Input Parameter:
619: .  pc - the KSP context

621:    Output Parameters:
622: +  mat - the matrix associated with the linear system was set
623: -  pmat - matrix associated with the preconditioner was set, usually the same

625:    Level: intermediate

627: .keywords: KSP, get, operators, matrix, linear system

629: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators(), PCGetOperatorsSet()
630: @*/
631: PetscErrorCode  KSPGetOperatorsSet(KSP ksp,PetscBool  *mat,PetscBool  *pmat)
632: {

637:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
638:   PCGetOperatorsSet(ksp->pc,mat,pmat);
639:   return(0);
640: }

644: /*@C
645:    KSPSetPreSolve - Sets a function that is called before every KSPSolve() is started

647:    Logically Collective on KSP

649:    Input Parameters:
650: +   ksp - the solver object
651: .   presolve - the function to call before the solve
652: -   prectx - any context needed by the function

654:    Level: developer

656: .keywords: KSP, create, context

658: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP, KSPSetPostSolve()
659: @*/
660: PetscErrorCode  KSPSetPreSolve(KSP ksp,PetscErrorCode (*presolve)(KSP,Vec,Vec,void*),void *prectx)
661: {
664:   ksp->presolve = presolve;
665:   ksp->prectx   = prectx;
666:   return(0);
667: }

671: /*@C
672:    KSPSetPostSolve - Sets a function that is called after every KSPSolve() completes (whether it converges or not)

674:    Logically Collective on KSP

676:    Input Parameters:
677: +   ksp - the solver object
678: .   postsolve - the function to call after the solve
679: -   postctx - any context needed by the function

681:    Level: developer

683: .keywords: KSP, create, context

685: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP, KSPSetPreSolve()
686: @*/
687: PetscErrorCode  KSPSetPostSolve(KSP ksp,PetscErrorCode (*postsolve)(KSP,Vec,Vec,void*),void *postctx)
688: {
691:   ksp->postsolve = postsolve;
692:   ksp->postctx   = postctx;
693:   return(0);
694: }

698: /*@
699:    KSPCreate - Creates the default KSP context.

701:    Collective on MPI_Comm

703:    Input Parameter:
704: .  comm - MPI communicator

706:    Output Parameter:
707: .  ksp - location to put the KSP context

709:    Notes:
710:    The default KSP type is GMRES with a restart of 30, using modified Gram-Schmidt
711:    orthogonalization.

713:    Level: beginner

715: .keywords: KSP, create, context

717: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP
718: @*/
719: PetscErrorCode  KSPCreate(MPI_Comm comm,KSP *inksp)
720: {
721:   KSP            ksp;
723:   void           *ctx;

727:   *inksp = 0;
728: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
729:   KSPInitializePackage();
730: #endif

732:   PetscHeaderCreate(ksp,_p_KSP,struct _KSPOps,KSP_CLASSID,"KSP","Krylov Method","KSP",comm,KSPDestroy,KSPView);

734:   ksp->max_it  = 10000;
735:   ksp->pc_side = PC_SIDE_DEFAULT;
736:   ksp->rtol    = 1.e-5;
737:   ksp->abstol  = 1.e-50;
738:   ksp->divtol  = 1.e4;

740:   ksp->chknorm        = -1;
741:   ksp->normtype       = KSP_NORM_DEFAULT;
742:   ksp->rnorm          = 0.0;
743:   ksp->its            = 0;
744:   ksp->guess_zero     = PETSC_TRUE;
745:   ksp->calc_sings     = PETSC_FALSE;
746:   ksp->res_hist       = NULL;
747:   ksp->res_hist_alloc = NULL;
748:   ksp->res_hist_len   = 0;
749:   ksp->res_hist_max   = 0;
750:   ksp->res_hist_reset = PETSC_TRUE;
751:   ksp->numbermonitors = 0;

753:   KSPDefaultConvergedCreate(&ctx);
754:   KSPSetConvergenceTest(ksp,KSPDefaultConverged,ctx,KSPDefaultConvergedDestroy);
755:   ksp->ops->buildsolution = KSPBuildSolutionDefault;
756:   ksp->ops->buildresidual = KSPBuildResidualDefault;

758:   ksp->vec_sol    = 0;
759:   ksp->vec_rhs    = 0;
760:   ksp->pc         = 0;
761:   ksp->data       = 0;
762:   ksp->nwork      = 0;
763:   ksp->work       = 0;
764:   ksp->reason     = KSP_CONVERGED_ITERATING;
765:   ksp->setupstage = KSP_SETUP_NEW;

767:   KSPNormSupportTableReset_Private(ksp);

769:   *inksp = ksp;
770:   return(0);
771: }

775: /*@C
776:    KSPSetType - Builds KSP for a particular solver.

778:    Logically Collective on KSP

780:    Input Parameters:
781: +  ksp      - the Krylov space context
782: -  type - a known method

784:    Options Database Key:
785: .  -ksp_type  <method> - Sets the method; use -help for a list
786:     of available methods (for instance, cg or gmres)

788:    Notes:
789:    See "petsc/include/petscksp.h" for available methods (for instance,
790:    KSPCG or KSPGMRES).

792:   Normally, it is best to use the KSPSetFromOptions() command and
793:   then set the KSP type from the options database rather than by using
794:   this routine.  Using the options database provides the user with
795:   maximum flexibility in evaluating the many different Krylov methods.
796:   The KSPSetType() routine is provided for those situations where it
797:   is necessary to set the iterative solver independently of the command
798:   line or options database.  This might be the case, for example, when
799:   the choice of iterative solver changes during the execution of the
800:   program, and the user's application is taking responsibility for
801:   choosing the appropriate method.  In other words, this routine is
802:   not for beginners.

804:   Level: intermediate

806:   Developer Note: KSPRegister() is used to add Krylov types to KSPList from which they
807:   are accessed by KSPSetType().

809: .keywords: KSP, set, method

811: .seealso: PCSetType(), KSPType, KSPRegister(), KSPCreate()

813: @*/
814: PetscErrorCode  KSPSetType(KSP ksp, KSPType type)
815: {
816:   PetscErrorCode ierr,(*r)(KSP);
817:   PetscBool      match;


823:   PetscObjectTypeCompare((PetscObject)ksp,type,&match);
824:   if (match) return(0);

826:    PetscFunctionListFind(KSPList,type,&r);
827:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested KSP type %s",type);
828:   /* Destroy the previous private KSP context */
829:   if (ksp->ops->destroy) {
830:     (*ksp->ops->destroy)(ksp);
831:     ksp->ops->destroy = NULL;
832:   }
833:   /* Reinitialize function pointers in KSPOps structure */
834:   PetscMemzero(ksp->ops,sizeof(struct _KSPOps));
835:   ksp->ops->buildsolution = KSPBuildSolutionDefault;
836:   ksp->ops->buildresidual = KSPBuildResidualDefault;
837:   KSPNormSupportTableReset_Private(ksp);
838:   /* Call the KSPCreate_XXX routine for this particular Krylov solver */
839:   ksp->setupstage = KSP_SETUP_NEW;
840:   PetscObjectChangeTypeName((PetscObject)ksp,type);
841:   (*r)(ksp);
842:   return(0);
843: }

847: /*@C
848:    KSPGetType - Gets the KSP type as a string from the KSP object.

850:    Not Collective

852:    Input Parameter:
853: .  ksp - Krylov context

855:    Output Parameter:
856: .  name - name of KSP method

858:    Level: intermediate

860: .keywords: KSP, get, method, name

862: .seealso: KSPSetType()
863: @*/
864: PetscErrorCode  KSPGetType(KSP ksp,KSPType *type)
865: {
869:   *type = ((PetscObject)ksp)->type_name;
870:   return(0);
871: }

875: /*@C
876:   KSPRegister -  Adds a method to the Krylov subspace solver package.

878:    Not Collective

880:    Input Parameters:
881: +  name_solver - name of a new user-defined solver
882: -  routine_create - routine to create method context

884:    Notes:
885:    KSPRegister() may be called multiple times to add several user-defined solvers.

887:    Sample usage:
888: .vb
889:    KSPRegister("my_solver",MySolverCreate);
890: .ve

892:    Then, your solver can be chosen with the procedural interface via
893: $     KSPSetType(ksp,"my_solver")
894:    or at runtime via the option
895: $     -ksp_type my_solver

897:    Level: advanced

899: .keywords: KSP, register

901: .seealso: KSPRegisterAll(), KSPRegisterDestroy()

903: @*/
904: PetscErrorCode  KSPRegister(const char sname[],PetscErrorCode (*function)(KSP))
905: {

909:   PetscFunctionListAdd(&KSPList,sname,function);
910:   return(0);
911: }

915: /*@
916:   KSPSetNullSpace - Sets the null space of the operator

918:   Logically Collective on KSP

920:   Input Parameters:
921: +  ksp - the Krylov space object
922: -  nullsp - the null space of the operator

924:   Notes: If the Mat provided to KSP has a nullspace added to it with MatSetNullSpace() then
925:          KSP will automatically use the MatNullSpace and you don't need to call KSPSetNullSpace().

927:   Level: advanced

929: .seealso: KSPSetOperators(), MatNullSpaceCreate(), KSPGetNullSpace(), MatSetNullSpace()
930: @*/
931: PetscErrorCode  KSPSetNullSpace(KSP ksp,MatNullSpace nullsp)
932: {

938:   PetscObjectReference((PetscObject)nullsp);
939:   if (ksp->nullsp) { MatNullSpaceDestroy(&ksp->nullsp); }
940:   ksp->nullsp = nullsp;
941:   return(0);
942: }

946: /*@
947:   KSPGetNullSpace - Gets the null space of the operator

949:   Not Collective

951:   Input Parameters:
952: +  ksp - the Krylov space object
953: -  nullsp - the null space of the operator

955:   Level: advanced

957: .seealso: KSPSetOperators(), MatNullSpaceCreate(), KSPSetNullSpace()
958: @*/
959: PetscErrorCode  KSPGetNullSpace(KSP ksp,MatNullSpace *nullsp)
960: {
964:   *nullsp = ksp->nullsp;
965:   return(0);
966: }