Actual source code: itcreate.c

  1: /*$Id: itcreate.c,v 1.206 2001/08/06 21:16:38 bsmith Exp $*/
  2: /*
  3:      The basic KSP routines, Create, View etc. are here.
  4: */
 5:  #include src/sles/ksp/kspimpl.h
 6:  #include petscsys.h

  8: /* Logging support */
  9: int KSP_COOKIE;
 10: int KSP_GMRESOrthogonalization;

 12: EXTERN int SLESInitializePackage(char *);

 14: PetscTruth KSPRegisterAllCalled = PETSC_FALSE;

 16: #undef __FUNCT__  
 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:    Note:
 28:    The available visualization contexts include
 29: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
 30: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
 31:          output where only the first processor opens
 32:          the file.  All other processors send their 
 33:          data to the first processor to print. 

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

 38:    Level: developer

 40: .keywords: KSP, view

 42: .seealso: PCView(), PetscViewerASCIIOpen()
 43: @*/
 44: int KSPView(KSP ksp,PetscViewer viewer)
 45: {
 46:   char        *type;
 47:   int         ierr;
 48:   PetscTruth  isascii;

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

 56:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
 57:   if (isascii) {
 58:     KSPGetType(ksp,&type);
 59:     if (ksp->prefix) {
 60:       PetscViewerASCIIPrintf(viewer,"KSP Object:(%s)n",ksp->prefix);
 61:     } else {
 62:       PetscViewerASCIIPrintf(viewer,"KSP Object:n");
 63:     }
 64:     if (type) {
 65:       PetscViewerASCIIPrintf(viewer,"  type: %sn",type);
 66:     } else {
 67:       PetscViewerASCIIPrintf(viewer,"  type: not yet setn");
 68:     }
 69:     if (ksp->ops->view) {
 70:       PetscViewerASCIIPushTab(viewer);
 71:       (*ksp->ops->view)(ksp,viewer);
 72:       PetscViewerASCIIPopTab(viewer);
 73:     }
 74:     if (ksp->guess_zero) {PetscViewerASCIIPrintf(viewer,"  maximum iterations=%d, initial guess is zeron",ksp->max_it);}
 75:     else                 {PetscViewerASCIIPrintf(viewer,"  maximum iterations=%dn", ksp->max_it);}
 76:     if (ksp->guess_knoll) {PetscViewerASCIIPrintf(viewer,"  using preconditioner applied to right hand side for initial guessn");}
 77:     PetscViewerASCIIPrintf(viewer,"  tolerances:  relative=%g, absolute=%g, divergence=%gn",ksp->rtol,ksp->atol,ksp->divtol);
 78:     if (ksp->pc_side == PC_RIGHT)          {PetscViewerASCIIPrintf(viewer,"  right preconditioningn");}
 79:     else if (ksp->pc_side == PC_SYMMETRIC) {PetscViewerASCIIPrintf(viewer,"  symmetric preconditioningn");}
 80:     else                                   {PetscViewerASCIIPrintf(viewer,"  left preconditioningn");}
 81:   } else {
 82:     if (ksp->ops->view) {
 83:       (*ksp->ops->view)(ksp,viewer);
 84:     }
 85:   }
 86:   return(0);
 87: }

 89: /*
 90:    Contains the list of registered KSP routines
 91: */
 92: PetscFList KSPList = 0;

 94: #undef __FUNCT__  
 96: /*@C
 97:    KSPSetNormType - Sets the norm that is used for convergence testing.

 99:    Collective on KSP

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


115:    Options Database Key:
116: .   -ksp_norm_type <none,preconditioned,unpreconditioned,natural>

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

121:    Level: advanced

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

125: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSPSkipConverged()                               
126: @*/
127: int KSPSetNormType(KSP ksp,KSPNormType normtype)
128: {

132:   ksp->normtype = normtype;
133:   if (normtype == KSP_NO_NORM) {
134:     PetscLogInfo(ksp,"KSPSetNormType:Warning seting KSPNormType to skip computing the normn
135:   make sure you set the KSP convergence test to KSPSkipConvergencen");
136:   }
137:   return(0);
138: }

140: #undef __FUNCT__  
142: static int KSPPublish_Petsc(PetscObject obj)
143: {
144: #if defined(PETSC_HAVE_AMS)
145:   KSP          v = (KSP) obj;
146:   int          ierr;
147: #endif


151: #if defined(PETSC_HAVE_AMS)
152:   /* if it is already published then return */
153:   if (v->amem >=0) return(0);

155:   PetscObjectPublishBaseBegin(obj);
156:   AMS_Memory_add_field((AMS_Memory)v->amem,"Iteration",&v->its,1,AMS_INT,AMS_READ,
157:                                 AMS_COMMON,AMS_REDUCT_UNDEF);
158:   AMS_Memory_add_field((AMS_Memory)v->amem,"Residual",&v->rnorm,1,AMS_DOUBLE,AMS_READ,
159:                                 AMS_COMMON,AMS_REDUCT_UNDEF);

161:   if (v->res_hist_max > 0) {
162:     AMS_Memory_add_field((AMS_Memory)v->amem,"ResidualNormsCount",&v->res_hist_len,1,AMS_INT,AMS_READ,
163:                                 AMS_COMMON,AMS_REDUCT_UNDEF);
164:     AMS_Memory_add_field((AMS_Memory)v->amem,"ResidualNormsCountMax",&v->res_hist_max,1,AMS_INT,AMS_READ,
165:                                 AMS_COMMON,AMS_REDUCT_UNDEF);
166:     AMS_Memory_add_field((AMS_Memory)v->amem,"ResidualNorms",v->res_hist,v->res_hist_max,AMS_DOUBLE,AMS_READ,
167:                                 AMS_COMMON,AMS_REDUCT_UNDEF);
168:   }

170:   PetscObjectPublishBaseEnd(obj);
171: #endif

173:   return(0);
174: }


177: #undef __FUNCT__  
179: /*@C
180:    KSPCreate - Creates the default KSP context.

182:    Collective on MPI_Comm

184:    Input Parameter:
185: .  comm - MPI communicator

187:    Output Parameter:
188: .  ksp - location to put the KSP context

190:    Notes:
191:    The default KSP type is GMRES with a restart of 30, using modified Gram-Schmidt
192:    orthogonalization.

194:    Level: developer

196: .keywords: KSP, create, context

198: .seealso: KSPSetUp(), KSPSolve(), KSPDestroy(), KSP
199: @*/
200: int KSPCreate(MPI_Comm comm,KSP *inksp)
201: {
202:   KSP ksp;

207:   *inksp = 0;
208: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
209:   SLESInitializePackage(PETSC_NULL);
210: #endif

212:   PetscHeaderCreate(ksp,_p_KSP,struct _KSPOps,KSP_COOKIE,-1,"KSP",comm,KSPDestroy,KSPView);
213:   PetscLogObjectCreate(ksp);
214:   *inksp             = ksp;
215:   ksp->bops->publish = KSPPublish_Petsc;

217:   ksp->type          = -1;
218:   ksp->max_it        = 10000;
219:   ksp->pc_side       = PC_LEFT;
220:   ksp->rtol          = 1.e-5;
221:   ksp->atol          = 1.e-50;
222:   ksp->divtol        = 1.e4;

224:   ksp->normtype            = KSP_PRECONDITIONED_NORM;
225:   ksp->rnorm               = 0.0;
226:   ksp->its                 = 0;
227:   ksp->guess_zero          = PETSC_TRUE;
228:   ksp->calc_sings          = PETSC_FALSE;
229:   ksp->res_hist            = PETSC_NULL;
230:   ksp->res_hist_len        = 0;
231:   ksp->res_hist_max        = 0;
232:   ksp->res_hist_reset      = PETSC_TRUE;
233:   ksp->numbermonitors      = 0;
234:   ksp->converged           = KSPDefaultConverged;
235:   ksp->ops->buildsolution  = KSPDefaultBuildSolution;
236:   ksp->ops->buildresidual  = KSPDefaultBuildResidual;

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

240:   ksp->vec_sol         = 0;
241:   ksp->vec_rhs         = 0;
242:   ksp->B               = 0;

244:   ksp->ops->solve      = 0;
245:   ksp->ops->setup      = 0;
246:   ksp->ops->destroy    = 0;

248:   ksp->data            = 0;
249:   ksp->nwork           = 0;
250:   ksp->work            = 0;

252:   ksp->cnvP            = 0;

254:   ksp->reason          = KSP_CONVERGED_ITERATING;

256:   ksp->setupcalled     = 0;
257:   PetscPublishAll(ksp);
258:   return(0);
259: }
260: 
261: #undef __FUNCT__  
263: /*@C
264:    KSPSetType - Builds KSP for a particular solver. 

266:    Collective on KSP

268:    Input Parameters:
269: +  ksp      - the Krylov space context
270: -  type - a known method

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

276:    Notes:  
277:    See "petsc/include/petscksp.h" for available methods (for instance,
278:    KSPCG or KSPGMRES).

280:   Normally, it is best to use the SLESSetFromOptions() command and
281:   then set the KSP type from the options database rather than by using
282:   this routine.  Using the options database provides the user with
283:   maximum flexibility in evaluating the many different Krylov methods.
284:   The KSPSetType() routine is provided for those situations where it
285:   is necessary to set the iterative solver independently of the command
286:   line or options database.  This might be the case, for example, when
287:   the choice of iterative solver changes during the execution of the
288:   program, and the user's application is taking responsibility for
289:   choosing the appropriate method.  In other words, this routine is
290:   not for beginners.

292:   Level: intermediate

294: .keywords: KSP, set, method

296: .seealso: PCSetType(), KSPType

298: @*/
299: int KSPSetType(KSP ksp,KSPType type)
300: {
301:   int        ierr,(*r)(KSP);
302:   PetscTruth match;


308:   PetscTypeCompare((PetscObject)ksp,type,&match);
309:   if (match) return(0);

311:   if (ksp->data) {
312:     /* destroy the old private KSP context */
313:     (*ksp->ops->destroy)(ksp);
314:     ksp->data = 0;
315:   }
316:   /* Get the function pointers for the iterative method requested */
317:   if (!KSPRegisterAllCalled) {KSPRegisterAll(PETSC_NULL);}

319:    PetscFListFind(ksp->comm,KSPList,type,(void (**)(void)) &r);

321:   if (!r) SETERRQ1(1,"Unknown KSP type given: %s",type);

323:   ksp->setupcalled = 0;
324:   (*r)(ksp);

326:   PetscObjectChangeTypeName((PetscObject)ksp,type);
327:   return(0);
328: }

330: #undef __FUNCT__  
332: /*@C
333:    KSPRegisterDestroy - Frees the list of KSP methods that were
334:    registered by KSPRegisterDynamic().

336:    Not Collective

338:    Level: advanced

340: .keywords: KSP, register, destroy

342: .seealso: KSPRegisterDynamic(), KSPRegisterAll()
343: @*/
344: int KSPRegisterDestroy(void)
345: {

349:   if (KSPList) {
350:     PetscFListDestroy(&KSPList);
351:     KSPList = 0;
352:   }
353:   KSPRegisterAllCalled = PETSC_FALSE;
354:   return(0);
355: }

357: #undef __FUNCT__  
359: /*@C
360:    KSPGetType - Gets the KSP type as a string from the KSP object.

362:    Not Collective

364:    Input Parameter:
365: .  ksp - Krylov context 

367:    Output Parameter:
368: .  name - name of KSP method 

370:    Level: intermediate

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

374: .seealso: KSPSetType()
375: @*/
376: int KSPGetType(KSP ksp,KSPType *type)
377: {
380:   *type = ksp->type_name;
381:   return(0);
382: }

384: #undef __FUNCT__  
386: /*@
387:    KSPSetFromOptions - Sets KSP options from the options database.
388:    This routine must be called before KSPSetUp() if the user is to be 
389:    allowed to set the Krylov type. 

391:    Collective on KSP

393:    Input Parameters:
394: .  ksp - the Krylov space context

396:    Options Database Keys:
397: +   -ksp_max_it - maximum number of linear iterations
398: .   -ksp_rtol rtol - relative tolerance used in default determination of convergence, i.e.
399:                 if residual norm decreases by this factor than convergence is declared
400: .   -ksp_atol atol - absolute tolerance used in default convergence test, i.e. if residual 
401:                 norm is less than this then convergence is declared
402: .   -ksp_divtol tol - if residual norm increases by this factor than divergence is declared
403: .   -ksp_norm_type - none - skip norms used in convergence tests (useful only when not using 
404: $                       convergence test (say you always want to run with 5 iterations) to 
405: $                       save on communication overhead
406: $                    preconditioned - default for left preconditioning 
407: $                    unpreconditioned - see KSPSetNormType()
408: $                    natural - see KSPSetNormType()
409: .   -ksp_knoll - compute initial guess by applying the preconditioner to the right hand side
410: .   -ksp_cancelmonitors - cancel all previous convergene monitor routines set
411: .   -ksp_monitor - print residual norm at each iteration
412: .   -ksp_xmonitor - plot residual norm at each iteration
413: .   -ksp_vecmonitor - plot solution at each iteration
414: -   -ksp_singmonitor - monitor extremem singular values at each iteration

416:    Notes:  
417:    To see all options, run your program with the -help option
418:    or consult the users manual.

420:    Level: developer

422: .keywords: KSP, set, from, options, database

424: .seealso: 
425: @*/
426: int KSPSetFromOptions(KSP ksp)
427: {
428:   int        ierr;
429:   char       type[256],*stype[] = {"none","preconditioned","unpreconditioned","natural"};
430:   PetscTruth flg;

434:   if (!KSPRegisterAllCalled) {KSPRegisterAll(PETSC_NULL);}
435:   PetscOptionsBegin(ksp->comm,ksp->prefix,"Krylov Method (KSP) Options","KSP");
436:     PetscOptionsList("-ksp_type","Krylov method","KSPSetType",KSPList,(char*)(ksp->type_name?ksp->type_name:KSPGMRES),type,256,&flg);
437:     if (flg) {
438:       KSPSetType(ksp,type);
439:     }
440:     /*
441:       Set the type if it was never set.
442:     */
443:     if (!ksp->type_name) {
444:       KSPSetType(ksp,KSPGMRES);
445:     }

447:     PetscOptionsInt("-ksp_max_it","Maximum number of iterations","KSPSetTolerances",ksp->max_it,&ksp->max_it,PETSC_NULL);
448:     PetscOptionsReal("-ksp_rtol","Relative decrease in residual norm","KSPSetTolerances",ksp->rtol,&ksp->rtol,PETSC_NULL);
449:     PetscOptionsReal("-ksp_atol","Absolute value of residual norm","KSPSetTolerances",ksp->atol,&ksp->atol,PETSC_NULL);
450:     PetscOptionsReal("-ksp_divtol","Residual norm increase cause divergence","KSPSetTolerances",ksp->divtol,&ksp->divtol,PETSC_NULL);
451:     PetscOptionsLogical("-ksp_knoll","Use preconditioner applied to b for initial guess","KSPSetInitialGuessKnoll",ksp->guess_knoll,
452:                                   &ksp->guess_knoll,PETSC_NULL);

454:     PetscOptionsEList("-ksp_norm_type","KSP Norm type","KSPSetNormType",stype,4,"preconditioned",type,256,&flg);
455:     if (flg) {
456:       PetscTruth isnone,ispreconditioned,isunpreconditioned,isnatural;

458:       PetscStrcmp(type,stype[0],&isnone);
459:       PetscStrcmp(type,stype[1],&ispreconditioned);
460:       PetscStrcmp(type,stype[2],&isunpreconditioned);
461:       PetscStrcmp(type,stype[3],&isnatural);

463:       if (isnone) {
464:         KSPSetNormType(ksp,KSP_NO_NORM);
465:         KSPSetConvergenceTest(ksp,KSPSkipConverged,0);
466:       } else if (ispreconditioned) {
467:         KSPSetNormType(ksp,KSP_PRECONDITIONED_NORM);
468:       } else if (isunpreconditioned) {
469:         KSPSetNormType(ksp,KSP_UNPRECONDITIONED_NORM);
470:       } else if (isnatural) {
471:         KSPSetNormType(ksp,KSP_NATURAL_NORM);
472:       } else {
473:         SETERRQ1(1,"Unknown KSP normtype %s",type);
474:       }
475:     }

477:     PetscOptionsName("-ksp_cancelmonitors","Remove any hardwired monitor routines","KSPClearMonitor",&flg);
478:     /* -----------------------------------------------------------------------*/
479:     /*
480:       Cancels all monitors hardwired into code before call to KSPSetFromOptions()
481:     */
482:     if (flg) {
483:       KSPClearMonitor(ksp);
484:     }
485:     /*
486:       Prints preconditioned residual norm at each iteration
487:     */
488:     PetscOptionsName("-ksp_monitor","Monitor preconditioned residual norm","KSPSetMonitor",&flg);
489:     if (flg) {
490:       KSPSetMonitor(ksp,KSPDefaultMonitor,PETSC_NULL,PETSC_NULL);
491:     }
492:     /*
493:       Plots the vector solution 
494:     */
495:     PetscOptionsName("-ksp_vecmonitor","Monitor solution graphically","KSPSetMonitor",&flg);
496:     if (flg) {
497:       KSPSetMonitor(ksp,KSPVecViewMonitor,PETSC_NULL,PETSC_NULL);
498:     }
499:     /*
500:       Prints preconditioned and true residual norm at each iteration
501:     */
502:     PetscOptionsName("-ksp_truemonitor","Monitor true (unpreconditioned) residual norm","KSPSetMonitor",&flg);
503:     if (flg) {
504:       KSPSetMonitor(ksp,KSPTrueMonitor,PETSC_NULL,PETSC_NULL);
505:     }
506:     /*
507:       Prints extreme eigenvalue estimates at each iteration
508:     */
509:     PetscOptionsName("-ksp_singmonitor","Monitor singular values","KSPSetMonitor",&flg);
510:     if (flg) {
511:       KSPSetComputeSingularValues(ksp,PETSC_TRUE);
512:       KSPSetMonitor(ksp,KSPSingularValueMonitor,PETSC_NULL,PETSC_NULL);
513:     }
514:     /*
515:       Prints preconditioned residual norm with fewer digits
516:     */
517:     PetscOptionsName("-ksp_smonitor","Monitor preconditioned residual norm with fewer digitis","KSPSetMonitor",&flg);
518:     if (flg) {
519:       KSPSetMonitor(ksp,KSPDefaultSMonitor,PETSC_NULL,PETSC_NULL);
520:     }
521:     /*
522:       Graphically plots preconditioned residual norm
523:     */
524:     PetscOptionsName("-ksp_xmonitor","Monitor graphically preconditioned residual norm","KSPSetMonitor",&flg);
525:     if (flg) {
526:       KSPSetMonitor(ksp,KSPLGMonitor,PETSC_NULL,PETSC_NULL);
527:     }
528:     /*
529:       Graphically plots preconditioned and true residual norm
530:     */
531:     PetscOptionsName("-ksp_xtruemonitor","Monitor graphically true residual norm","KSPSetMonitor",&flg);
532:     if (flg){
533:       KSPSetMonitor(ksp,KSPLGTrueMonitor,PETSC_NULL,PETSC_NULL);
534:     }

536:     /* -----------------------------------------------------------------------*/

538:     PetscOptionsLogicalGroupBegin("-ksp_left_pc","Use left preconditioning","KSPSetPreconditionerSide",&flg);
539:     if (flg) { KSPSetPreconditionerSide(ksp,PC_LEFT); }
540:     PetscOptionsLogicalGroup("-ksp_right_pc","Use right preconditioning","KSPSetPreconditionerSide",&flg);
541:     if (flg) { KSPSetPreconditionerSide(ksp,PC_RIGHT);}
542:     PetscOptionsLogicalGroupEnd("-ksp_symmetric_pc","Use symmetric (factorized) preconditioning","KSPSetPreconditionerSide",&flg);
543:     if (flg) { KSPSetPreconditionerSide(ksp,PC_SYMMETRIC);}

545:     PetscOptionsName("-ksp_compute_singularvalues","Compute singular values of preconditioned operator","KSPSetComputeSingularValues",&flg);
546:     if (flg) { KSPSetComputeSingularValues(ksp,PETSC_TRUE); }
547:     PetscOptionsName("-ksp_compute_eigenvalues","Compute eigenvalues of preconditioned operator","KSPSetComputeSingularValues",&flg);
548:     if (flg) { KSPSetComputeSingularValues(ksp,PETSC_TRUE); }
549:     PetscOptionsName("-ksp_plot_eigenvalues","Scatter plot extreme eigenvalues","KSPSetComputeSingularValues",&flg);
550:     if (flg) { KSPSetComputeSingularValues(ksp,PETSC_TRUE); }

552:     if (ksp->ops->setfromoptions) {
553:       (*ksp->ops->setfromoptions)(ksp);
554:     }
555:   PetscOptionsEnd();


558:   return(0);
559: }

561: /*MC
562:    KSPRegisterDynamic - Adds a method to the Krylov subspace solver package.

564:    Synopsis:
565:    int KSPRegisterDynamic(char *name_solver,char *path,char *name_create,int (*routine_create)(KSP))

567:    Not Collective

569:    Input Parameters:
570: +  name_solver - name of a new user-defined solver
571: .  path - path (either absolute or relative) the library containing this solver
572: .  name_create - name of routine to create method context
573: -  routine_create - routine to create method context

575:    Notes:
576:    KSPRegisterDynamic() may be called multiple times to add several user-defined solvers.

578:    If dynamic libraries are used, then the fourth input argument (routine_create)
579:    is ignored.

581:    Sample usage:
582: .vb
583:    KSPRegisterDynamic("my_solver",/home/username/my_lib/lib/libO/solaris/mylib.a,
584:                "MySolverCreate",MySolverCreate);
585: .ve

587:    Then, your solver can be chosen with the procedural interface via
588: $     KSPSetType(ksp,"my_solver")
589:    or at runtime via the option
590: $     -ksp_type my_solver

592:    Level: advanced

594:    Notes: Environmental variables such as ${PETSC_ARCH}, ${PETSC_DIR}, ${PETSC_LIB_DIR}, ${BOPT},
595:           and others of the form ${any_environmental_variable} occuring in pathname will be 
596:           replaced with appropriate values.
597:          If your function is not being put into a shared library then use KSPRegister() instead

599: .keywords: KSP, register

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

603: M*/

605: #undef __FUNCT__  
607: /*@C
608:   KSPRegister - See KSPRegisterDynamic()

610:   Level: advanced
611: @*/
612: int KSPRegister(char *sname,char *path,char *name,int (*function)(KSP))
613: {
614:   int  ierr;
615:   char fullname[256];

618:   PetscFListConcat(path,name,fullname);
619:   PetscFListAdd(&KSPList,sname,fullname,(void (*)(void))function);
620:   return(0);
621: }