Actual source code: itfunc.c

petsc-3.6.1 2015-07-22
Report Typos and Errors
  2: /*
  3:       Interface KSP routines that the user calls.
  4: */

  6: #include <petsc/private/kspimpl.h>   /*I "petscksp.h" I*/
  7: #include <petscdm.h>

 11: /*@
 12:    KSPComputeExtremeSingularValues - Computes the extreme singular values
 13:    for the preconditioned operator. Called after or during KSPSolve().

 15:    Not Collective

 17:    Input Parameter:
 18: .  ksp - iterative context obtained from KSPCreate()

 20:    Output Parameters:
 21: .  emin, emax - extreme singular values

 23:    Options Database Keys:
 24: .  -ksp_compute_singularvalues - compute extreme singular values and print when KSPSolve completes.

 26:    Notes:
 27:    One must call KSPSetComputeSingularValues() before calling KSPSetUp()
 28:    (or use the option -ksp_compute_eigenvalues) in order for this routine to work correctly.

 30:    Many users may just want to use the monitoring routine
 31:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
 32:    to print the extreme singular values at each iteration of the linear solve.

 34:    Estimates of the smallest singular value may be very inaccurate, especially if the Krylov method has not converged.
 35:    The largest singular value is usually accurate to within a few percent if the method has converged, but is still not
 36:    intended for eigenanalysis.

 38:    Disable restarts if using KSPGMRES, otherwise this estimate will only be using those iterations after the last
 39:    restart. See KSPGMRESSetRestart() for more details.

 41:    Level: advanced

 43: .keywords: KSP, compute, extreme, singular, values

 45: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeEigenvalues()
 46: @*/
 47: PetscErrorCode  KSPComputeExtremeSingularValues(KSP ksp,PetscReal *emax,PetscReal *emin)
 48: {

 55:   if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Singular values not requested before KSPSetUp()");

 57:   if (ksp->ops->computeextremesingularvalues) {
 58:     (*ksp->ops->computeextremesingularvalues)(ksp,emax,emin);
 59:   } else {
 60:     *emin = -1.0;
 61:     *emax = -1.0;
 62:   }
 63:   return(0);
 64: }

 68: /*@
 69:    KSPComputeEigenvalues - Computes the extreme eigenvalues for the
 70:    preconditioned operator. Called after or during KSPSolve().

 72:    Not Collective

 74:    Input Parameter:
 75: +  ksp - iterative context obtained from KSPCreate()
 76: -  n - size of arrays r and c. The number of eigenvalues computed (neig) will, in
 77:        general, be less than this.

 79:    Output Parameters:
 80: +  r - real part of computed eigenvalues, provided by user with a dimension of at least n
 81: .  c - complex part of computed eigenvalues, provided by user with a dimension of at least n
 82: -  neig - actual number of eigenvalues computed (will be less than or equal to n)

 84:    Options Database Keys:
 85: +  -ksp_compute_eigenvalues - Prints eigenvalues to stdout
 86: -  -ksp_plot_eigenvalues - Plots eigenvalues in an x-window display

 88:    Notes:
 89:    The number of eigenvalues estimated depends on the size of the Krylov space
 90:    generated during the KSPSolve() ; for example, with
 91:    CG it corresponds to the number of CG iterations, for GMRES it is the number
 92:    of GMRES iterations SINCE the last restart. Any extra space in r[] and c[]
 93:    will be ignored.

 95:    KSPComputeEigenvalues() does not usually provide accurate estimates; it is
 96:    intended only for assistance in understanding the convergence of iterative
 97:    methods, not for eigenanalysis. For accurate computation of eigenvalues we recommend using
 98:    the excellent package SLEPc.

100:    One must call KSPSetComputeEigenvalues() before calling KSPSetUp()
101:    in order for this routine to work correctly.

103:    Many users may just want to use the monitoring routine
104:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
105:    to print the singular values at each iteration of the linear solve.

107:    Level: advanced

109: .keywords: KSP, compute, extreme, singular, values

111: .seealso: KSPSetComputeSingularValues(), KSPMonitorSingularValue(), KSPComputeExtremeSingularValues()
112: @*/
113: PetscErrorCode  KSPComputeEigenvalues(KSP ksp,PetscInt n,PetscReal r[],PetscReal c[],PetscInt *neig)
114: {

121:   if (n<0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Requested < 0 Eigenvalues");
123:   if (!ksp->calc_sings) SETERRQ(PetscObjectComm((PetscObject)ksp),4,"Eigenvalues not requested before KSPSetUp()");

125:   if (ksp->ops->computeeigenvalues) {
126:     (*ksp->ops->computeeigenvalues)(ksp,n,r,c,neig);
127:   } else {
128:     *neig = 0;
129:   }
130:   return(0);
131: }

135: /*@
136:    KSPSetUpOnBlocks - Sets up the preconditioner for each block in
137:    the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
138:    methods.

140:    Collective on KSP

142:    Input Parameter:
143: .  ksp - the KSP context

145:    Notes:
146:    KSPSetUpOnBlocks() is a routine that the user can optinally call for
147:    more precise profiling (via -log_summary) of the setup phase for these
148:    block preconditioners.  If the user does not call KSPSetUpOnBlocks(),
149:    it will automatically be called from within KSPSolve().

151:    Calling KSPSetUpOnBlocks() is the same as calling PCSetUpOnBlocks()
152:    on the PC context within the KSP context.

154:    Level: advanced

156: .keywords: KSP, setup, blocks

158: .seealso: PCSetUpOnBlocks(), KSPSetUp(), PCSetUp()
159: @*/
160: PetscErrorCode  KSPSetUpOnBlocks(KSP ksp)
161: {

166:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
167:   PCSetUpOnBlocks(ksp->pc);
168:   return(0);
169: }

173: /*@
174:    KSPSetReusePreconditioner - reuse the current preconditioner, do not construct a new one even if the operator changes

176:    Collective on KSP

178:    Input Parameters:
179: +  ksp   - iterative context obtained from KSPCreate()
180: -  flag - PETSC_TRUE to reuse the current preconditioner

182:    Level: intermediate

184: .keywords: KSP, setup

186: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner()
187: @*/
188: PetscErrorCode  KSPSetReusePreconditioner(KSP ksp,PetscBool flag)
189: {

194:   PCSetReusePreconditioner(ksp->pc,flag);
195:   return(0);
196: }

200: /*@
201:    KSPSetSkipPCSetFromOptions - prevents KSPSetFromOptions() from call PCSetFromOptions(). This is used if the same PC is shared by more than one KSP so its options are not resetable for each KSP

203:    Collective on KSP

205:    Input Parameters:
206: +  ksp   - iterative context obtained from KSPCreate()
207: -  flag - PETSC_TRUE to skip calling the PCSetFromOptions()

209:    Level: intermediate

211: .keywords: KSP, setup

213: .seealso: KSPCreate(), KSPSolve(), KSPDestroy(), PCSetReusePreconditioner()
214: @*/
215: PetscErrorCode  KSPSetSkipPCSetFromOptions(KSP ksp,PetscBool flag)
216: {
219:   ksp->skippcsetfromoptions = flag;
220:   return(0);
221: }

225: /*@
226:    KSPSetUp - Sets up the internal data structures for the
227:    later use of an iterative solver.

229:    Collective on KSP

231:    Input Parameter:
232: .  ksp   - iterative context obtained from KSPCreate()

234:    Level: developer

236: .keywords: KSP, setup

238: .seealso: KSPCreate(), KSPSolve(), KSPDestroy()
239: @*/
240: PetscErrorCode  KSPSetUp(KSP ksp)
241: {
243:   Mat            A,B;
244:   Mat            mat,pmat;
245:   MatNullSpace   nullsp;
246: 

250:   /* reset the convergence flag from the previous solves */
251:   ksp->reason = KSP_CONVERGED_ITERATING;

253:   if (!((PetscObject)ksp)->type_name) {
254:     KSPSetType(ksp,KSPGMRES);
255:   }
256:   KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);

258:   if (ksp->dmActive && !ksp->setupstage) {
259:     /* first time in so build matrix and vector data structures using DM */
260:     if (!ksp->vec_rhs) {DMCreateGlobalVector(ksp->dm,&ksp->vec_rhs);}
261:     if (!ksp->vec_sol) {DMCreateGlobalVector(ksp->dm,&ksp->vec_sol);}
262:     DMCreateMatrix(ksp->dm,&A);
263:     KSPSetOperators(ksp,A,A);
264:     PetscObjectDereference((PetscObject)A);
265:   }

267:   if (ksp->dmActive) {
268:     DMKSP kdm;
269:     DMGetDMKSP(ksp->dm,&kdm);

271:     if (kdm->ops->computeinitialguess && ksp->setupstage != KSP_SETUP_NEWRHS) {
272:       /* only computes initial guess the first time through */
273:       (*kdm->ops->computeinitialguess)(ksp,ksp->vec_sol,kdm->initialguessctx);
274:       KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
275:     }
276:     if (kdm->ops->computerhs) {
277:       (*kdm->ops->computerhs)(ksp,ksp->vec_rhs,kdm->rhsctx);
278:     }

280:     if (ksp->setupstage != KSP_SETUP_NEWRHS) {
281:       if (kdm->ops->computeoperators) {
282:         KSPGetOperators(ksp,&A,&B);
283:         (*kdm->ops->computeoperators)(ksp,A,B,kdm->operatorsctx);
284:         KSPSetOperators(ksp,A,B);
285:       } else SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONGSTATE,"You called KSPSetDM() but did not use DMKSPSetComputeOperators() or KSPSetDMActive(dm,PETSC_FALSE);");
286:     }
287:   }

289:   if (ksp->setupstage == KSP_SETUP_NEWRHS) return(0);
290:   PetscLogEventBegin(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);

