Actual source code: itcreate.c

  1: #define PETSCKSP_DLL

  3: /*
  4:      The basic KSP routines, Create, View etc. are here.
  5: */
 6:  #include include/private/kspimpl.h
 7:  #include petscsys.h

  9: /* Logging support */
 10: PetscCookie  KSP_COOKIE = 0;
 11: PetscEvent  KSP_GMRESOrthogonalization = 0, KSP_SetUp = 0, KSP_Solve = 0;


 14: PetscTruth KSPRegisterAllCalled = PETSC_FALSE;

 18: /*@C 
 19:    KSPView - Prints the KSP data structure.

 21:    Collective on KSP

 23:    Input Parameters:
 24: +  ksp - the Krylov space context
 25: -  viewer - visualization context

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

 30:    Note:
 31:    The available visualization contexts include
 32: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 33: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 34:          output where only the first processor opens
 35:          the file.  All other processors send their 
 36:          data to the first processor to print. 

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

 41:    Level: beginner

 43: .keywords: KSP, view

 45: .seealso: PCView(), PetscViewerASCIIOpen()
 46: @*/
 47: PetscErrorCode  KSPView(KSP ksp,PetscViewer viewer)
 48: {
 49:   const char     *type;
 51:   PetscTruth     iascii;

 55:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(ksp->comm);

 59:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
 60:   if (iascii) {
 61:     KSPGetType(ksp,&type);
 62:     if (ksp->prefix) {
 63:       PetscViewerASCIIPrintf(viewer,"KSP Object:(%s)\n",ksp->prefix);
 64:     } else {
 65:       PetscViewerASCIIPrintf(viewer,"KSP Object:\n");
 66:     }
 67:     if (type) {
 68:       PetscViewerASCIIPrintf(viewer,"  type: %s\n",type);
 69:     } else {
 70:       PetscViewerASCIIPrintf(viewer,"  type: not yet set\n");
 71:     }
 72:     if (ksp->ops->view) {
 73:       PetscViewerASCIIPushTab(viewer);
 74:       (*ksp->ops->view)(ksp,viewer);
 75:       PetscViewerASCIIPopTab(viewer);
 76:     }
 77:     if (ksp->guess_zero) {PetscViewerASCIIPrintf(viewer,"  maximum iterations=%D, initial guess is zero\n",ksp->max_it);}
 78:     else                 {PetscViewerASCIIPrintf(viewer,"  maximum iterations=%D\n", ksp->max_it);}
 79:     if (ksp->guess_knoll) {PetscViewerASCIIPrintf(viewer,"  using preconditioner applied to right hand side for initial guess\n");}
 80:     PetscViewerASCIIPrintf(viewer,"  tolerances:  relative=%G, absolute=%G, divergence=%G\n",ksp->rtol,ksp->abstol,ksp->divtol);
 81:     if (ksp->pc_side == PC_RIGHT)          {PetscViewerASCIIPrintf(viewer,"  right preconditioning\n");}
 82:     else if (ksp->pc_side == PC_SYMMETRIC) {PetscViewerASCIIPrintf(viewer,"  symmetric preconditioning\n");}
 83:     else                                   {PetscViewerASCIIPrintf(viewer,"  left preconditioning\n");}
 84:   } else {
 85:     if (ksp->ops->view) {
 86:       (*ksp->ops->view)(ksp,viewer);
 87:     }
 88:   }
 89:   PCView(ksp->pc,viewer);
 90:   return(0);
 91: }

 93: /*
 94:    Contains the list of registered KSP routines
 95: */
 96: PetscFList KSPList = 0;

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

103:    Collective on KSP

105:    Input Parameter:
106: +  ksp - Krylov solver context
107: -  normtype - one of 
108: $   KSP_NO_NORM - skips computing the norm, this should only be used if you are using
109: $                 the Krylov method as a smoother with a fixed small number of iterations.
110: $                 You must also call KSPSetConvergenceTest(ksp,KSPSkipConverged,PETSC_NULL);
111: $                 supported only by CG, Richardson, Bi-CG-stab, CR, and CGS methods.
112: $   KSP_PRECONDITIONED_NORM - the default for left preconditioned solves, uses the l2 norm
113: $                 of the preconditioned residual
114: $   KSP_UNPRECONDITIONED_NORM - uses the l2 norm of the true b - Ax residual, supported only by
115: $                 CG, CHEBYCHEV, and RICHARDSON, automatically true for right (see KSPSetPreconditioningSide) 
116: $                 preconditioning..
117: $   KSP_NATURAL_NORM - supported  by cg, cr, and cgs 


