Actual source code: mfnsolve.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */
 10: /*
 11:    MFN routines related to the solution process
 12: */

 14: #include <slepc/private/mfnimpl.h>

 16: static PetscErrorCode MFNSolve_Private(MFN mfn,Vec b,Vec x)
 17: {

 21:   VecSetErrorIfLocked(x,3);

 23:   /* call setup */
 24:   MFNSetUp(mfn);
 25:   mfn->its = 0;

 27:   MFNViewFromOptions(mfn,NULL,"-mfn_view_pre");

 29:   /* check nonzero right-hand side */
 30:   VecNorm(b,NORM_2,&mfn->bnorm);
 31:   if (!mfn->bnorm) SETERRQ(PetscObjectComm((PetscObject)mfn),PETSC_ERR_ARG_WRONG,"Cannot pass a zero b vector to MFNSolve()");

 33:   /* call solver */
 34:   PetscLogEventBegin(MFN_Solve,mfn,b,x,0);
 35:   if (b!=x) { VecLockReadPush(b); }
 36:   (*mfn->ops->solve)(mfn,b,x);
 37:   if (b!=x) { VecLockReadPop(b); }
 38:   PetscLogEventEnd(MFN_Solve,mfn,b,x,0);

 40:   if (!mfn->reason) SETERRQ(PetscObjectComm((PetscObject)mfn),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");

 42:   if (mfn->errorifnotconverged && mfn->reason < 0) SETERRQ(PetscObjectComm((PetscObject)mfn),PETSC_ERR_NOT_CONVERGED,"MFNSolve has not converged");

 44:   /* various viewers */
 45:   MFNViewFromOptions(mfn,NULL,"-mfn_view");
 46:   MFNConvergedReasonViewFromOptions(mfn);
 47:   MatViewFromOptions(mfn->A,(PetscObject)mfn,"-mfn_view_mat");
 48:   VecViewFromOptions(b,(PetscObject)mfn,"-mfn_view_rhs");
 49:   VecViewFromOptions(x,(PetscObject)mfn,"-mfn_view_solution");
 50:   return(0);
 51: }

 53: /*@
 54:    MFNSolve - Solves the matrix function problem. Given a vector b, the
 55:    vector x = f(A)*b is returned.

 57:    Collective on mfn

 59:    Input Parameters:
 60: +  mfn - matrix function context obtained from MFNCreate()
 61: -  b   - the right hand side vector

 63:    Output Parameter:
 64: .  x   - the solution (this may be the same vector as b, then b will be
 65:          overwritten with the answer)

 67:    Options Database Keys:
 68: +  -mfn_view - print information about the solver used
 69: .  -mfn_view_mat binary - save the matrix to the default binary viewer
 70: .  -mfn_view_rhs binary - save right hand side vector to the default binary viewer
 71: .  -mfn_view_solution binary - save computed solution vector to the default binary viewer
 72: -  -mfn_converged_reason - print reason for convergence, and number of iterations

 74:    Notes:
 75:    The matrix A is specified with MFNSetOperator().
 76:    The function f is specified with MFNSetFN().

 78:    Level: beginner

 80: .seealso: MFNCreate(), MFNSetUp(), MFNDestroy(), MFNSetTolerances(),
 81:           MFNSetOperator(), MFNSetFN()
 82: @*/
 83: PetscErrorCode MFNSolve(MFN mfn,Vec b,Vec x)
 84: {

 93:   mfn->transpose_solve = PETSC_FALSE;
 94:   MFNSolve_Private(mfn,b,x);
 95:   return(0);
 96: }

 98: /*@
 99:    MFNSolveTranspose - Solves the transpose matrix function problem. Given a vector b,
100:    the vector x = f(A^T)*b is returned.

102:    Collective on mfn

104:    Input Parameters:
105: +  mfn - matrix function context obtained from MFNCreate()
106: -  b   - the right hand side vector

108:    Output Parameter:
109: .  x   - the solution (this may be the same vector as b, then b will be
110:          overwritten with the answer)

112:    Note:
113:    See available options at MFNSolve().

115:    Level: beginner

117: .seealso: MFNSolve()
118: @*/
119: PetscErrorCode MFNSolveTranspose(MFN mfn,Vec b,Vec x)
120: {

129:   mfn->transpose_solve = PETSC_TRUE;
130:   if (!mfn->AT) { MatCreateTranspose(mfn->A,&mfn->AT); }
131:   MFNSolve_Private(mfn,b,x);
132:   return(0);
133: }

135: /*@
136:    MFNGetIterationNumber - Gets the current iteration number. If the
137:    call to MFNSolve() is complete, then it returns the number of iterations
138:    carried out by the solution method.

140:    Not Collective

142:    Input Parameter:
143: .  mfn - the matrix function context

145:    Output Parameter:
146: .  its - number of iterations

148:    Note:
149:    During the i-th iteration this call returns i-1. If MFNSolve() is
150:    complete, then parameter "its" contains either the iteration number at
151:    which convergence was successfully reached, or failure was detected.
152:    Call MFNGetConvergedReason() to determine if the solver converged or
153:    failed and why.

155:    Level: intermediate

157: .seealso: MFNGetConvergedReason(), MFNSetTolerances()
158: @*/
159: PetscErrorCode MFNGetIterationNumber(MFN mfn,PetscInt *its)
160: {
164:   *its = mfn->its;
165:   return(0);
166: }

168: /*@
169:    MFNGetConvergedReason - Gets the reason why the MFNSolve() iteration was
170:    stopped.

172:    Not Collective

174:    Input Parameter:
175: .  mfn - the matrix function context

177:    Output Parameter:
178: .  reason - negative value indicates diverged, positive value converged

180:    Notes:

182:    Possible values for reason are
183: +  MFN_CONVERGED_TOL - converged up to tolerance
184: .  MFN_CONVERGED_ITS - solver completed the requested number of steps
185: .  MFN_DIVERGED_ITS - required more than max_it iterations to reach convergence
186: -  MFN_DIVERGED_BREAKDOWN - generic breakdown in method

188:    Can only be called after the call to MFNSolve() is complete.

190:    Basic solvers (e.g. unrestarted Krylov iterations) cannot determine if the
191:    computation is accurate up to the requested tolerance. In that case, the
192:    converged reason is set to MFN_CONVERGED_ITS if the requested number of steps
193:    (for instance, the ncv value in unrestarted Krylov methods) have been
194:    completed successfully.

196:    Level: intermediate

198: .seealso: MFNSetTolerances(), MFNSolve(), MFNConvergedReason, MFNSetErrorIfNotConverged()
199: @*/
200: PetscErrorCode MFNGetConvergedReason(MFN mfn,MFNConvergedReason *reason)
201: {
205:   *reason = mfn->reason;
206:   return(0);
207: }