292:   switch (ksp->setupstage) {
293:   case KSP_SETUP_NEW:
294:     (*ksp->ops->setup)(ksp);
295:     break;
296:   case KSP_SETUP_NEWMATRIX: {   /* This should be replaced with a more general mechanism */
297:   } break;
298:   default: break;
299:   }

301:   PCGetOperators(ksp->pc,&mat,&pmat);
302:   /* scale the matrix if requested */
303:   if (ksp->dscale) {
304:     PetscScalar *xx;
305:     PetscInt    i,n;
306:     PetscBool   zeroflag = PETSC_FALSE;
307:     if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
308:     if (!ksp->diagonal) { /* allocate vector to hold diagonal */
309:       MatCreateVecs(pmat,&ksp->diagonal,0);
310:     }
311:     MatGetDiagonal(pmat,ksp->diagonal);
312:     VecGetLocalSize(ksp->diagonal,&n);
313:     VecGetArray(ksp->diagonal,&xx);
314:     for (i=0; i<n; i++) {
315:       if (xx[i] != 0.0) xx[i] = 1.0/PetscSqrtReal(PetscAbsScalar(xx[i]));
316:       else {
317:         xx[i]    = 1.0;
318:         zeroflag = PETSC_TRUE;
319:       }
320:     }
321:     VecRestoreArray(ksp->diagonal,&xx);
322:     if (zeroflag) {
323:       PetscInfo(ksp,"Zero detected in diagonal of matrix, using 1 at those locations\n");
324:     }
325:     MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
326:     if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
327:     ksp->dscalefix2 = PETSC_FALSE;
328:   }
329:   PetscLogEventEnd(KSP_SetUp,ksp,ksp->vec_rhs,ksp->vec_sol,0);
330:   if (!ksp->pc) {KSPGetPC(ksp,&ksp->pc);}
331:   PCSetErrorIfFailure(ksp->pc,ksp->errorifnotconverged);
332:   PCSetUp(ksp->pc);
333:   MatGetNullSpace(mat,&nullsp);
334:   if (nullsp) {
335:     PetscBool test = PETSC_FALSE;
336:     PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_test_null_space",&test,NULL);
337:     if (test) {
338:       MatNullSpaceTest(nullsp,mat,NULL);
339:     }
340:   }
341:   ksp->setupstage = KSP_SETUP_NEWRHS;
342:   return(0);
343: }

347: /*@
348:    KSPReasonView - Displays the reason a KSP solve converged or diverged to a viewer

350:    Collective on KSP

352:    Parameter:
353: +  ksp - iterative context obtained from KSPCreate()
354: -  viewer - the viewer to display the reason


357:    Options Database Keys:
358: .  -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations

360:    Level: beginner

362: .keywords: KSP, solve, linear system

364: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
365:           KSPSolveTranspose(), KSPGetIterationNumber()
366: @*/
367: PetscErrorCode  KSPReasonView(KSP ksp,PetscViewer viewer)
368: {
370:   PetscBool      isAscii;

373:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isAscii);
374:   if (isAscii) {
375:     PetscViewerASCIIAddTab(viewer,((PetscObject)ksp)->tablevel);
376:     if (ksp->reason > 0) {
377:       if (((PetscObject) ksp)->prefix) {
378:         PetscViewerASCIIPrintf(viewer,"Linear %s solve converged due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);
379:       } else {
380:         PetscViewerASCIIPrintf(viewer,"Linear solve converged due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
381:       }
382:     } else {
383:       if (((PetscObject) ksp)->prefix) {
384:         PetscViewerASCIIPrintf(viewer,"Linear %s solve did not converge due to %s iterations %D\n",((PetscObject) ksp)->prefix,KSPConvergedReasons[ksp->reason],ksp->its);
385:       } else {
386:         PetscViewerASCIIPrintf(viewer,"Linear solve did not converge due to %s iterations %D\n",KSPConvergedReasons[ksp->reason],ksp->its);
387:       }
388:     }
389:     PetscViewerASCIISubtractTab(viewer,((PetscObject)ksp)->tablevel);
390:   }
391:   return(0);
392: }

394: #if defined(PETSC_HAVE_THREADSAFETY)
395: #define KSPReasonViewFromOptions KSPReasonViewFromOptionsUnsafe
396: #endif

400: /*@C
401:   KSPReasonViewFromOptions - Processes command line options to determine if/how a KSPReason is to be viewed. 

403:   Collective on KSP

405:   Input Parameters:
406: . ksp   - the KSP object

408:   Level: intermediate

410: @*/
411: PetscErrorCode KSPReasonViewFromOptions(KSP ksp)
412: {
413:   PetscErrorCode    ierr;
414:   PetscViewer       viewer;
415:   PetscBool         flg;
416:   static PetscBool  incall = PETSC_FALSE;
417:   PetscViewerFormat format;

420:   if (incall) return(0);
421:   incall = PETSC_TRUE;
422:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)ksp),((PetscObject)ksp)->prefix,"-ksp_converged_reason",&viewer,&format,&flg);
423:   if (flg) {
424:     PetscViewerPushFormat(viewer,format);
425:     KSPReasonView(ksp,viewer);
426:     PetscViewerPopFormat(viewer);
427:     PetscViewerDestroy(&viewer);
428:   }
429:   incall = PETSC_FALSE;
430:   return(0);
431: }

433: #if defined(PETSC_HAVE_THREADSAFETY)
434: #undef KSPReasonViewFromOptions
435: PetscErrorCode KSPReasonViewFromOptions(KSP ksp)
436: {
438: #pragma omp critical
439:   KSPReasonViewFromOptionsUnsafe(ksp);
440:   return ierr;
441: }
442: #endif

444: #include <petscdraw.h>
447: /*@
448:    KSPSolve - Solves linear system.

450:    Collective on KSP

452:    Parameter:
453: +  ksp - iterative context obtained from KSPCreate()
454: .  b - the right hand side vector
455: -  x - the solution  (this may be the same vector as b, then b will be overwritten with answer)

457:    Options Database Keys:
458: +  -ksp_compute_eigenvalues - compute preconditioned operators eigenvalues
459: .  -ksp_plot_eigenvalues - plot the computed eigenvalues in an X-window
460: .  -ksp_compute_eigenvalues_explicitly - compute the eigenvalues by forming the dense operator and useing LAPACK
461: .  -ksp_plot_eigenvalues_explicitly - plot the explicitly computing eigenvalues
462: .  -ksp_view_mat binary - save matrix to the default binary viewer
463: .  -ksp_view_pmat binary - save matrix used to build preconditioner to the default binary viewer
464: .  -ksp_view_rhs binary - save right hand side vector to the default binary viewer
465: .  -ksp_view_solution binary - save computed solution vector to the default binary viewer
466:            (can be read later with src/ksp/examples/tutorials/ex10.c for testing solvers)
467: .  -ksp_view_mat_explicit - for matrix-free operators, computes the matrix entries and views them
468: .  -ksp_view_preconditioned_operator_explicit - computes the product of the preconditioner and matrix as an explicit matrix and views it
469: .  -ksp_converged_reason - print reason for converged or diverged, also prints number of iterations
470: .  -ksp_final_residual - print 2-norm of true linear system residual at the end of the solution process
471: -  -ksp_view - print the ksp data structure at the end of the system solution

473:    Notes:

475:    If one uses KSPSetDM() then x or b need not be passed. Use KSPGetSolution() to access the solution in this case.

477:    The operator is specified with KSPSetOperators().

479:    Call KSPGetConvergedReason() to determine if the solver converged or failed and
480:    why. The number of iterations can be obtained from KSPGetIterationNumber().

482:    If using a direct method (e.g., via the KSP solver
483:    KSPPREONLY and a preconditioner such as PCLU/PCILU),
484:    then its=1.  See KSPSetTolerances() and KSPConvergedDefault()
485:    for more details.

487:    Understanding Convergence:
488:    The routines KSPMonitorSet(), KSPComputeEigenvalues(), and
489:    KSPComputeEigenvaluesExplicitly() provide information on additional
490:    options to monitor convergence and print eigenvalue information.

492:    Level: beginner

494: .keywords: KSP, solve, linear system

496: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
497:           KSPSolveTranspose(), KSPGetIterationNumber()
498: @*/
499: PetscErrorCode  KSPSolve(KSP ksp,Vec b,Vec x)
500: {
501:   PetscErrorCode    ierr;
502:   PetscMPIInt       rank;
503:   PetscBool         flag1,flag2,flag3,flg = PETSC_FALSE,inXisinB=PETSC_FALSE,guess_zero;
504:   Mat               mat,pmat;
505:   MPI_Comm          comm;
506:   PetscInt          pcreason;
507:   MatNullSpace      nullsp;
508:   Vec               btmp,vec_rhs=0;

514:   comm = PetscObjectComm((PetscObject)ksp);
515:   if (x && x == b) {
516:     if (!ksp->guess_zero) SETERRQ(comm,PETSC_ERR_ARG_INCOMP,"Cannot use x == b with nonzero initial guess");
517:     VecDuplicate(b,&x);
518:     inXisinB = PETSC_TRUE;
519:   }
520:   if (b) {
521:     PetscObjectReference((PetscObject)b);
522:     VecDestroy(&ksp->vec_rhs);
523:     ksp->vec_rhs = b;
524:   }
525:   if (x) {
526:     PetscObjectReference((PetscObject)x);
527:     VecDestroy(&ksp->vec_sol);
528:     ksp->vec_sol = x;
529:   }
530:   KSPViewFromOptions(ksp,NULL,"-ksp_view_pre");

532:   if (ksp->presolve) {
533:     (*ksp->presolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->prectx);
534:   }
535:   PetscLogEventBegin(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);

537:   /* reset the residual history list if requested */
538:   if (ksp->res_hist_reset) ksp->res_hist_len = 0;
539:   ksp->transpose_solve = PETSC_FALSE;

541:   if (ksp->guess) {
542:     KSPFischerGuessFormGuess(ksp->guess,ksp->vec_rhs,ksp->vec_sol);
543:     ksp->guess_zero = PETSC_FALSE;
544:   }
545:   /* KSPSetUp() scales the matrix if needed */
546:   KSPSetUp(ksp);
547:   PCGetSetUpFailedReason(ksp->pc,&pcreason);
548:   if (pcreason) {
549:     ksp->reason = KSP_DIVERGED_PCSETUP_FAILED;
550:     goto skipsolve;
551:   }
552:   KSPSetUpOnBlocks(ksp);
553:   VecLocked(ksp->vec_sol,3);

555:   PCGetOperators(ksp->pc,&mat,&pmat);
556:   /* diagonal scale RHS if called for */
557:   if (ksp->dscale) {
558:     VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
559:     /* second time in, but matrix was scaled back to original */
560:     if (ksp->dscalefix && ksp->dscalefix2) {
561:       Mat mat,pmat;

563:       PCGetOperators(ksp->pc,&mat,&pmat);
564:       MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
565:       if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
566:     }

568:     /*  scale initial guess */
569:     if (!ksp->guess_zero) {
570:       if (!ksp->truediagonal) {
571:         VecDuplicate(ksp->diagonal,&ksp->truediagonal);
572:         VecCopy(ksp->diagonal,ksp->truediagonal);
573:         VecReciprocal(ksp->truediagonal);
574:       }
575:       VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->truediagonal);
576:     }
577:   }
578:   PCPreSolve(ksp->pc,ksp);