120:    Options Database Key:
121: .   -ksp_norm_type <none,preconditioned,unpreconditioned,natural>

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

126:    Level: advanced

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

130: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged()                               
131: @*/
132: PetscErrorCode  KSPSetNormType(KSP ksp,KSPNormType normtype)
133: {

138:   ksp->normtype = normtype;
139:   if (normtype == KSP_NO_NORM) {
140:     PetscInfo(ksp,"Warning seting KSPNormType to skip computing the norm\n\
141:   make sure you set the KSP convergence test to KSPSkipConvergence\n");
142:   }
143:   return(0);
144: }

148: static PetscErrorCode KSPPublish_Petsc(PetscObject obj)
149: {
151:   return(0);
152: }

156: /*@
157:    KSPSetOperators - Sets the matrix associated with the linear system
158:    and a (possibly) different one associated with the preconditioner. 

160:    Collective on KSP and Mat

162:    Input Parameters:
163: +  ksp - the KSP context
164: .  Amat - the matrix associated with the linear system
165: .  Pmat - the matrix to be used in constructing the preconditioner, usually the
166:           same as Amat. 
167: -  flag - flag indicating information about the preconditioner matrix structure
168:    during successive linear solves.  This flag is ignored the first time a
169:    linear system is solved, and thus is irrelevant when solving just one linear
170:    system.

172:    Notes: 
173:    The flag can be used to eliminate unnecessary work in the preconditioner 
174:    during the repeated solution of linear systems of the same size.  The
175:    available options are
176: $    SAME_PRECONDITIONER -
177: $      Pmat is identical during successive linear solves.
178: $      This option is intended for folks who are using
179: $      different Amat and Pmat matrices and want to reuse the
180: $      same preconditioner matrix.  For example, this option
181: $      saves work by not recomputing incomplete factorization
182: $      for ILU/ICC preconditioners.
183: $    SAME_NONZERO_PATTERN -
184: $      Pmat has the same nonzero structure during
185: $      successive linear solves. 
186: $    DIFFERENT_NONZERO_PATTERN -
187: $      Pmat does not have the same nonzero structure.

189:     Caution:
190:     If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion
191:     and does not check the structure of the matrix.  If you erroneously
192:     claim that the structure is the same when it actually is not, the new
193:     preconditioner will not function correctly.  Thus, use this optimization
194:     feature carefully!

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

199:     Level: beginner

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

209: $         KSP/PCGetOperators(ksp/pc,&mat,PETSC_NULL,PETSC_NULL); is equivalent to
210: $           set size, type, etc of mat

212: $         MatCreate(comm,&mat);
213: $         KSP/PCSetOperators(ksp/pc,mat,mat,SAME_NONZERO_PATTERN);
214: $         PetscObjectDereference((PetscObject)mat);
215: $           set size, type, etc of mat

217:      and

219: $         KSP/PCGetOperators(ksp/pc,&mat,&pmat,PETSC_NULL); is equivalent to
220: $           set size, type, etc of mat and pmat

222: $         MatCreate(comm,&mat);
223: $         MatCreate(comm,&pmat);
224: $         KSP/PCSetOperators(ksp/pc,mat,pmat,SAME_NONZERO_PATTERN);
225: $         PetscObjectDereference((PetscObject)mat);
226: $         PetscObjectDereference((PetscObject)pmat);
227: $           set size, type, etc of mat and pmat

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

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

240: .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPGetOperators()
241: @*/
242: PetscErrorCode  KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat,MatStructure flag)
243: {

252:   PCSetOperators(ksp->pc,Amat,Pmat,flag);
253:   if (ksp->setupcalled > 1) ksp->setupcalled = 1;  /* so that next solve call will call setup */
254:   return(0);
255: }

