1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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: {
 18:   PetscFunctionBegin;
 19:   PetscCall(VecSetErrorIfLocked(x,3));

 21:   /* call setup */
 22:   PetscCall(MFNSetUp(mfn));
 23:   mfn->its = 0;

 25:   PetscCall(MFNViewFromOptions(mfn,NULL,"-mfn_view_pre"));

 27:   /* check nonzero right-hand side */
 28:   PetscCall(VecNorm(b,NORM_2,&mfn->bnorm));
 29:   PetscCheck(mfn->bnorm,PetscObjectComm((PetscObject)mfn),PETSC_ERR_ARG_WRONG,"Cannot pass a zero b vector to MFNSolve()");

 31:   /* call solver */
 32:   PetscCall(PetscLogEventBegin(MFN_Solve,mfn,b,x,0));
 33:   if (b!=x) PetscCall(VecLockReadPush(b));
 34:   PetscUseTypeMethod(mfn,solve,b,x);
 35:   if (b!=x) PetscCall(VecLockReadPop(b));
 36:   PetscCall(PetscLogEventEnd(MFN_Solve,mfn,b,x,0));

 38:   PetscCheck(mfn->reason,PetscObjectComm((PetscObject)mfn),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");

 40:   PetscCheck(!mfn->errorifnotconverged || mfn->reason>=0,PetscObjectComm((PetscObject)mfn),PETSC_ERR_NOT_CONVERGED,"MFNSolve has not converged");

 42:   /* various viewers */
 43:   PetscCall(MFNViewFromOptions(mfn,NULL,"-mfn_view"));
 44:   PetscCall(MFNConvergedReasonViewFromOptions(mfn));
 45:   PetscCall(MatViewFromOptions(mfn->A,(PetscObject)mfn,"-mfn_view_mat"));
 46:   PetscCall(VecViewFromOptions(b,(PetscObject)mfn,"-mfn_view_rhs"));
 47:   PetscCall(VecViewFromOptions(x,(PetscObject)mfn,"-mfn_view_solution"));
 48:   PetscFunctionReturn(PETSC_SUCCESS);
 49: }

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

 55:    Collective

 57:    Input Parameters:
 58: +  mfn - the matrix function solver context
 59: -  b   - the right hand side vector

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

 65:    Options Database Keys:
 66: +  -mfn_view             - print information about the solver used
 67: .  -mfn_view_pre         - print information about the solver before the solve starts
 68: .  -mfn_view_mat         - view the matrix
 69: .  -mfn_view_rhs         - view right hand side vector
 70: .  -mfn_view_solution    - view computed solution vector
 71: -  -mfn_converged_reason - print reason for convergence/divergence, and number of iterations

 73:    Notes:
 74:    The matrix $A$ is specified with `MFNSetOperator()`. The function $f$ is
 75:    specified via the `FN` object obtained with `MFNGetFN()` or set with `MFNSetFN()`.

 77:    All the command-line options listed above admit an optional argument specifying
 78:    the viewer type and options. For instance, use `-mfn_view_mat binary:amatrix.bin`
 79:    to save the matrix to a binary file, or `-mfn_view_solution :sol.m:ascii_matlab`
 80:    to save the solution in a file that can be executed in Matlab.

 82:    Level: beginner

 84: .seealso: [](ch:mfn), `MFNCreate()`, `MFNSetUp()`, `MFNDestroy()`, `MFNSetTolerances()`, `MFNSetOperator()`, `MFNSetFN()`
 85: @*/
 86: PetscErrorCode MFNSolve(MFN mfn,Vec b,Vec x)
 87: {
 88:   PetscFunctionBegin;
 91:   PetscCheckSameComm(mfn,1,b,2);
 93:   if (b!=x) PetscCheckSameComm(mfn,1,x,3);
 94:   mfn->transpose_solve = PETSC_FALSE;
 95:   PetscCall(MFNSolve_Private(mfn,b,x));
 96:   PetscFunctionReturn(PETSC_SUCCESS);
 97: }

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

103:    Collective

105:    Input Parameters:
106: +  mfn - the matrix function solver context
107: -  b   - the right hand side vector

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

113:    Notes:
114:    The matrix $A$ is specified with `MFNSetOperator()`. The function $f$ is
115:    specified via the `FN` object obtained with `MFNGetFN()` or set with `MFNSetFN()`.

117:    See available command-line options at `MFNSolve()`.

119:    Developer Notes:
120:    This is currently implemented with an explicit transpose matrix created
121:    with `MatCreateTranspose()`.

123:    Currently there is no conjugate-transpose version.

125:    Level: beginner

127: .seealso: [](ch:mfn), `MFNSolve()`, `MatCreateTranspose()`
128: @*/
129: PetscErrorCode MFNSolveTranspose(MFN mfn,Vec b,Vec x)
130: {
131:   PetscFunctionBegin;
134:   PetscCheckSameComm(mfn,1,b,2);
136:   if (b!=x) PetscCheckSameComm(mfn,1,x,3);
137:   mfn->transpose_solve = PETSC_TRUE;
138:   if (!mfn->AT) PetscCall(MatCreateTranspose(mfn->A,&mfn->AT));
139:   PetscCall(MFNSolve_Private(mfn,b,x));
140:   PetscFunctionReturn(PETSC_SUCCESS);
141: }

143: /*@
144:    MFNGetIterationNumber - Gets the current iteration number. If the
145:    call to `MFNSolve()` is complete, then it returns the number of iterations
146:    carried out by the solution method.

148:    Not Collective

150:    Input Parameter:
151: .  mfn - the matrix function solver context

153:    Output Parameter:
154: .  its - number of iterations

156:    Note:
157:    During the $i$-th iteration this call returns $i-1$. If `MFNSolve()` is
158:    complete, then parameter `its` contains either the iteration number at
159:    which convergence was successfully reached, or failure was detected.
160:    Call `MFNGetConvergedReason()` to determine if the solver converged or
161:    failed and why.

163:    Level: intermediate

165: .seealso: [](ch:mfn), `MFNGetConvergedReason()`, `MFNSetTolerances()`
166: @*/
167: PetscErrorCode MFNGetIterationNumber(MFN mfn,PetscInt *its)
168: {
169:   PetscFunctionBegin;
171:   PetscAssertPointer(its,2);
172:   *its = mfn->its;
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

176: /*@
177:    MFNGetConvergedReason - Gets the reason why the `MFNSolve()` iteration was
178:    stopped.

180:    Not Collective

182:    Input Parameter:
183: .  mfn - the matrix function solver context

185:    Output Parameter:
186: .  reason - negative value indicates diverged, positive value converged, see
187:    `MFNConvergedReason` for the possible values

189:    Options Database Key:
190: .  -mfn_converged_reason - print reason for convergence/divergence, and number of iterations

192:    Notes:
193:    If this routine is called before or doing the `MFNSolve()` the value of
194:    `MFN_CONVERGED_ITERATING` is returned.

196:    Basic solvers (e.g., unrestarted Krylov iterations) cannot determine if the
197:    computation is accurate up to the requested tolerance. In that case, the
198:    converged reason is set to `MFN_CONVERGED_ITS` if the requested number of steps
199:    (for instance, the `ncv` value in unrestarted Krylov methods) have been
200:    completed successfully.

202:    Level: intermediate

204: .seealso: [](ch:mfn), `MFNSetTolerances()`, `MFNSolve()`, `MFNConvergedReason`, `MFNSetErrorIfNotConverged()`
205: @*/
206: PetscErrorCode MFNGetConvergedReason(MFN mfn,MFNConvergedReason *reason)
207: {
208:   PetscFunctionBegin;
210:   PetscAssertPointer(reason,2);
211:   *reason = mfn->reason;
212:   PetscFunctionReturn(PETSC_SUCCESS);
213: }