580:   if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
581:   if (ksp->guess_knoll) {
582:     PCApply(ksp->pc,ksp->vec_rhs,ksp->vec_sol);
583:     KSP_RemoveNullSpace(ksp,ksp->vec_sol);
584:     ksp->guess_zero = PETSC_FALSE;
585:   }

587:   /* can we mark the initial guess as zero for this solve? */
588:   guess_zero = ksp->guess_zero;
589:   if (!ksp->guess_zero) {
590:     PetscReal norm;

592:     VecNormAvailable(ksp->vec_sol,NORM_2,&flg,&norm);
593:     if (flg && !norm) ksp->guess_zero = PETSC_TRUE;
594:   }
595:   MatGetTransposeNullSpace(pmat,&nullsp);
596:   if (nullsp) {
597:     VecDuplicate(ksp->vec_rhs,&btmp);
598:     VecCopy(ksp->vec_rhs,btmp);
599:     MatNullSpaceRemove(nullsp,btmp);
600:     vec_rhs      = ksp->vec_rhs;
601:     ksp->vec_rhs = btmp;
602:   }
603:   VecLockPush(ksp->vec_rhs);
604:   (*ksp->ops->solve)(ksp);
605:   VecLockPop(ksp->vec_rhs);
606:   if (nullsp) {
607:     ksp->vec_rhs = vec_rhs;
608:     VecDestroy(&btmp);
609:   }

611:   ksp->guess_zero = guess_zero;


614:   if (!ksp->reason) SETERRQ(comm,PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
615:   ksp->totalits += ksp->its;

617:   skipsolve:
618:   KSPReasonViewFromOptions(ksp);
619:   PCPostSolve(ksp->pc,ksp);

621:   /* diagonal scale solution if called for */
622:   if (ksp->dscale) {
623:     VecPointwiseMult(ksp->vec_sol,ksp->vec_sol,ksp->diagonal);
624:     /* unscale right hand side and matrix */
625:     if (ksp->dscalefix) {
626:       Mat mat,pmat;

628:       VecReciprocal(ksp->diagonal);
629:       VecPointwiseMult(ksp->vec_rhs,ksp->vec_rhs,ksp->diagonal);
630:       PCGetOperators(ksp->pc,&mat,&pmat);
631:       MatDiagonalScale(pmat,ksp->diagonal,ksp->diagonal);
632:       if (mat != pmat) {MatDiagonalScale(mat,ksp->diagonal,ksp->diagonal);}
633:       VecReciprocal(ksp->diagonal);
634:       ksp->dscalefix2 = PETSC_TRUE;
635:     }
636:   }
637:   PetscLogEventEnd(KSP_Solve,ksp,ksp->vec_rhs,ksp->vec_sol,0);
638:   if (ksp->postsolve) {
639:     (*ksp->postsolve)(ksp,ksp->vec_rhs,ksp->vec_sol,ksp->postctx);
640:   }

642:   if (ksp->guess) {
643:     KSPFischerGuessUpdate(ksp->guess,ksp->vec_sol);
644:   }

646:   MPI_Comm_rank(comm,&rank);

648:   MatViewFromOptions(mat,(PetscObject)ksp,"-ksp_view_mat");
649:   MatViewFromOptions(pmat,(PetscObject)ksp,"-ksp_view_pmat");
650:   VecViewFromOptions(ksp->vec_rhs,(PetscObject)ksp,"-ksp_view_rhs");

652:   flag1 = PETSC_FALSE;
653:   flag2 = PETSC_FALSE;
654:   flag3 = PETSC_FALSE;
655:   PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues",&flag1,NULL);
656:   PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues",&flag2,NULL);
657:   PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigencontours",&flag3,NULL);
658:   if (flag1 || flag2 || flag3) {
659:     PetscInt  nits,n,i,neig;
660:     PetscReal *r,*c;

662:     KSPGetIterationNumber(ksp,&nits);
663:     n    = nits+2;

665:     if (!nits) {
666:       PetscPrintf(comm,"Zero iterations in solver, cannot approximate any eigenvalues\n");
667:     } else {
668:       PetscMalloc2(n,&r,n,&c);
669:       KSPComputeEigenvalues(ksp,n,r,c,&neig);
670:       if (flag1) {
671:         PetscPrintf(comm,"Iteratively computed eigenvalues\n");
672:         for (i=0; i<neig; i++) {
673:           if (c[i] >= 0.0) {
674:             PetscPrintf(comm,"%g + %gi\n",(double)r[i],(double)c[i]);
675:           } else {
676:             PetscPrintf(comm,"%g - %gi\n",(double)r[i],-(double)c[i]);
677:           }
678:         }
679:       }
680:       if (flag2 && !rank) {
681:         PetscDraw   draw;
682:         PetscDrawSP drawsp;

684:         if (!ksp->eigviewer) {
685:           PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Iteratively Computed Eigenvalues",PETSC_DECIDE,PETSC_DECIDE,400,400,&ksp->eigviewer);
686:         }
687:         PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);
688:         PetscDrawSPCreate(draw,1,&drawsp);
689:         PetscDrawSPReset(drawsp);
690:         for (i=0; i<neig; i++) {
691:           PetscDrawSPAddPoint(drawsp,r+i,c+i);
692:         }
693:         PetscDrawSPDraw(drawsp,PETSC_TRUE);
694:         PetscDrawSPDestroy(&drawsp);
695:       }
696:       if (flag3 && !rank) {
697:         KSPPlotEigenContours_Private(ksp,neig,r,c);
698:       }
699:       PetscFree2(r,c);
700:     }
701:   }

703:   flag1 = PETSC_FALSE;
704:   PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_singularvalues",&flag1,NULL);
705:   if (flag1) {
706:     PetscInt nits;

708:     KSPGetIterationNumber(ksp,&nits);
709:     if (!nits) {
710:       PetscPrintf(comm,"Zero iterations in solver, cannot approximate any singular values\n");
711:     } else {
712:       PetscReal emax,emin;

714:       KSPComputeExtremeSingularValues(ksp,&emax,&emin);
715:       PetscPrintf(comm,"Iteratively computed extreme singular values: max %g min %g max/min %g\n",(double)emax,(double)emin,(double)(emax/emin));
716:     }
717:   }


720:   flag1 = PETSC_FALSE;
721:   flag2 = PETSC_FALSE;
722:   PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_compute_eigenvalues_explicitly",&flag1,NULL);
723:   PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_plot_eigenvalues_explicitly",&flag2,NULL);
724:   if (flag1 || flag2) {
725:     PetscInt  n,i;
726:     PetscReal *r,*c;
727:     VecGetSize(ksp->vec_sol,&n);
728:     PetscMalloc2(n,&r,n,&c);
729:     KSPComputeEigenvaluesExplicitly(ksp,n,r,c);
730:     if (flag1) {
731:       PetscPrintf(comm,"Explicitly computed eigenvalues\n");
732:       for (i=0; i<n; i++) {
733:         if (c[i] >= 0.0) {
734:           PetscPrintf(comm,"%g + %gi\n",(double)r[i],(double)c[i]);
735:         } else {
736:           PetscPrintf(comm,"%g - %gi\n",(double)r[i],-(double)c[i]);
737:         }
738:       }
739:     }
740:     if (flag2 && !rank) {
741:       PetscDraw   draw;
742:       PetscDrawSP drawsp;

744:       if (!ksp->eigviewer) {
745:         PetscViewerDrawOpen(PETSC_COMM_SELF,0,"Explicitly Computed Eigenvalues",0,320,400,400,&ksp->eigviewer);
746:       }
747:       PetscViewerDrawGetDraw(ksp->eigviewer,0,&draw);
748:       PetscDrawSPCreate(draw,1,&drawsp);
749:       PetscDrawSPReset(drawsp);
750:       for (i=0; i<n; i++) {
751:         PetscDrawSPAddPoint(drawsp,r+i,c+i);
752:       }
753:       PetscDrawSPDraw(drawsp,PETSC_TRUE);
754:       PetscDrawSPDestroy(&drawsp);
755:     }
756:     PetscFree2(r,c);
757:   }