259: /*@
260:    KSPGetOperators - Gets the matrix associated with the linear system
261:    and a (possibly) different one associated with the preconditioner. 

263:    Collective on KSP and Mat

265:    Input Parameter:
266: .  ksp - the KSP context

268:    Output Parameters:
269: +  Amat - the matrix associated with the linear system
270: .  Pmat - the matrix to be used in constructing the preconditioner, usually the
271:           same as Amat. 
272: -  flag - flag indicating information about the preconditioner matrix structure
273:    during successive linear solves.  This flag is ignored the first time a
274:    linear system is solved, and thus is irrelevant when solving just one linear
275:    system.

277:     Level: intermediate

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

281: .seealso: KSPSolve(), KSPGetPC(), PCGetOperators(), PCSetOperators(), KSPSetOperators(), KSPGetOperatorsSet()
282: @*/
283: PetscErrorCode  KSPGetOperators(KSP ksp,Mat *Amat,Mat *Pmat,MatStructure *flag)
284: {

289:   PCGetOperators(ksp->pc,Amat,Pmat,flag);
290:   return(0);
291: }

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

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

301:    Input Parameter:
302: .  pc - the preconditioner context

304:    Output Parameters:
305: +  mat - the matrix associated with the linear system was set
306: -  pmat - matrix associated with the preconditioner was set, usually the same

308:    Level: intermediate

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

312: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators(), PCGetOperatorsSet()
313: @*/
314: PetscErrorCode  KSPGetOperatorsSet(KSP ksp,PetscTruth *mat,PetscTruth *pmat)
315: {

320:   PCGetOperatorsSet(ksp->pc,mat,pmat);
321:   return(0);
322: }

