Actual source code: mfnsolve.c

      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.
     81:    See `PetscObjectViewFromOptions()` for more details.

     83:    Level: beginner

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

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

    104:    Collective

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

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

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

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

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

    124:    Currently there is no conjugate-transpose version.

    126:    Level: beginner

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

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

    149:    Not Collective

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

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

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

    164:    Level: intermediate

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

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

    181:    Not Collective

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

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

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

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

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

    203:    Level: intermediate

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