759:   PetscOptionsHasName(((PetscObject)ksp)->prefix,"-ksp_view_mat_explicit",&flag2);
760:   if (flag2) {
761:     Mat A,B;
762:     PCGetOperators(ksp->pc,&A,NULL);
763:     MatComputeExplicitOperator(A,&B);
764:     MatViewFromOptions(B,(PetscObject)ksp,"-ksp_view_mat_explicit");
765:     MatDestroy(&B);
766:   }
767:   PetscOptionsHasName(((PetscObject)ksp)->prefix,"-ksp_view_preconditioned_operator_explicit",&flag2);
768:   if (flag2) {
769:     Mat B;
770:     KSPComputeExplicitOperator(ksp,&B);
771:     MatViewFromOptions(B,(PetscObject)ksp,"-ksp_view_preconditioned_operator_explicit");
772:     MatDestroy(&B);
773:   }
774:   KSPViewFromOptions(ksp,NULL,"-ksp_view");

776:   flg  = PETSC_FALSE;
777:   PetscOptionsGetBool(((PetscObject)ksp)->prefix,"-ksp_final_residual",&flg,NULL);
778:   if (flg) {
779:     Mat       A;
780:     Vec       t;
781:     PetscReal norm;
782:     if (ksp->dscale && !ksp->dscalefix) SETERRQ(comm,PETSC_ERR_ARG_WRONGSTATE,"Cannot compute final scale with -ksp_diagonal_scale except also with -ksp_diagonal_scale_fix");
783:     PCGetOperators(ksp->pc,&A,NULL);
784:     VecDuplicate(ksp->vec_rhs,&t);
785:     KSP_MatMult(ksp,A,ksp->vec_sol,t);
786:     VecAYPX(t, -1.0, ksp->vec_rhs);
787:     VecNorm(t,NORM_2,&norm);
788:     VecDestroy(&t);
789:     PetscPrintf(comm,"KSP final norm of residual %g\n",(double)norm);
790:   }
791:   VecViewFromOptions(ksp->vec_sol,(PetscObject)ksp,"-ksp_view_solution");

793:   if (inXisinB) {
794:     VecCopy(x,b);
795:     VecDestroy(&x);
796:   }
797:   PetscObjectSAWsBlock((PetscObject)ksp);
798:   if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(comm,PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
799:   return(0);
800: }

804: /*@
805:    KSPSolveTranspose - Solves the transpose of a linear system.

807:    Collective on KSP

809:    Input Parameter:
810: +  ksp - iterative context obtained from KSPCreate()
811: .  b - right hand side vector
812: -  x - solution vector

814:    Notes: For complex numbers this solve the non-Hermitian transpose system.

816:    Developer Notes: We need to implement a KSPSolveHermitianTranspose()

818:    Level: developer

820: .keywords: KSP, solve, linear system

822: .seealso: KSPCreate(), KSPSetUp(), KSPDestroy(), KSPSetTolerances(), KSPConvergedDefault(),
823:           KSPSolve()
824: @*/

826: PetscErrorCode  KSPSolveTranspose(KSP ksp,Vec b,Vec x)
827: {
829:   PetscBool      inXisinB=PETSC_FALSE;

835:   if (x == b) {
836:     VecDuplicate(b,&x);
837:     inXisinB = PETSC_TRUE;
838:   }
839:   PetscObjectReference((PetscObject)b);
840:   PetscObjectReference((PetscObject)x);
841:   VecDestroy(&ksp->vec_rhs);
842:   VecDestroy(&ksp->vec_sol);

844:   ksp->vec_rhs         = b;
845:   ksp->vec_sol         = x;
846:   ksp->transpose_solve = PETSC_TRUE;

848:   KSPSetUp(ksp);
849:   if (ksp->guess_zero) { VecSet(ksp->vec_sol,0.0);}
850:   (*ksp->ops->solve)(ksp);
851:   if (!ksp->reason) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
852:   KSPReasonViewFromOptions(ksp);
853:   if (inXisinB) {
854:     VecCopy(x,b);
855:     VecDestroy(&x);
856:   }
857:   if (ksp->errorifnotconverged && ksp->reason < 0) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_NOT_CONVERGED,"KSPSolve has not converged");
858:   return(0);
859: }

863: /*@
864:    KSPReset - Resets a KSP context to the kspsetupcalled = 0 state and removes any allocated Vecs and Mats

866:    Collective on KSP

868:    Input Parameter:
869: .  ksp - iterative context obtained from KSPCreate()

871:    Level: beginner

873: .keywords: KSP, destroy

875: .seealso: KSPCreate(), KSPSetUp(), KSPSolve()
876: @*/
877: PetscErrorCode  KSPReset(KSP ksp)
878: {

883:   if (!ksp) return(0);
884:   if (ksp->ops->reset) {
885:     (*ksp->ops->reset)(ksp);
886:   }
887:   if (ksp->pc) {PCReset(ksp->pc);}
888:   KSPFischerGuessDestroy(&ksp->guess);
889:   VecDestroyVecs(ksp->nwork,&ksp->work);
890:   VecDestroy(&ksp->vec_rhs);
891:   VecDestroy(&ksp->vec_sol);
892:   VecDestroy(&ksp->diagonal);
893:   VecDestroy(&ksp->truediagonal);

895:   ksp->setupstage = KSP_SETUP_NEW;
896:   return(0);
897: }

901: /*@
902:    KSPDestroy - Destroys KSP context.

904:    Collective on KSP

906:    Input Parameter:
907: .  ksp - iterative context obtained from KSPCreate()

909:    Level: beginner

911: .keywords: KSP, destroy

913: .seealso: KSPCreate(), KSPSetUp(), KSPSolve()
914: @*/
915: PetscErrorCode  KSPDestroy(KSP *ksp)
916: {
918:   PC             pc;

921:   if (!*ksp) return(0);
923:   if (--((PetscObject)(*ksp))->refct > 0) {*ksp = 0; return(0);}

925:   PetscObjectSAWsViewOff((PetscObject)*ksp);
926:   /*
927:    Avoid a cascading call to PCReset(ksp->pc) from the following call:
928:    PCReset() shouldn't be called from KSPDestroy() as it is unprotected by pc's
929:    refcount (and may be shared, e.g., by other ksps).
930:    */
931:   pc         = (*ksp)->pc;
932:   (*ksp)->pc = NULL;
933:   KSPReset((*ksp));
934:   (*ksp)->pc = pc;
935:     if ((*ksp)->ops->destroy) {(*(*ksp)->ops->destroy)(*ksp);}

937:   DMDestroy(&(*ksp)->dm);
938:   PCDestroy(&(*ksp)->pc);
939:   PetscFree((*ksp)->res_hist_alloc);
940:   if ((*ksp)->convergeddestroy) {
941:     (*(*ksp)->convergeddestroy)((*ksp)->cnvP);
942:   }
943:   KSPMonitorCancel((*ksp));
944:   PetscViewerDestroy(&(*ksp)->eigviewer);
945:   PetscHeaderDestroy(ksp);
946:   return(0);
947: }

951: /*@
952:     KSPSetPCSide - Sets the preconditioning side.

954:     Logically Collective on KSP

956:     Input Parameter:
957: .   ksp - iterative context obtained from KSPCreate()

959:     Output Parameter:
960: .   side - the preconditioning side, where side is one of
961: .vb
962:       PC_LEFT - left preconditioning (default)
963:       PC_RIGHT - right preconditioning
964:       PC_SYMMETRIC - symmetric preconditioning
965: .ve

967:     Options Database Keys:
968: .   -ksp_pc_side <right,left,symmetric>

970:     Notes:
971:     Left preconditioning is used by default for most Krylov methods except KSPFGMRES which only supports right preconditioning.

973:     For methods changing the side of the preconditioner changes the norm type that is used, see KSPSetNormType().

975:     Symmetric preconditioning is currently available only for the KSPQCG method. Note, however, that
976:     symmetric preconditioning can be emulated by using either right or left
977:     preconditioning and a pre or post processing step.

979:     Setting the PC side often affects the default norm type.  See KSPSetNormType() for details.

981:     Level: intermediate

983: .keywords: KSP, set, right, left, symmetric, side, preconditioner, flag

985: .seealso: KSPGetPCSide(), KSPSetNormType(), KSPGetNormType()
986: @*/
987: PetscErrorCode  KSPSetPCSide(KSP ksp,PCSide side)
988: {
992:   ksp->pc_side = ksp->pc_side_set = side;
993:   return(0);
994: }

998: /*@
999:     KSPGetPCSide - Gets the preconditioning side.

1001:     Not Collective

1003:     Input Parameter:
1004: .   ksp - iterative context obtained from KSPCreate()

1006:     Output Parameter:
1007: .   side - the preconditioning side, where side is one of
1008: .vb
1009:       PC_LEFT - left preconditioning (default)
1010:       PC_RIGHT - right preconditioning
1011:       PC_SYMMETRIC - symmetric preconditioning
1012: .ve

1014:     Level: intermediate

1016: .keywords: KSP, get, right, left, symmetric, side, preconditioner, flag

1018: .seealso: KSPSetPCSide()
1019: @*/
1020: PetscErrorCode  KSPGetPCSide(KSP ksp,PCSide *side)
1021: {

1027:   KSPSetUpNorms_Private(ksp,&ksp->normtype,&ksp->pc_side);
1028:   *side = ksp->pc_side;
1029:   return(0);
1030: }