326: /*@
327:    KSPCreate - Creates the default KSP context.

329:    Collective on MPI_Comm

331:    Input Parameter:
332: .  comm - MPI communicator

334:    Output Parameter:
335: .  ksp - location to put the KSP context

337:    Notes:
338:    The default KSP type is GMRES with a restart of 30, using modified Gram-Schmidt
339:    orthogonalization.

341:    Level: beginner

343: .keywords: KSP, create, context

345: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP
346: @*/
347: PetscErrorCode  KSPCreate(MPI_Comm comm,KSP *inksp)
348: {
349:   KSP            ksp;

354:   *inksp = 0;
355: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
356:   KSPInitializePackage(PETSC_NULL);
357: #endif

359:   PetscHeaderCreate(ksp,_p_KSP,struct _KSPOps,KSP_COOKIE,-1,"KSP",comm,KSPDestroy,KSPView);
360:   *inksp             = ksp;
361:   ksp->bops->publish = KSPPublish_Petsc;

363:   ksp->type          = -1;
364:   ksp->max_it        = 10000;
365:   ksp->pc_side       = PC_LEFT;
366:   ksp->rtol          = 1.e-5;
367:   ksp->abstol          = 1.e-50;
368:   ksp->divtol        = 1.e4;

370:   ksp->normtype            = KSP_PRECONDITIONED_NORM;
371:   ksp->rnorm               = 0.0;
372:   ksp->its                 = 0;
373:   ksp->guess_zero          = PETSC_TRUE;
374:   ksp->calc_sings          = PETSC_FALSE;
375:   ksp->res_hist            = PETSC_NULL;
376:   ksp->res_hist_len        = 0;
377:   ksp->res_hist_max        = 0;
378:   ksp->res_hist_reset      = PETSC_TRUE;
379:   ksp->numbermonitors      = 0;
380:   ksp->converged           = KSPDefaultConverged;
381:   ksp->ops->buildsolution  = KSPDefaultBuildSolution;
382:   ksp->ops->buildresidual  = KSPDefaultBuildResidual;

384:   ksp->ops->setfromoptions = 0;

386:   ksp->vec_sol         = 0;
387:   ksp->vec_rhs         = 0;
388:   ksp->pc              = 0;

390:   ksp->ops->solve      = 0;
391:   ksp->ops->setup      = 0;
392:   ksp->ops->destroy    = 0;

394:   ksp->data            = 0;
395:   ksp->nwork           = 0;
396:   ksp->work            = 0;

398:   ksp->cnvP            = 0;

400:   ksp->reason          = KSP_CONVERGED_ITERATING;

402:   ksp->setupcalled     = 0;
403:   PetscPublishAll(ksp);
404:   PCCreate(comm,&ksp->pc);
405:   return(0);
406: }
407: 
410: /*@C
411:    KSPSetType - Builds KSP for a particular solver. 

413:    Collective on KSP

415:    Input Parameters:
416: +  ksp      - the Krylov space context
417: -  type - a known method

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

423:    Notes:  
424:    See "petsc/include/petscksp.h" for available methods (for instance,
425:    KSPCG or KSPGMRES).

427:   Normally, it is best to use the KSPSetFromOptions() command and
428:   then set the KSP type from the options database rather than by using
429:   this routine.  Using the options database provides the user with
430:   maximum flexibility in evaluating the many different Krylov methods.
431:   The KSPSetType() routine is provided for those situations where it
432:   is necessary to set the iterative solver independently of the command
433:   line or options database.  This might be the case, for example, when
434:   the choice of iterative solver changes during the execution of the
435:   program, and the user's application is taking responsibility for
436:   choosing the appropriate method.  In other words, this routine is
437:   not for beginners.

439:   Level: intermediate

441: .keywords: KSP, set, method

443: .seealso: PCSetType(), KSPType

445: @*/
446: PetscErrorCode  KSPSetType(KSP ksp, KSPType type)
447: {
448:   PetscErrorCode ierr,(*r)(KSP);
449:   PetscTruth     match;


455:   PetscTypeCompare((PetscObject)ksp,type,&match);
456:   if (match) return(0);

458:   if (ksp->data) {
459:     /* destroy the old private KSP context */
460:     (*ksp->ops->destroy)(ksp);
461:     ksp->data = 0;
462:   }
463:   /* Get the function pointers for the iterative method requested */
464:   if (!KSPRegisterAllCalled) {KSPRegisterAll(PETSC_NULL);}
465:    PetscFListFind(ksp->comm,KSPList,type,(void (**)(void)) &r);
466:   if (!r) SETERRQ1(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown KSP type given: %s",type);
467:   ksp->setupcalled = 0;
468:   (*r)(ksp);
469:   PetscObjectChangeTypeName((PetscObject)ksp,type);
470:   return(0);
471: }

475: /*@
476:    KSPRegisterDestroy - Frees the list of KSP methods that were
477:    registered by KSPRegisterDynamic().

479:    Not Collective

481:    Level: advanced

483: .keywords: KSP, register, destroy

485: .seealso: KSPRegisterDynamic(), KSPRegisterAll()
486: @*/
487: PetscErrorCode  KSPRegisterDestroy(void)
488: {

492:   if (KSPList) {
493:     PetscFListDestroy(&KSPList);
494:     KSPList = 0;
495:   }
496:   KSPRegisterAllCalled = PETSC_FALSE;
497:   return(0);
498: }

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

505:    Not Collective

507:    Input Parameter:
508: .  ksp - Krylov context 

510:    Output Parameter:
511: .  name - name of KSP method 

513:    Level: intermediate

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

517: .seealso: KSPSetType()
518: @*/
519: PetscErrorCode  KSPGetType(KSP ksp,KSPType *type)
520: {
524:   *type = ksp->type_name;
525:   return(0);
526: }

530: /*@C
531:   KSPRegister - See KSPRegisterDynamic()

533:   Level: advanced
534: @*/
535: PetscErrorCode  KSPRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(KSP))
536: {
538:   char           fullname[PETSC_MAX_PATH_LEN];

541:   PetscFListConcat(path,name,fullname);
542:   PetscFListAdd(&KSPList,sname,fullname,(void (*)(void))function);
543:   return(0);
544: }

548: /*@
549:   KSPSetNullSpace - Sets the null space of the operator

551:   Collective on KSP

553:   Input Parameters:
554: +  ksp - the Krylov space object
555: -  nullsp - the null space of the operator

557:   Level: advanced

559: .seealso: KSPSetOperators(), MatNullSpaceCreate(), KSPGetNullSpace()
560: @*/
561: PetscErrorCode  KSPSetNullSpace(KSP ksp,MatNullSpace nullsp)
562: {

566:   PetscObjectReference((PetscObject)nullsp);
567:   if (ksp->nullsp) { MatNullSpaceDestroy(ksp->nullsp); }
568:   ksp->nullsp = nullsp;
569:   return(0);
570: }

574: /*@
575:   KSPGetNullSpace - Gets the null space of the operator

577:   Collective on KSP

579:   Input Parameters:
580: +  ksp - the Krylov space object
581: -  nullsp - the null space of the operator

583:   Level: advanced

585: .seealso: KSPSetOperators(), MatNullSpaceCreate(), KSPSetNullSpace()
586: @*/
587: PetscErrorCode  KSPGetNullSpace(KSP ksp,MatNullSpace *nullsp)
588: {
590:   *nullsp = ksp->nullsp;
591:   return(0);
592: }