Actual source code: nepsolve.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:    NEP routines related to the solution process

     13:    References:

     15:        [1] C. Campos and J.E. Roman, "NEP: a module for the parallel solution
     16:            of nonlinear eigenvalue problems in SLEPc", ACM Trans. Math. Soft.
     17:            47(3), 23:1--23:29, 2021.
     18: */

     20: #include <slepc/private/nepimpl.h>
     21: #include <slepc/private/bvimpl.h>
     22: #include <petscdraw.h>

     24: static PetscBool  cited = PETSC_FALSE;
     25: static const char citation[] =
     26:   "@Article{slepc-nep,\n"
     27:   "   author = \"C. Campos and J. E. Roman\",\n"
     28:   "   title = \"{NEP}: a module for the parallel solution of nonlinear eigenvalue problems in {SLEPc}\",\n"
     29:   "   journal = \"{ACM} Trans. Math. Software\",\n"
     30:   "   volume = \"47\",\n"
     31:   "   number = \"3\",\n"
     32:   "   pages = \"23:1--23:29\",\n"
     33:   "   year = \"2021\",\n"
     34:   "   doi = \"10.1145/3447544\"\n"
     35:   "}\n";

     37: PetscErrorCode NEPComputeVectors(NEP nep)
     38: {
     39:   PetscFunctionBegin;
     40:   NEPCheckSolved(nep,1);
     41:   if (nep->state==NEP_STATE_SOLVED) PetscTryTypeMethod(nep,computevectors);
     42:   nep->state = NEP_STATE_EIGENVECTORS;
     43:   PetscFunctionReturn(PETSC_SUCCESS);
     44: }

     46: /*@
     47:    NEPSolve - Solves the nonlinear eigenproblem.

     49:    Collective

     51:    Input Parameter:
     52: .  nep - the nonlinear eigensolver context

     54:    Options Database Keys:
     55: +  -nep_view             - print information about the solver once the solve is complete
     56: .  -nep_view_pre         - print information about the solver before the solve starts
     57: .  -nep_view_matk        - view the split form matrix $A_k$ (replace `k` by an integer from 0 to `nt`-1)
     58: .  -nep_view_fnk         - view the split form function $f_k$ (replace `k` by an integer from 0 to `nt`-1)
     59: .  -nep_view_vectors     - view the computed eigenvectors
     60: .  -nep_view_values      - view the computed eigenvalues
     61: .  -nep_converged_reason - print reason for convergence/divergence, and number of iterations
     62: .  -nep_error_absolute   - print absolute errors of each eigenpair
     63: .  -nep_error_relative   - print relative errors of each eigenpair
     64: -  -nep_error_backward   - print backward errors of each eigenpair

     66:    Notes:
     67:    `NEPSolve()` will return without generating an error regardless of whether
     68:    all requested solutions were computed or not. Call `NEPGetConverged()` to get the
     69:    actual number of computed solutions, and `NEPGetConvergedReason()` to determine if
     70:    the solver converged or failed and why.

     72:    All the command-line options listed above admit an optional argument specifying
     73:    the viewer type and options. For instance, use `-nep_view_vectors binary:myvecs.bin`
     74:    to save the eigenvectors to a binary file, `-nep_view_values draw` to draw the computed
     75:    eigenvalues graphically, or `-nep_error_relative :myerr.m:ascii_matlab` to save
     76:    the errors in a file that can be executed in Matlab.
     77:    See `PetscObjectViewFromOptions()` for more details.

     79:    Level: beginner

     81: .seealso: [](ch:nep), `NEPCreate()`, `NEPSetUp()`, `NEPDestroy()`, `NEPSetTolerances()`, `NEPGetConverged()`, `NEPGetConvergedReason()`
     82: @*/
     83: PetscErrorCode NEPSolve(NEP nep)
     84: {
     85:   PetscInt       i;
     86:   char           str[16];

     88:   PetscFunctionBegin;
     90:   if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
     91:   PetscCall(PetscCitationsRegister(citation,&cited));
     92:   PetscCall(PetscLogEventBegin(NEP_Solve,nep,0,0,0));

     94:   /* call setup */
     95:   PetscCall(NEPSetUp(nep));
     96:   nep->nconv = 0;
     97:   nep->its = 0;
     98:   for (i=0;i<nep->ncv;i++) {
     99:     nep->eigr[i]   = 0.0;
    100:     nep->eigi[i]   = 0.0;
    101:     nep->errest[i] = 0.0;
    102:     nep->perm[i]   = i;
    103:   }
    104:   PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view_pre"));
    105:   PetscCall(RGViewFromOptions(nep->rg,NULL,"-rg_view"));

    107:   /* call solver */
    108:   PetscUseTypeMethod(nep,solve);
    109:   PetscCheck(nep->reason,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
    110:   nep->state = NEP_STATE_SOLVED;

    112:   /* Only the first nconv columns contain useful information */
    113:   PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
    114:   if (nep->twosided) PetscCall(BVSetActiveColumns(nep->W,0,nep->nconv));

    116:   if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
    117:     PetscCall(NEPComputeVectors(nep));
    118:     PetscCall(NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv));
    119:     nep->state = NEP_STATE_EIGENVECTORS;
    120:   }

    122:   /* sort eigenvalues according to nep->which parameter */
    123:   PetscCall(SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm));
    124:   PetscCall(PetscLogEventEnd(NEP_Solve,nep,0,0,0));

    126:   /* various viewers */
    127:   PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
    128:   PetscCall(NEPConvergedReasonViewFromOptions(nep));
    129:   PetscCall(NEPErrorViewFromOptions(nep));
    130:   PetscCall(NEPValuesViewFromOptions(nep));
    131:   PetscCall(NEPVectorsViewFromOptions(nep));
    132:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
    133:     for (i=0;i<nep->nt;i++) {
    134:       PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_mat%" PetscInt_FMT,i));
    135:       PetscCall(MatViewFromOptions(nep->A[i],(PetscObject)nep,str));
    136:       PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_fn%" PetscInt_FMT,i));
    137:       PetscCall(FNViewFromOptions(nep->f[i],(PetscObject)nep,str));
    138:     }
    139:   }

    141:   /* Remove the initial subspace */
    142:   nep->nini = 0;

    144:   /* Reset resolvent information */
    145:   PetscCall(MatDestroy(&nep->resolvent));
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

    149: /*@
    150:    NEPProjectOperator - Computes the projection of the nonlinear operator.

    152:    Collective

    154:    Input Parameters:
    155: +  nep - the nonlinear eigensolver context
    156: .  j0  - initial index
    157: -  j1  - final index

    159:    Notes:
    160:    This is available for split operator only.

    162:    The nonlinear operator $T(\lambda)$ is projected onto $\operatorname{span}(V)$, where $V$ is
    163:    an orthonormal basis built internally by the solver. The projected
    164:    operator is equal to $\sum_i V^* A_i V f_i(\lambda)$, so this function
    165:    computes all matrices $E_i = V^* A_i V$, and stores them in the extra
    166:    matrices inside `DS`, see `DSMatType`. Only rows/columns in the range [`j0`,`j1`-1] are computed,
    167:    the previous ones are assumed to be available already.

    169:    Level: developer

    171: .seealso: [](ch:nep), `NEPSetSplitOperator()`
    172: @*/
    173: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
    174: {
    175:   PetscInt       k;
    176:   Mat            G;

    178:   PetscFunctionBegin;
    182:   NEPCheckProblem(nep,1);
    183:   NEPCheckSplit(nep,1);
    184:   PetscCall(BVSetActiveColumns(nep->V,j0,j1));
    185:   for (k=0;k<nep->nt;k++) {
    186:     PetscCall(DSGetMat(nep->ds,DSMatExtra[k],&G));
    187:     PetscCall(BVMatProject(nep->V,nep->A[k],nep->V,G));
    188:     PetscCall(DSRestoreMat(nep->ds,DSMatExtra[k],&G));
    189:   }
    190:   PetscFunctionReturn(PETSC_SUCCESS);
    191: }

    193: /*@
    194:    NEPApplyFunction - Applies the nonlinear function $T(\lambda)$ to a given vector.

    196:    Collective

    198:    Input Parameters:
    199: +  nep    - the nonlinear eigensolver context
    200: .  lambda - scalar argument
    201: .  x      - vector to be multiplied against
    202: -  v      - workspace vector (used only in the case of split form)

    204:    Output Parameters:
    205: +  y   - result vector
    206: .  A   - (optional) Function matrix, for callback interface only
    207: -  B   - (unused) preconditioning matrix

    209:    Note:
    210:    If the nonlinear operator is represented in split form, the result
    211:    $y = T(\lambda)x$ is computed without building $T(\lambda)$ explicitly. In
    212:    that case, parameters `A` and `B` are not used. Otherwise, the matrix
    213:    $T(\lambda)$ is built and the effect is the same as a call to
    214:    `NEPComputeFunction()` followed by a `MatMult()`.

    216:    Level: developer

    218: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeFunction()`, `NEPApplyAdjoint()`
    219: @*/
    220: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
    221: {
    222:   PetscInt       i;
    223:   PetscScalar    alpha;

    225:   PetscFunctionBegin;

    234:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
    235:     PetscCall(VecSet(y,0.0));
    236:     for (i=0;i<nep->nt;i++) {
    237:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
    238:       PetscCall(MatMult(nep->A[i],x,v));
    239:       PetscCall(VecAXPY(y,alpha,v));
    240:     }
    241:   } else {
    242:     if (!A) A = nep->function;
    243:     PetscCall(NEPComputeFunction(nep,lambda,A,A));
    244:     PetscCall(MatMult(A,x,y));
    245:   }
    246:   PetscFunctionReturn(PETSC_SUCCESS);
    247: }

    249: /*@
    250:    NEPApplyAdjoint - Applies the adjoint nonlinear function $T^*(\lambda)$ to a given vector.

    252:    Collective

    254:    Input Parameters:
    255: +  nep    - the nonlinear eigensolver context
    256: .  lambda - scalar argument
    257: .  x      - vector to be multiplied against
    258: -  v      - workspace vector (used only in the case of split form)

    260:    Output Parameters:
    261: +  y   - result vector
    262: .  A   - (optional) Function matrix, for callback interface only
    263: -  B   - (unused) preconditioning matrix

    265:    Note:
    266:    If the nonlinear operator is represented in split form, the result
    267:    $y = T^*(\lambda)x$ is computed without building $T(\lambda)$ explicitly. In
    268:    that case, parameters `A` and `B` are not used. Otherwise, the matrix
    269:    $T(\lambda)$ is built and the effect is the same as a call to
    270:    `NEPComputeFunction()` followed by a `MatMultHermitianTranspose()`.

    272:    Level: developer

    274: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeFunction()`, `NEPApplyFunction()`
    275: @*/
    276: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
    277: {
    278:   PetscInt       i;
    279:   PetscScalar    alpha;
    280:   Vec            w;

    282:   PetscFunctionBegin;

    291:   PetscCall(VecDuplicate(x,&w));
    292:   PetscCall(VecCopy(x,w));
    293:   PetscCall(VecConjugate(w));
    294:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
    295:     PetscCall(VecSet(y,0.0));
    296:     for (i=0;i<nep->nt;i++) {
    297:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
    298:       PetscCall(MatMultHermitianTranspose(nep->A[i],w,v));
    299:       PetscCall(VecAXPY(y,alpha,v));
    300:     }
    301:   } else {
    302:     if (!A) A = nep->function;
    303:     PetscCall(NEPComputeFunction(nep,lambda,A,A));
    304:     PetscCall(MatMultHermitianTranspose(A,w,y));
    305:   }
    306:   PetscCall(VecDestroy(&w));
    307:   PetscCall(VecConjugate(y));
    308:   PetscFunctionReturn(PETSC_SUCCESS);
    309: }

    311: /*@
    312:    NEPApplyJacobian - Applies the nonlinear Jacobian $T'(\lambda)$ to a given vector.

    314:    Collective

    316:    Input Parameters:
    317: +  nep    - the nonlinear eigensolver context
    318: .  lambda - scalar argument
    319: .  x      - vector to be multiplied against
    320: -  v      - workspace vector (used only in the case of split form)

    322:    Output Parameters:
    323: +  y   - result vector
    324: -  A   - (optional) Jacobian matrix, for callback interface only

    326:    Note:
    327:    If the nonlinear operator is represented in split form, the result
    328:    $y = T'(\lambda)x$ is computed without building $T'(\lambda)$ explicitly. In
    329:    that case, parameter `A` is not used. Otherwise, the matrix
    330:    $T'(\lambda)$ is built and the effect is the same as a call to
    331:    `NEPComputeJacobian()` followed by a `MatMult()`.

    333:    Level: developer

    335: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeJacobian()`
    336: @*/
    337: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
    338: {
    339:   PetscInt       i;
    340:   PetscScalar    alpha;

    342:   PetscFunctionBegin;

    350:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
    351:     PetscCall(VecSet(y,0.0));
    352:     for (i=0;i<nep->nt;i++) {
    353:       PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
    354:       PetscCall(MatMult(nep->A[i],x,v));
    355:       PetscCall(VecAXPY(y,alpha,v));
    356:     }
    357:   } else {
    358:     if (!A) A = nep->jacobian;
    359:     PetscCall(NEPComputeJacobian(nep,lambda,A));
    360:     PetscCall(MatMult(A,x,y));
    361:   }
    362:   PetscFunctionReturn(PETSC_SUCCESS);
    363: }

    365: /*@
    366:    NEPGetIterationNumber - Gets the current iteration number. If the
    367:    call to `NEPSolve()` is complete, then it returns the number of iterations
    368:    carried out by the solution method.

    370:    Not Collective

    372:    Input Parameter:
    373: .  nep - the nonlinear eigensolver context

    375:    Output Parameter:
    376: .  its - number of iterations

    378:    Note:
    379:    During the $i$-th iteration this call returns $i-1$. If `NEPSolve()` is
    380:    complete, then parameter `its` contains either the iteration number at
    381:    which convergence was successfully reached, or failure was detected.
    382:    Call `NEPGetConvergedReason()` to determine if the solver converged or
    383:    failed and why.

    385:    Level: intermediate

    387: .seealso: [](ch:nep), `NEPGetConvergedReason()`, `NEPSetTolerances()`
    388: @*/
    389: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
    390: {
    391:   PetscFunctionBegin;
    393:   PetscAssertPointer(its,2);
    394:   *its = nep->its;
    395:   PetscFunctionReturn(PETSC_SUCCESS);
    396: }

    398: /*@
    399:    NEPGetConverged - Gets the number of converged eigenpairs.

    401:    Not Collective

    403:    Input Parameter:
    404: .  nep - the nonlinear eigensolver context

    406:    Output Parameter:
    407: .  nconv - number of converged eigenpairs

    409:    Note:
    410:    This function should be called after `NEPSolve()` has finished.

    412:    The value `nconv` may be different from the number of requested solutions
    413:    `nev`, but not larger than `ncv`, see `NEPSetDimensions()`.

    415:    Level: beginner

    417: .seealso: [](ch:nep), `NEPSetDimensions()`, `NEPSolve()`, `NEPGetEigenpair()`
    418: @*/
    419: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
    420: {
    421:   PetscFunctionBegin;
    423:   PetscAssertPointer(nconv,2);
    424:   NEPCheckSolved(nep,1);
    425:   *nconv = nep->nconv;
    426:   PetscFunctionReturn(PETSC_SUCCESS);
    427: }

    429: /*@
    430:    NEPGetConvergedReason - Gets the reason why the `NEPSolve()` iteration was
    431:    stopped.

    433:    Not Collective

    435:    Input Parameter:
    436: .  nep - the nonlinear eigensolver context

    438:    Output Parameter:
    439: .  reason - negative value indicates diverged, positive value converged, see
    440:    `NEPConvergedReason` for the possible values

    442:    Options Database Key:
    443: .  -nep_converged_reason - print reason for convergence/divergence, and number of iterations

    445:    Note:
    446:    If this routine is called before or doing the `NEPSolve()` the value of
    447:    `NEP_CONVERGED_ITERATING` is returned.

    449:    Level: intermediate

    451: .seealso: [](ch:nep), `NEPSetTolerances()`, `NEPSolve()`, `NEPConvergedReason`
    452: @*/
    453: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
    454: {
    455:   PetscFunctionBegin;
    457:   PetscAssertPointer(reason,2);
    458:   NEPCheckSolved(nep,1);
    459:   *reason = nep->reason;
    460:   PetscFunctionReturn(PETSC_SUCCESS);
    461: }

    463: /*@
    464:    NEPGetEigenpair - Gets the `i`-th solution of the eigenproblem as computed by
    465:    `NEPSolve()`. The solution consists in both the eigenvalue and the eigenvector.

    467:    Collective

    469:    Input Parameters:
    470: +  nep - the nonlinear eigensolver context
    471: -  i   - index of the solution

    473:    Output Parameters:
    474: +  eigr - real part of eigenvalue
    475: .  eigi - imaginary part of eigenvalue
    476: .  Vr   - real part of eigenvector
    477: -  Vi   - imaginary part of eigenvector

    479:    Notes:
    480:    It is allowed to pass `NULL` for `Vr` and `Vi`, if the eigenvector is not
    481:    required. Otherwise, the caller must provide valid `Vec` objects, i.e.,
    482:    they must be created by the calling program with e.g. `MatCreateVecs()`.

    484:    If the eigenvalue is real, then `eigi` and `Vi` are set to zero. If PETSc is
    485:    configured with complex scalars the eigenvalue is stored
    486:    directly in `eigr` (`eigi` is set to zero) and the eigenvector in `Vr` (`Vi`
    487:    is set to zero). In any case, the user can pass `NULL` in any argument that
    488:    is not required.

    490:    The index `i` should be a value between 0 and `nconv`-1 (see `NEPGetConverged()`).
    491:    Eigenpairs are indexed according to the ordering criterion established
    492:    with `NEPSetWhichEigenpairs()`.

    494:    The eigenvector is normalized to have unit norm.

    496:    Level: beginner

    498: .seealso: [](ch:nep), `NEPSolve()`, `NEPGetConverged()`, `NEPSetWhichEigenpairs()`, `NEPGetLeftEigenvector()`
    499: @*/
    500: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
    501: {
    502:   PetscInt       k;

    504:   PetscFunctionBegin;
    509:   NEPCheckSolved(nep,1);
    510:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    511:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");

    513:   PetscCall(NEPComputeVectors(nep));
    514:   k = nep->perm[i];

    516:   /* eigenvalue */
    517: #if defined(PETSC_USE_COMPLEX)
    518:   if (eigr) *eigr = nep->eigr[k];
    519:   if (eigi) *eigi = 0;
    520: #else
    521:   if (eigr) *eigr = nep->eigr[k];
    522:   if (eigi) *eigi = nep->eigi[k];
    523: #endif

    525:   /* eigenvector */
    526:   PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
    527:   PetscFunctionReturn(PETSC_SUCCESS);
    528: }

    530: /*@
    531:    NEPGetLeftEigenvector - Gets the `i`-th left eigenvector as computed by `NEPSolve()`.

    533:    Collective

    535:    Input Parameters:
    536: +  nep - the nonlinear eigensolver context
    537: -  i   - index of the solution

    539:    Output Parameters:
    540: +  Wr   - real part of left eigenvector
    541: -  Wi   - imaginary part of left eigenvector

    543:    Notes:
    544:    The caller must provide valid `Vec` objects, i.e., they must be created
    545:    by the calling program with e.g. `MatCreateVecs()`.

    547:    If the corresponding eigenvalue is real, then `Wi` is set to zero. If PETSc is
    548:    configured with complex scalars the eigenvector is stored directly in `Wr`
    549:    (`Wi` is set to zero). In any case, the user can pass `NULL` in `Wr` or `Wi` if
    550:    one of them is not required.

    552:    The index `i` should be a value between 0 and `nconv`-1 (see `NEPGetConverged()`).
    553:    Eigensolutions are indexed according to the ordering criterion established
    554:    with `NEPSetWhichEigenpairs()`.

    556:    Left eigenvectors are available only if the twosided flag was set, see
    557:    `NEPSetTwoSided()`.

    559:    Level: intermediate

    561: .seealso: [](ch:nep), `NEPGetEigenpair()`, `NEPGetConverged()`, `NEPSetWhichEigenpairs()`, `NEPSetTwoSided()`
    562: @*/
    563: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
    564: {
    565:   PetscInt       k;

    567:   PetscFunctionBegin;
    572:   NEPCheckSolved(nep,1);
    573:   PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
    574:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    575:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
    576:   PetscCall(NEPComputeVectors(nep));
    577:   k = nep->perm[i];
    578:   PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
    579:   PetscFunctionReturn(PETSC_SUCCESS);
    580: }

    582: /*@
    583:    NEPGetErrorEstimate - Returns the error estimate associated to the `i`-th
    584:    computed eigenpair.

    586:    Not Collective

    588:    Input Parameters:
    589: +  nep - the nonlinear eigensolver context
    590: -  i   - index of eigenpair

    592:    Output Parameter:
    593: .  errest - the error estimate

    595:    Note:
    596:    This is the error estimate used internally by the eigensolver. The actual
    597:    error bound can be computed with `NEPComputeError()`.

    599:    Level: advanced

    601: .seealso: [](ch:nep), `NEPComputeError()`
    602: @*/
    603: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
    604: {
    605:   PetscFunctionBegin;
    607:   PetscAssertPointer(errest,3);
    608:   NEPCheckSolved(nep,1);
    609:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
    610:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
    611:   *errest = nep->errest[nep->perm[i]];
    612:   PetscFunctionReturn(PETSC_SUCCESS);
    613: }

    615: /*
    616:    NEPComputeResidualNorm_Private - Computes the norm of the residual vector
    617:    associated with an eigenpair.

    619:    Input Parameters:
    620:      adj    - whether the adjoint T^* must be used instead of T
    621:      lambda - eigenvalue
    622:      x      - eigenvector
    623:      w      - array of work vectors (two vectors in split form, one vector otherwise)
    624: */
    625: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
    626: {
    627:   Vec            y,z=NULL;

    629:   PetscFunctionBegin;
    630:   y = w[0];
    631:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
    632:   if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
    633:   else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
    634:   PetscCall(VecNorm(y,NORM_2,norm));
    635:   PetscFunctionReturn(PETSC_SUCCESS);
    636: }

    638: /*@
    639:    NEPComputeError - Computes the error (based on the residual norm) associated
    640:    with the `i`-th computed eigenpair.

    642:    Collective

    644:    Input Parameters:
    645: +  nep  - the nonlinear eigensolver context
    646: .  i    - the solution index
    647: -  type - the type of error to compute, see `NEPErrorType`

    649:    Output Parameter:
    650: .  error - the error

    652:    Notes:
    653:    The error can be computed in various ways, all of them based on the residual
    654:    norm computed as $\|T(\lambda)x\|_2$ where $(\lambda,x)$ is the approximate eigenpair.

    656:    If the computation of left eigenvectors was enabled with `NEPSetTwoSided()`,
    657:    then the error will be computed using the maximum of the value above and
    658:    the left residual norm $\|y^*T(\lambda)\|_2$, where $y$ is the approximate left
    659:    eigenvector.

    661:    Level: beginner

    663: .seealso: [](ch:nep), `NEPErrorType`, `NEPSolve()`, `NEPGetErrorEstimate()`, `NEPSetTwoSided()`
    664: @*/
    665: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
    666: {
    667:   Vec            xr,xi=NULL;
    668:   PetscInt       j,nwork,issplit=0;
    669:   PetscScalar    kr,ki,s;
    670:   PetscReal      er,z=0.0,errorl,nrm;
    671:   PetscBool      flg;

    673:   PetscFunctionBegin;
    677:   PetscAssertPointer(error,4);
    678:   NEPCheckSolved(nep,1);

    680:   /* allocate work vectors */
    681: #if defined(PETSC_USE_COMPLEX)
    682:   nwork = 2;
    683: #else
    684:   nwork = 3;
    685: #endif
    686:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
    687:     issplit = 1;
    688:     nwork++;  /* need an extra work vector for NEPComputeResidualNorm_Private */
    689:   }
    690:   PetscCall(NEPSetWorkVecs(nep,nwork));
    691:   xr = nep->work[issplit+1];
    692: #if !defined(PETSC_USE_COMPLEX)
    693:   xi = nep->work[issplit+2];
    694: #endif

    696:   /* compute residual norms */
    697:   PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
    698: #if !defined(PETSC_USE_COMPLEX)
    699:   PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
    700: #endif
    701:   PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
    702:   PetscCall(VecNorm(xr,NORM_2,&er));

    704:   /* if two-sided, compute left residual norm and take the maximum */
    705:   if (nep->twosided) {
    706:     PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
    707:     PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
    708:     *error = PetscMax(*error,errorl);
    709:   }

    711:   /* compute error */
    712:   switch (type) {
    713:     case NEP_ERROR_ABSOLUTE:
    714:       break;
    715:     case NEP_ERROR_RELATIVE:
    716:       *error /= PetscAbsScalar(kr)*er;
    717:       break;
    718:     case NEP_ERROR_BACKWARD:
    719:       if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
    720:         PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
    721:         PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
    722:         PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
    723:         PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
    724:         *error /= nrm*er;
    725:         break;
    726:       }
    727:       /* initialization of matrix norms */
    728:       if (!nep->nrma[0]) {
    729:         for (j=0;j<nep->nt;j++) {
    730:           PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
    731:           PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
    732:           PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
    733:         }
    734:       }
    735:       for (j=0;j<nep->nt;j++) {
    736:         PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
    737:         z = z + nep->nrma[j]*PetscAbsScalar(s);
    738:       }
    739:       *error /= z*er;
    740:       break;
    741:     default:
    742:       SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
    743:   }
    744:   PetscFunctionReturn(PETSC_SUCCESS);
    745: }

    747: /*@
    748:    NEPComputeFunction - Computes the function matrix $T(\lambda)$ that has been
    749:    set with `NEPSetFunction()`.

    751:    Collective

    753:    Input Parameters:
    754: +  nep    - the nonlinear eigensolver context
    755: -  lambda - the scalar argument

    757:    Output Parameters:
    758: +  A   - Function matrix
    759: -  B   - optional preconditioning matrix

    761:    Note:
    762:    `NEPComputeFunction()` is typically used within nonlinear eigensolvers
    763:    implementations, so most users would not generally call this routine
    764:    themselves.

    766:    Level: developer

    768: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetFunction()`, `NEPComputeJacobian()`
    769: @*/
    770: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
    771: {
    772:   PetscInt       i;
    773:   PetscScalar    alpha;

    775:   PetscFunctionBegin;
    777:   NEPCheckProblem(nep,1);
    778:   switch (nep->fui) {
    779:   case NEP_USER_INTERFACE_CALLBACK:
    780:     PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
    781:     PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
    782:     PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
    783:     PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
    784:     break;
    785:   case NEP_USER_INTERFACE_SPLIT:
    786:     PetscCall(MatZeroEntries(A));
    787:     if (A != B) PetscCall(MatZeroEntries(B));
    788:     for (i=0;i<nep->nt;i++) {
    789:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
    790:       PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
    791:       if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
    792:     }
    793:     break;
    794:   }
    795:   PetscFunctionReturn(PETSC_SUCCESS);
    796: }

    798: /*@
    799:    NEPComputeJacobian - Computes the Jacobian matrix $T'(\lambda)$ that has been
    800:    set with `NEPSetJacobian()`.

    802:    Collective

    804:    Input Parameters:
    805: +  nep    - the nonlinear eigensolver context
    806: -  lambda - the scalar argument

    808:    Output Parameter:
    809: .  A   - Jacobian matrix

    811:    Notes:
    812:    Most users should not need to explicitly call this routine, as it
    813:    is used internally within the nonlinear eigensolvers.

    815:    Level: developer

    817: .seealso: [](ch:nep), `NEPSetJacobian()`, `NEPGetJacobian()`, `NEPComputeFunction()`
    818: @*/
    819: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
    820: {
    821:   PetscInt       i;
    822:   PetscScalar    alpha;

    824:   PetscFunctionBegin;
    826:   NEPCheckProblem(nep,1);
    827:   switch (nep->fui) {
    828:   case NEP_USER_INTERFACE_CALLBACK:
    829:     PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
    830:     PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
    831:     PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
    832:     PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
    833:     break;
    834:   case NEP_USER_INTERFACE_SPLIT:
    835:     PetscCall(MatZeroEntries(A));
    836:     for (i=0;i<nep->nt;i++) {
    837:       PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
    838:       PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
    839:     }
    840:     break;
    841:   }
    842:   PetscFunctionReturn(PETSC_SUCCESS);
    843: }