1034: /*@
1035:    KSPGetTolerances - Gets the relative, absolute, divergence, and maximum
1036:    iteration tolerances used by the default KSP convergence tests.

1038:    Not Collective

1040:    Input Parameter:
1041: .  ksp - the Krylov subspace context

1043:    Output Parameters:
1044: +  rtol - the relative convergence tolerance
1045: .  abstol - the absolute convergence tolerance
1046: .  dtol - the divergence tolerance
1047: -  maxits - maximum number of iterations

1049:    Notes:
1050:    The user can specify NULL for any parameter that is not needed.

1052:    Level: intermediate

1054: .keywords: KSP, get, tolerance, absolute, relative, divergence, convergence,
1055:            maximum, iterations

1057: .seealso: KSPSetTolerances()
1058: @*/
1059: PetscErrorCode  KSPGetTolerances(KSP ksp,PetscReal *rtol,PetscReal *abstol,PetscReal *dtol,PetscInt *maxits)
1060: {
1063:   if (abstol) *abstol = ksp->abstol;
1064:   if (rtol) *rtol = ksp->rtol;
1065:   if (dtol) *dtol = ksp->divtol;
1066:   if (maxits) *maxits = ksp->max_it;
1067:   return(0);
1068: }

1072: /*@
1073:    KSPSetTolerances - Sets the relative, absolute, divergence, and maximum
1074:    iteration tolerances used by the default KSP convergence testers.

1076:    Logically Collective on KSP

1078:    Input Parameters:
1079: +  ksp - the Krylov subspace context
1080: .  rtol - the relative convergence tolerance, relative decrease in the (possibly preconditioned) residual norm
1081: .  abstol - the absolute convergence tolerance   absolute size of the (possibly preconditioned) residual norm
1082: .  dtol - the divergence tolerance,   amount (possibly preconditioned) residual norm can increase before KSPConvergedDefault() concludes that the method is diverging
1083: -  maxits - maximum number of iterations to use

1085:    Options Database Keys:
1086: +  -ksp_atol <abstol> - Sets abstol
1087: .  -ksp_rtol <rtol> - Sets rtol
1088: .  -ksp_divtol <dtol> - Sets dtol
1089: -  -ksp_max_it <maxits> - Sets maxits

1091:    Notes:
1092:    Use PETSC_DEFAULT to retain the default value of any of the tolerances.

1094:    See KSPConvergedDefault() for details how these parameters are used in the default convergence test.  See also KSPSetConvergenceTest()
1095:    for setting user-defined stopping criteria.

1097:    Level: intermediate

1099: .keywords: KSP, set, tolerance, absolute, relative, divergence,
1100:            convergence, maximum, iterations

1102: .seealso: KSPGetTolerances(), KSPConvergedDefault(), KSPSetConvergenceTest()
1103: @*/
1104: PetscErrorCode  KSPSetTolerances(KSP ksp,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt maxits)
1105: {

1113:   if (rtol != PETSC_DEFAULT) {
1114:     if (rtol < 0.0 || 1.0 <= rtol) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Relative tolerance %g must be non-negative and less than 1.0",(double)rtol);
1115:     ksp->rtol = rtol;
1116:   }
1117:   if (abstol != PETSC_DEFAULT) {
1118:     if (abstol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Absolute tolerance %g must be non-negative",(double)abstol);
1119:     ksp->abstol = abstol;
1120:   }
1121:   if (dtol != PETSC_DEFAULT) {
1122:     if (dtol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Divergence tolerance %g must be larger than 1.0",(double)dtol);
1123:     ksp->divtol = dtol;
1124:   }
1125:   if (maxits != PETSC_DEFAULT) {
1126:     if (maxits < 0) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Maximum number of iterations %D must be non-negative",maxits);
1127:     ksp->max_it = maxits;
1128:   }
1129:   return(0);
1130: }

1134: /*@
1135:    KSPSetInitialGuessNonzero - Tells the iterative solver that the
1136:    initial guess is nonzero; otherwise KSP assumes the initial guess
1137:    is to be zero (and thus zeros it out before solving).

1139:    Logically Collective on KSP

1141:    Input Parameters:
1142: +  ksp - iterative context obtained from KSPCreate()
1143: -  flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero

1145:    Options database keys:
1146: .  -ksp_initial_guess_nonzero : use nonzero initial guess; this takes an optional truth value (0/1/no/yes/true/false)

1148:    Level: beginner

1150:    Notes:
1151:     If this is not called the X vector is zeroed in the call to KSPSolve().

1153: .keywords: KSP, set, initial guess, nonzero

1155: .seealso: KSPGetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll()
1156: @*/
1157: PetscErrorCode  KSPSetInitialGuessNonzero(KSP ksp,PetscBool flg)
1158: {
1162:   ksp->guess_zero = (PetscBool) !(int)flg;
1163:   return(0);
1164: }

1168: /*@
1169:    KSPGetInitialGuessNonzero - Determines whether the KSP solver is using
1170:    a zero initial guess.

1172:    Not Collective

1174:    Input Parameter:
1175: .  ksp - iterative context obtained from KSPCreate()

1177:    Output Parameter:
1178: .  flag - PETSC_TRUE if guess is nonzero, else PETSC_FALSE

1180:    Level: intermediate

1182: .keywords: KSP, set, initial guess, nonzero

1184: .seealso: KSPSetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll()
1185: @*/
1186: PetscErrorCode  KSPGetInitialGuessNonzero(KSP ksp,PetscBool  *flag)
1187: {
1191:   if (ksp->guess_zero) *flag = PETSC_FALSE;
1192:   else *flag = PETSC_TRUE;
1193:   return(0);
1194: }

1198: /*@
1199:    KSPSetErrorIfNotConverged - Causes KSPSolve() to generate an error if the solver has not converged.

1201:    Logically Collective on KSP

1203:    Input Parameters:
1204: +  ksp - iterative context obtained from KSPCreate()
1205: -  flg - PETSC_TRUE indicates you want the error generated

1207:    Options database keys:
1208: .  -ksp_error_if_not_converged : this takes an optional truth value (0/1/no/yes/true/false)

1210:    Level: intermediate

1212:    Notes:
1213:     Normally PETSc continues if a linear solver fails to converge, you can call KSPGetConvergedReason() after a KSPSolve()
1214:     to determine if it has converged.

1216: .keywords: KSP, set, initial guess, nonzero

1218: .seealso: KSPGetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll(), KSPGetErrorIfNotConverged()
1219: @*/
1220: PetscErrorCode  KSPSetErrorIfNotConverged(KSP ksp,PetscBool flg)
1221: {
1225:   ksp->errorifnotconverged = flg;
1226:   return(0);
1227: }

1231: /*@
1232:    KSPGetErrorIfNotConverged - Will KSPSolve() generate an error if the solver does not converge?

1234:    Not Collective

1236:    Input Parameter:
1237: .  ksp - iterative context obtained from KSPCreate()

1239:    Output Parameter:
1240: .  flag - PETSC_TRUE if it will generate an error, else PETSC_FALSE

1242:    Level: intermediate

1244: .keywords: KSP, set, initial guess, nonzero

1246: .seealso: KSPSetInitialGuessNonzero(), KSPSetInitialGuessKnoll(), KSPGetInitialGuessKnoll(), KSPSetErrorIfNotConverged()
1247: @*/
1248: PetscErrorCode  KSPGetErrorIfNotConverged(KSP ksp,PetscBool  *flag)
1249: {
1253:   *flag = ksp->errorifnotconverged;
1254:   return(0);
1255: }

1259: /*@
1260:    KSPSetInitialGuessKnoll - Tells the iterative solver to use PCApply(pc,b,..) to compute the initial guess (The Knoll trick)

1262:    Logically Collective on KSP

1264:    Input Parameters:
1265: +  ksp - iterative context obtained from KSPCreate()
1266: -  flg - PETSC_TRUE or PETSC_FALSE

1268:    Level: advanced


1271: .keywords: KSP, set, initial guess, nonzero

1273: .seealso: KSPGetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero()
1274: @*/
1275: PetscErrorCode  KSPSetInitialGuessKnoll(KSP ksp,PetscBool flg)
1276: {
1280:   ksp->guess_knoll = flg;
1281:   return(0);
1282: }

1286: /*@
1287:    KSPGetInitialGuessKnoll - Determines whether the KSP solver is using the Knoll trick (using PCApply(pc,b,...) to compute
1288:      the initial guess

1290:    Not Collective

1292:    Input Parameter:
1293: .  ksp - iterative context obtained from KSPCreate()

1295:    Output Parameter:
1296: .  flag - PETSC_TRUE if using Knoll trick, else PETSC_FALSE

1298:    Level: advanced

1300: .keywords: KSP, set, initial guess, nonzero

1302: .seealso: KSPSetInitialGuessKnoll(), KSPSetInitialGuessNonzero(), KSPGetInitialGuessNonzero()
1303: @*/
1304: PetscErrorCode  KSPGetInitialGuessKnoll(KSP ksp,PetscBool  *flag)
1305: {
1309:   *flag = ksp->guess_knoll;
1310:   return(0);
1311: }

1315: /*@
1316:    KSPGetComputeSingularValues - Gets the flag indicating whether the extreme singular
1317:    values will be calculated via a Lanczos or Arnoldi process as the linear
1318:    system is solved.

1320:    Not Collective

1322:    Input Parameter:
1323: .  ksp - iterative context obtained from KSPCreate()

1325:    Output Parameter:
1326: .  flg - PETSC_TRUE or PETSC_FALSE

1328:    Options Database Key:
1329: .  -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()

1331:    Notes:
1332:    Currently this option is not valid for all iterative methods.

1334:    Many users may just want to use the monitoring routine
1335:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1336:    to print the singular values at each iteration of the linear solve.

1338:    Level: advanced

1340: .keywords: KSP, set, compute, singular values

1342: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue()
1343: @*/
1344: PetscErrorCode  KSPGetComputeSingularValues(KSP ksp,PetscBool  *flg)
1345: {
1349:   *flg = ksp->calc_sings;
1350:   return(0);
1351: }

1355: /*@
1356:    KSPSetComputeSingularValues - Sets a flag so that the extreme singular
1357:    values will be calculated via a Lanczos or Arnoldi process as the linear
1358:    system is solved.

1360:    Logically Collective on KSP

1362:    Input Parameters:
1363: +  ksp - iterative context obtained from KSPCreate()
1364: -  flg - PETSC_TRUE or PETSC_FALSE

1366:    Options Database Key:
1367: .  -ksp_monitor_singular_value - Activates KSPSetComputeSingularValues()

1369:    Notes:
1370:    Currently this option is not valid for all iterative methods.

1372:    Many users may just want to use the monitoring routine
1373:    KSPMonitorSingularValue() (which can be set with option -ksp_monitor_singular_value)
1374:    to print the singular values at each iteration of the linear solve.

1376:    Level: advanced

1378: .keywords: KSP, set, compute, singular values

1380: .seealso: KSPComputeExtremeSingularValues(), KSPMonitorSingularValue()
1381: @*/
1382: PetscErrorCode  KSPSetComputeSingularValues(KSP ksp,PetscBool flg)
1383: {
1387:   ksp->calc_sings = flg;
1388:   return(0);
1389: }

1393: /*@
1394:    KSPGetComputeEigenvalues - Gets the flag indicating that the extreme eigenvalues
1395:    values will be calculated via a Lanczos or Arnoldi process as the linear
1396:    system is solved.

1398:    Not Collective

1400:    Input Parameter:
1401: .  ksp - iterative context obtained from KSPCreate()

1403:    Output Parameter:
1404: .  flg - PETSC_TRUE or PETSC_FALSE

1406:    Notes:
1407:    Currently this option is not valid for all iterative methods.

1409:    Level: advanced

1411: .keywords: KSP, set, compute, eigenvalues

1413: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly()
1414: @*/
1415: PetscErrorCode  KSPGetComputeEigenvalues(KSP ksp,PetscBool  *flg)
1416: {
1420:   *flg = ksp->calc_sings;
1421:   return(0);
1422: }

1426: /*@
1427:    KSPSetComputeEigenvalues - Sets a flag so that the extreme eigenvalues
1428:    values will be calculated via a Lanczos or Arnoldi process as the linear
1429:    system is solved.

1431:    Logically Collective on KSP

1433:    Input Parameters:
1434: +  ksp - iterative context obtained from KSPCreate()
1435: -  flg - PETSC_TRUE or PETSC_FALSE

1437:    Notes:
1438:    Currently this option is not valid for all iterative methods.

1440:    Level: advanced

1442: .keywords: KSP, set, compute, eigenvalues

1444: .seealso: KSPComputeEigenvalues(), KSPComputeEigenvaluesExplicitly()
1445: @*/
1446: PetscErrorCode  KSPSetComputeEigenvalues(KSP ksp,PetscBool flg)
1447: {
1451:   ksp->calc_sings = flg;
1452:   return(0);
1453: }

1457: /*@
1458:    KSPGetRhs - Gets the right-hand-side vector for the linear system to
1459:    be solved.

1461:    Not Collective

1463:    Input Parameter:
1464: .  ksp - iterative context obtained from KSPCreate()

1466:    Output Parameter:
1467: .  r - right-hand-side vector

1469:    Level: developer

1471: .keywords: KSP, get, right-hand-side, rhs

1473: .seealso: KSPGetSolution(), KSPSolve()
1474: @*/
1475: PetscErrorCode  KSPGetRhs(KSP ksp,Vec *r)
1476: {
1480:   *r = ksp->vec_rhs;
1481:   return(0);
1482: }

1486: /*@
1487:    KSPGetSolution - Gets the location of the solution for the
1488:    linear system to be solved.  Note that this may not be where the solution
1489:    is stored during the iterative process; see KSPBuildSolution().

1491:    Not Collective

1493:    Input Parameters:
1494: .  ksp - iterative context obtained from KSPCreate()

1496:    Output Parameters:
1497: .  v - solution vector

1499:    Level: developer

1501: .keywords: KSP, get, solution

1503: .seealso: KSPGetRhs(),  KSPBuildSolution(), KSPSolve()
1504: @*/
1505: PetscErrorCode  KSPGetSolution(KSP ksp,Vec *v)
1506: {
1510:   *v = ksp->vec_sol;
1511:   return(0);
1512: }

1516: /*@
1517:    KSPSetPC - Sets the preconditioner to be used to calculate the
1518:    application of the preconditioner on a vector.

1520:    Collective on KSP

1522:    Input Parameters:
1523: +  ksp - iterative context obtained from KSPCreate()
1524: -  pc   - the preconditioner object

1526:    Notes:
1527:    Use KSPGetPC() to retrieve the preconditioner context (for example,
1528:    to free it at the end of the computations).

1530:    Level: developer

1532: .keywords: KSP, set, precondition, Binv

1534: .seealso: KSPGetPC()
1535: @*/
1536: PetscErrorCode  KSPSetPC(KSP ksp,PC pc)
1537: {

1544:   PetscObjectReference((PetscObject)pc);
1545:   PCDestroy(&ksp->pc);
1546:   ksp->pc = pc;
1547:   PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);
1548:   return(0);
1549: }

1553: /*@
1554:    KSPGetPC - Returns a pointer to the preconditioner context
1555:    set with KSPSetPC().

1557:    Not Collective

1559:    Input Parameters:
1560: .  ksp - iterative context obtained from KSPCreate()

1562:    Output Parameter:
1563: .  pc - preconditioner context

1565:    Level: developer

1567: .keywords: KSP, get, preconditioner, Binv

1569: .seealso: KSPSetPC()
1570: @*/
1571: PetscErrorCode  KSPGetPC(KSP ksp,PC *pc)
1572: {

1578:   if (!ksp->pc) {
1579:     PCCreate(PetscObjectComm((PetscObject)ksp),&ksp->pc);
1580:     PetscObjectIncrementTabLevel((PetscObject)ksp->pc,(PetscObject)ksp,0);
1581:     PetscLogObjectParent((PetscObject)ksp,(PetscObject)ksp->pc);
1582:   }
1583:   *pc = ksp->pc;
1584:   return(0);
1585: }

1589: /*@
1590:    KSPMonitor - runs the user provided monitor routines, if they exist

1592:    Collective on KSP

1594:    Input Parameters:
1595: +  ksp - iterative context obtained from KSPCreate()
1596: .  it - iteration number
1597: -  rnorm - relative norm of the residual

1599:    Notes:
1600:    This routine is called by the KSP implementations.
1601:    It does not typically need to be called by the user.

1603:    Level: developer

1605: .seealso: KSPMonitorSet()
1606: @*/
1607: PetscErrorCode KSPMonitor(KSP ksp,PetscInt it,PetscReal rnorm)
1608: {
1609:   PetscInt       i, n = ksp->numbermonitors;

1613:   for (i=0; i<n; i++) {
1614:     (*ksp->monitor[i])(ksp,it,rnorm,ksp->monitorcontext[i]);
1615:   }
1616:   return(0);
1617: }

1621: /*@C
1622:    KSPMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
1623:    the residual/error etc.

1625:    Logically Collective on KSP

1627:    Input Parameters:
1628: +  ksp - iterative context obtained from KSPCreate()
1629: .  monitor - pointer to function (if this is NULL, it turns off monitoring
1630: .  mctx    - [optional] context for private data for the
1631:              monitor routine (use NULL if no context is desired)
1632: -  monitordestroy - [optional] routine that frees monitor context
1633:           (may be NULL)

1635:    Calling Sequence of monitor:
1636: $     monitor (KSP ksp, int it, PetscReal rnorm, void *mctx)

1638: +  ksp - iterative context obtained from KSPCreate()
1639: .  it - iteration number
1640: .  rnorm - (estimated) 2-norm of (preconditioned) residual
1641: -  mctx  - optional monitoring context, as set by KSPMonitorSet()

1643:    Options Database Keys:
1644: +    -ksp_monitor        - sets KSPMonitorDefault()
1645: .    -ksp_monitor_true_residual    - sets KSPMonitorTrueResidualNorm()
1646: .    -ksp_monitor_max    - sets KSPMonitorTrueResidualMaxNorm()
1647: .    -ksp_monitor_lg_residualnorm    - sets line graph monitor,
1648:                            uses KSPMonitorLGResidualNormCreate()
1649: .    -ksp_monitor_lg_true_residualnorm   - sets line graph monitor,
1650:                            uses KSPMonitorLGResidualNormCreate()
1651: .    -ksp_monitor_singular_value    - sets KSPMonitorSingularValue()
1652: -    -ksp_monitor_cancel - cancels all monitors that have
1653:                           been hardwired into a code by
1654:                           calls to KSPMonitorSet(), but
1655:                           does not cancel those set via
1656:                           the options database.

1658:    Notes:
1659:    The default is to do nothing.  To print the residual, or preconditioned
1660:    residual if KSPSetNormType(ksp,KSP_NORM_PRECONDITIONED) was called, use
1661:    KSPMonitorDefault() as the monitoring routine, with a null monitoring
1662:    context.

1664:    Several different monitoring routines may be set by calling
1665:    KSPMonitorSet() multiple times; all will be called in the
1666:    order in which they were set.

1668:    Fortran notes: Only a single monitor function can be set for each KSP object

1670:    Level: beginner

1672: .keywords: KSP, set, monitor

1674: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorCancel()
1675: @*/
1676: PetscErrorCode  KSPMonitorSet(KSP ksp,PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,void*),void *mctx,PetscErrorCode (*monitordestroy)(void**))
1677: {
1678:   PetscInt       i;

1683:   if (ksp->numbermonitors >= MAXKSPMONITORS) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_OUTOFRANGE,"Too many KSP monitors set");
1684:   for (i=0; i<ksp->numbermonitors;i++) {
1685:     if (monitor == ksp->monitor[i] && monitordestroy == ksp->monitordestroy[i] && mctx == ksp->monitorcontext[i]) {
1686:       if (monitordestroy) {
1687:         (*monitordestroy)(&mctx);
1688:       }
1689:       return(0);
1690:     }
1691:   }
1692:   ksp->monitor[ksp->numbermonitors]          = monitor;
1693:   ksp->monitordestroy[ksp->numbermonitors]   = monitordestroy;
1694:   ksp->monitorcontext[ksp->numbermonitors++] = (void*)mctx;
1695:   return(0);
1696: }

1700: /*@
1701:    KSPMonitorCancel - Clears all monitors for a KSP object.

1703:    Logically Collective on KSP

1705:    Input Parameters:
1706: .  ksp - iterative context obtained from KSPCreate()

1708:    Options Database Key:
1709: .  -ksp_monitor_cancel - Cancels all monitors that have
1710:     been hardwired into a code by calls to KSPMonitorSet(),
1711:     but does not cancel those set via the options database.

1713:    Level: intermediate

1715: .keywords: KSP, set, monitor

1717: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate(), KSPMonitorSet()
1718: @*/
1719: PetscErrorCode  KSPMonitorCancel(KSP ksp)
1720: {
1722:   PetscInt       i;

1726:   for (i=0; i<ksp->numbermonitors; i++) {
1727:     if (ksp->monitordestroy[i]) {
1728:       (*ksp->monitordestroy[i])(&ksp->monitorcontext[i]);
1729:     }
1730:   }
1731:   ksp->numbermonitors = 0;
1732:   return(0);
1733: }

1737: /*@C
1738:    KSPGetMonitorContext - Gets the monitoring context, as set by
1739:    KSPMonitorSet() for the FIRST monitor only.

1741:    Not Collective

1743:    Input Parameter:
1744: .  ksp - iterative context obtained from KSPCreate()

1746:    Output Parameter:
1747: .  ctx - monitoring context

1749:    Level: intermediate

1751: .keywords: KSP, get, monitor, context

1753: .seealso: KSPMonitorDefault(), KSPMonitorLGResidualNormCreate()
1754: @*/
1755: PetscErrorCode  KSPGetMonitorContext(KSP ksp,void **ctx)
1756: {
1759:   *ctx =      (ksp->monitorcontext[0]);
1760:   return(0);
1761: }

1765: /*@
1766:    KSPSetResidualHistory - Sets the array used to hold the residual history.
1767:    If set, this array will contain the residual norms computed at each
1768:    iteration of the solver.

1770:    Not Collective

1772:    Input Parameters:
1773: +  ksp - iterative context obtained from KSPCreate()
1774: .  a   - array to hold history
1775: .  na  - size of a
1776: -  reset - PETSC_TRUE indicates the history counter is reset to zero
1777:            for each new linear solve

1779:    Level: advanced

1781:    Notes: The array is NOT freed by PETSc so the user needs to keep track of
1782:            it and destroy once the KSP object is destroyed.

1784:    If 'a' is NULL then space is allocated for the history. If 'na' PETSC_DECIDE or PETSC_DEFAULT then a
1785:    default array of length 10000 is allocated.

1787: .keywords: KSP, set, residual, history, norm

1789: .seealso: KSPGetResidualHistory()

1791: @*/
1792: PetscErrorCode  KSPSetResidualHistory(KSP ksp,PetscReal a[],PetscInt na,PetscBool reset)
1793: {


1799:   PetscFree(ksp->res_hist_alloc);
1800:   if (na != PETSC_DECIDE && na != PETSC_DEFAULT && a) {
1801:     ksp->res_hist     = a;
1802:     ksp->res_hist_max = na;
1803:   } else {
1804:     if (na != PETSC_DECIDE && na != PETSC_DEFAULT) ksp->res_hist_max = na;
1805:     else                                           ksp->res_hist_max = 10000; /* like default ksp->max_it */
1806:     PetscCalloc1(ksp->res_hist_max,&ksp->res_hist_alloc);

1808:     ksp->res_hist = ksp->res_hist_alloc;
1809:   }
1810:   ksp->res_hist_len   = 0;
1811:   ksp->res_hist_reset = reset;
1812:   return(0);
1813: }

1817: /*@C
1818:    KSPGetResidualHistory - Gets the array used to hold the residual history
1819:    and the number of residuals it contains.

1821:    Not Collective

1823:    Input Parameter:
1824: .  ksp - iterative context obtained from KSPCreate()

1826:    Output Parameters:
1827: +  a   - pointer to array to hold history (or NULL)
1828: -  na  - number of used entries in a (or NULL)

1830:    Level: advanced

1832:    Notes:
1833:      Can only be called after a KSPSetResidualHistory() otherwise a and na are set to zero

1835:      The Fortran version of this routine has a calling sequence
1836: $   call KSPGetResidualHistory(KSP ksp, integer na, integer ierr)
1837:     note that you have passed a Fortran array into KSPSetResidualHistory() and you need
1838:     to access the residual values from this Fortran array you provided. Only the na (number of
1839:     residual norms currently held) is set.

1841: .keywords: KSP, get, residual, history, norm

1843: .seealso: KSPGetResidualHistory()

1845: @*/
1846: PetscErrorCode  KSPGetResidualHistory(KSP ksp,PetscReal *a[],PetscInt *na)
1847: {
1850:   if (a) *a = ksp->res_hist;
1851:   if (na) *na = ksp->res_hist_len;
1852:   return(0);
1853: }

1857: /*@C
1858:    KSPSetConvergenceTest - Sets the function to be used to determine
1859:    convergence.

1861:    Logically Collective on KSP

1863:    Input Parameters:
1864: +  ksp - iterative context obtained from KSPCreate()
1865: .  converge - pointer to int function
1866: .  cctx    - context for private data for the convergence routine (may be null)
1867: -  destroy - a routine for destroying the context (may be null)

1869:    Calling sequence of converge:
1870: $     converge (KSP ksp, int it, PetscReal rnorm, KSPConvergedReason *reason,void *mctx)

1872: +  ksp - iterative context obtained from KSPCreate()
1873: .  it - iteration number
1874: .  rnorm - (estimated) 2-norm of (preconditioned) residual
1875: .  reason - the reason why it has converged or diverged
1876: -  cctx  - optional convergence context, as set by KSPSetConvergenceTest()


1879:    Notes:
1880:    Must be called after the KSP type has been set so put this after
1881:    a call to KSPSetType(), or KSPSetFromOptions().

1883:    The default convergence test, KSPConvergedDefault(), aborts if the
1884:    residual grows to more than 10000 times the initial residual.

1886:    The default is a combination of relative and absolute tolerances.
1887:    The residual value that is tested may be an approximation; routines
1888:    that need exact values should compute them.

1890:    In the default PETSc convergence test, the precise values of reason
1891:    are macros such as KSP_CONVERGED_RTOL, which are defined in petscksp.h.

1893:    Level: advanced

1895: .keywords: KSP, set, convergence, test, context

1897: .seealso: KSPConvergedDefault(), KSPGetConvergenceContext(), KSPSetTolerances()
1898: @*/
1899: PetscErrorCode  KSPSetConvergenceTest(KSP ksp,PetscErrorCode (*converge)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*),void *cctx,PetscErrorCode (*destroy)(void*))
1900: {

1905:   if (ksp->convergeddestroy) {
1906:     (*ksp->convergeddestroy)(ksp->cnvP);
1907:   }
1908:   ksp->converged        = converge;
1909:   ksp->convergeddestroy = destroy;
1910:   ksp->cnvP             = (void*)cctx;
1911:   return(0);
1912: }

1916: /*@C
1917:    KSPGetConvergenceContext - Gets the convergence context set with
1918:    KSPSetConvergenceTest().

1920:    Not Collective

1922:    Input Parameter:
1923: .  ksp - iterative context obtained from KSPCreate()

1925:    Output Parameter:
1926: .  ctx - monitoring context

1928:    Level: advanced

1930: .keywords: KSP, get, convergence, test, context

1932: .seealso: KSPConvergedDefault(), KSPSetConvergenceTest()
1933: @*/
1934: PetscErrorCode  KSPGetConvergenceContext(KSP ksp,void **ctx)
1935: {
1938:   *ctx = ksp->cnvP;
1939:   return(0);
1940: }

1944: /*@C
1945:    KSPBuildSolution - Builds the approximate solution in a vector provided.
1946:    This routine is NOT commonly needed (see KSPSolve()).

1948:    Collective on KSP

1950:    Input Parameter:
1951: .  ctx - iterative context obtained from KSPCreate()

1953:    Output Parameter:
1954:    Provide exactly one of
1955: +  v - location to stash solution.
1956: -  V - the solution is returned in this location. This vector is created
1957:        internally. This vector should NOT be destroyed by the user with
1958:        VecDestroy().

1960:    Notes:
1961:    This routine can be used in one of two ways
1962: .vb
1963:       KSPBuildSolution(ksp,NULL,&V);
1964:    or
1965:       KSPBuildSolution(ksp,v,NULL); or KSPBuildSolution(ksp,v,&v);
1966: .ve
1967:    In the first case an internal vector is allocated to store the solution
1968:    (the user cannot destroy this vector). In the second case the solution
1969:    is generated in the vector that the user provides. Note that for certain
1970:    methods, such as KSPCG, the second case requires a copy of the solution,
1971:    while in the first case the call is essentially free since it simply
1972:    returns the vector where the solution already is stored. For some methods
1973:    like GMRES this is a reasonably expensive operation and should only be
1974:    used in truly needed.

1976:    Level: advanced

1978: .keywords: KSP, build, solution

1980: .seealso: KSPGetSolution(), KSPBuildResidual()
1981: @*/
1982: PetscErrorCode  KSPBuildSolution(KSP ksp,Vec v,Vec *V)
1983: {

1988:   if (!V && !v) SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_ARG_WRONG,"Must provide either v or V");
1989:   if (!V) V = &v;
1990:   (*ksp->ops->buildsolution)(ksp,v,V);
1991:   return(0);
1992: }

1996: /*@C
1997:    KSPBuildResidual - Builds the residual in a vector provided.

1999:    Collective on KSP

2001:    Input Parameter:
2002: .  ksp - iterative context obtained from KSPCreate()

2004:    Output Parameters:
2005: +  v - optional location to stash residual.  If v is not provided,
2006:        then a location is generated.
2007: .  t - work vector.  If not provided then one is generated.
2008: -  V - the residual

2010:    Notes:
2011:    Regardless of whether or not v is provided, the residual is
2012:    returned in V.

2014:    Level: advanced

2016: .keywords: KSP, build, residual

2018: .seealso: KSPBuildSolution()
2019: @*/
2020: PetscErrorCode  KSPBuildResidual(KSP ksp,Vec t,Vec v,Vec *V)
2021: {
2023:   PetscBool      flag = PETSC_FALSE;
2024:   Vec            w    = v,tt = t;

2028:   if (!w) {
2029:     VecDuplicate(ksp->vec_rhs,&w);
2030:     PetscLogObjectParent((PetscObject)ksp,(PetscObject)w);
2031:   }
2032:   if (!tt) {
2033:     VecDuplicate(ksp->vec_sol,&tt); flag = PETSC_TRUE;
2034:     PetscLogObjectParent((PetscObject)ksp,(PetscObject)tt);
2035:   }
2036:   (*ksp->ops->buildresidual)(ksp,tt,w,V);
2037:   if (flag) {VecDestroy(&tt);}
2038:   return(0);
2039: }

2043: /*@
2044:    KSPSetDiagonalScale - Tells KSP to symmetrically diagonally scale the system
2045:      before solving. This actually CHANGES the matrix (and right hand side).

2047:    Logically Collective on KSP

2049:    Input Parameter:
2050: +  ksp - the KSP context
2051: -  scale - PETSC_TRUE or PETSC_FALSE

2053:    Options Database Key:
2054: +   -ksp_diagonal_scale -
2055: -   -ksp_diagonal_scale_fix - scale the matrix back AFTER the solve


2058:     Notes: Scales the matrix by  D^(-1/2)  A  D^(-1/2)  [D^(1/2) x ] = D^(-1/2) b
2059:        where D_{ii} is 1/abs(A_{ii}) unless A_{ii} is zero and then it is 1.

2061:     BE CAREFUL with this routine: it actually scales the matrix and right
2062:     hand side that define the system. After the system is solved the matrix
2063:     and right hand side remain scaled unless you use KSPSetDiagonalScaleFix()

2065:     This should NOT be used within the SNES solves if you are using a line
2066:     search.

2068:     If you use this with the PCType Eisenstat preconditioner than you can
2069:     use the PCEisenstatSetNoDiagonalScaling() option, or -pc_eisenstat_no_diagonal_scaling
2070:     to save some unneeded, redundant flops.

2072:    Level: intermediate

2074: .keywords: KSP, set, options, prefix, database

2076: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScaleFix()
2077: @*/
2078: PetscErrorCode  KSPSetDiagonalScale(KSP ksp,PetscBool scale)
2079: {
2083:   ksp->dscale = scale;
2084:   return(0);
2085: }

2089: /*@
2090:    KSPGetDiagonalScale - Checks if KSP solver scales the matrix and
2091:                           right hand side

2093:    Not Collective

2095:    Input Parameter:
2096: .  ksp - the KSP context

2098:    Output Parameter:
2099: .  scale - PETSC_TRUE or PETSC_FALSE

2101:    Notes:
2102:     BE CAREFUL with this routine: it actually scales the matrix and right
2103:     hand side that define the system. After the system is solved the matrix
2104:     and right hand side remain scaled  unless you use KSPSetDiagonalScaleFix()

2106:    Level: intermediate

2108: .keywords: KSP, set, options, prefix, database

2110: .seealso: KSPSetDiagonalScale(), KSPSetDiagonalScaleFix()
2111: @*/
2112: PetscErrorCode  KSPGetDiagonalScale(KSP ksp,PetscBool  *scale)
2113: {
2117:   *scale = ksp->dscale;
2118:   return(0);
2119: }

2123: /*@
2124:    KSPSetDiagonalScaleFix - Tells KSP to diagonally scale the system
2125:      back after solving.

2127:    Logically Collective on KSP

2129:    Input Parameter:
2130: +  ksp - the KSP context
2131: -  fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2132:          rescale (default)

2134:    Notes:
2135:      Must be called after KSPSetDiagonalScale()

2137:      Using this will slow things down, because it rescales the matrix before and
2138:      after each linear solve. This is intended mainly for testing to allow one
2139:      to easily get back the original system to make sure the solution computed is
2140:      accurate enough.

2142:    Level: intermediate

2144: .keywords: KSP, set, options, prefix, database

2146: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPGetDiagonalScaleFix()
2147: @*/
2148: PetscErrorCode  KSPSetDiagonalScaleFix(KSP ksp,PetscBool fix)
2149: {
2153:   ksp->dscalefix = fix;
2154:   return(0);
2155: }

2159: /*@
2160:    KSPGetDiagonalScaleFix - Determines if KSP diagonally scales the system
2161:      back after solving.

2163:    Not Collective

2165:    Input Parameter:
2166: .  ksp - the KSP context

2168:    Output Parameter:
2169: .  fix - PETSC_TRUE to scale back after the system solve, PETSC_FALSE to not
2170:          rescale (default)

2172:    Notes:
2173:      Must be called after KSPSetDiagonalScale()

2175:      If PETSC_TRUE will slow things down, because it rescales the matrix before and
2176:      after each linear solve. This is intended mainly for testing to allow one
2177:      to easily get back the original system to make sure the solution computed is
2178:      accurate enough.

2180:    Level: intermediate

2182: .keywords: KSP, set, options, prefix, database

2184: .seealso: KSPGetDiagonalScale(), KSPSetDiagonalScale(), KSPSetDiagonalScaleFix()
2185: @*/
2186: PetscErrorCode  KSPGetDiagonalScaleFix(KSP ksp,PetscBool  *fix)
2187: {
2191:   *fix = ksp->dscalefix;
2192:   return(0);
2193: }

2197: /*@C
2198:    KSPSetComputeOperators - set routine to compute the linear operators

2200:    Logically Collective

2202:    Input Arguments:
2203: +  ksp - the KSP context
2204: .  func - function to compute the operators
2205: -  ctx - optional context

2207:    Calling sequence of func:
2208: $  func(KSP ksp,Mat A,Mat B,void *ctx)

2210: +  ksp - the KSP context
2211: .  A - the linear operator
2212: .  B - preconditioning matrix
2213: -  ctx - optional user-provided context

2215:    Notes: The user provided func() will be called automatically at the very next call to KSPSolve(). It will not be called at future KSPSolve() calls
2216:           unless either KSPSetComputeOperators() or KSPSetOperators() is called before that KSPSolve() is called.

2218:           To reuse the same preconditioner for the next KSPSolve() and not compute a new one based on the most recently computed matrix call KSPSetReusePreconditioner()

2220:    Level: beginner

2222: .seealso: KSPSetOperators(), KSPSetComputeRHS(), DMKSPSetComputeOperators(), KSPSetComputeInitialGuess()
2223: @*/
2224: PetscErrorCode KSPSetComputeOperators(KSP ksp,PetscErrorCode (*func)(KSP,Mat,Mat,void*),void *ctx)
2225: {
2227:   DM             dm;

2231:   KSPGetDM(ksp,&dm);
2232:   DMKSPSetComputeOperators(dm,func,ctx);
2233:   if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX;
2234:   return(0);
2235: }

2239: /*@C
2240:    KSPSetComputeRHS - set routine to compute the right hand side of the linear system

2242:    Logically Collective

2244:    Input Arguments:
2245: +  ksp - the KSP context
2246: .  func - function to compute the right hand side
2247: -  ctx - optional context

2249:    Calling sequence of func:
2250: $  func(KSP ksp,Vec b,void *ctx)

2252: +  ksp - the KSP context
2253: .  b - right hand side of linear system
2254: -  ctx - optional user-provided context

2256:    Notes: The routine you provide will be called EACH you call KSPSolve() to prepare the new right hand side for that solve

2258:    Level: beginner

2260: .seealso: KSPSolve(), DMKSPSetComputeRHS(), KSPSetComputeOperators()
2261: @*/
2262: PetscErrorCode KSPSetComputeRHS(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2263: {
2265:   DM             dm;

2269:   KSPGetDM(ksp,&dm);
2270:   DMKSPSetComputeRHS(dm,func,ctx);
2271:   return(0);
2272: }

2276: /*@C
2277:    KSPSetComputeInitialGuess - set routine to compute the initial guess of the linear system

2279:    Logically Collective

2281:    Input Arguments:
2282: +  ksp - the KSP context
2283: .  func - function to compute the initial guess
2284: -  ctx - optional context

2286:    Calling sequence of func:
2287: $  func(KSP ksp,Vec x,void *ctx)

2289: +  ksp - the KSP context
2290: .  x - solution vector
2291: -  ctx - optional user-provided context

2293:    Level: beginner

2295: .seealso: KSPSolve(), KSPSetComputeRHS(), KSPSetComputeOperators(), DMKSPSetComputeInitialGuess()
2296: @*/
2297: PetscErrorCode KSPSetComputeInitialGuess(KSP ksp,PetscErrorCode (*func)(KSP,Vec,void*),void *ctx)
2298: {
2300:   DM             dm;

2304:   KSPGetDM(ksp,&dm);
2305:   DMKSPSetComputeInitialGuess(dm,func,ctx);
2306:   return(0);
2307: }