Actual source code: nepsolve.c

slepc-3.20.1 2023-11-27
Report Typos and Errors
  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 eigensystem.

 49:    Collective

 51:    Input Parameter:
 52: .  nep - eigensolver context obtained from NEPCreate()

 54:    Options Database Keys:
 55: +  -nep_view - print information about the solver used
 56: .  -nep_view_vectors - view the computed eigenvectors
 57: .  -nep_view_values - view the computed eigenvalues
 58: .  -nep_converged_reason - print reason for convergence, and number of iterations
 59: .  -nep_error_absolute - print absolute errors of each eigenpair
 60: .  -nep_error_relative - print relative errors of each eigenpair
 61: -  -nep_error_backward - print backward errors of each eigenpair

 63:    Notes:
 64:    All the command-line options listed above admit an optional argument specifying
 65:    the viewer type and options. For instance, use '-nep_view_vectors binary:myvecs.bin'
 66:    to save the eigenvectors to a binary file, '-nep_view_values draw' to draw the computed
 67:    eigenvalues graphically, or '-nep_error_relative :myerr.m:ascii_matlab' to save
 68:    the errors in a file that can be executed in Matlab.

 70:    Level: beginner

 72: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
 73: @*/
 74: PetscErrorCode NEPSolve(NEP nep)
 75: {
 76:   PetscInt       i;

 78:   PetscFunctionBegin;
 80:   if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
 81:   PetscCall(PetscCitationsRegister(citation,&cited));
 82:   PetscCall(PetscLogEventBegin(NEP_Solve,nep,0,0,0));

 84:   /* call setup */
 85:   PetscCall(NEPSetUp(nep));
 86:   nep->nconv = 0;
 87:   nep->its = 0;
 88:   for (i=0;i<nep->ncv;i++) {
 89:     nep->eigr[i]   = 0.0;
 90:     nep->eigi[i]   = 0.0;
 91:     nep->errest[i] = 0.0;
 92:     nep->perm[i]   = i;
 93:   }
 94:   PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view_pre"));
 95:   PetscCall(RGViewFromOptions(nep->rg,NULL,"-rg_view"));

 97:   /* call solver */
 98:   PetscUseTypeMethod(nep,solve);
 99:   PetscCheck(nep->reason,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
100:   nep->state = NEP_STATE_SOLVED;

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

106:   if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
107:     PetscCall(NEPComputeVectors(nep));
108:     PetscCall(NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv));
109:     nep->state = NEP_STATE_EIGENVECTORS;
110:   }

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

116:   /* various viewers */
117:   PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
118:   PetscCall(NEPConvergedReasonViewFromOptions(nep));
119:   PetscCall(NEPErrorViewFromOptions(nep));
120:   PetscCall(NEPValuesViewFromOptions(nep));
121:   PetscCall(NEPVectorsViewFromOptions(nep));

123:   /* Remove the initial subspace */
124:   nep->nini = 0;

126:   /* Reset resolvent information */
127:   PetscCall(MatDestroy(&nep->resolvent));
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: /*@
132:    NEPProjectOperator - Computes the projection of the nonlinear operator.

134:    Collective

136:    Input Parameters:
137: +  nep - the nonlinear eigensolver context
138: .  j0  - initial index
139: -  j1  - final index

141:    Notes:
142:    This is available for split operator only.

144:    The nonlinear operator T(lambda) is projected onto span(V), where V is
145:    an orthonormal basis built internally by the solver. The projected
146:    operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
147:    computes all matrices Ei = V'*A_i*V, and stores them in the extra
148:    matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
149:    the previous ones are assumed to be available already.

151:    Level: developer

153: .seealso: NEPSetSplitOperator()
154: @*/
155: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
156: {
157:   PetscInt       k;
158:   Mat            G;

160:   PetscFunctionBegin;
164:   NEPCheckProblem(nep,1);
165:   NEPCheckSplit(nep,1);
166:   PetscCall(BVSetActiveColumns(nep->V,j0,j1));
167:   for (k=0;k<nep->nt;k++) {
168:     PetscCall(DSGetMat(nep->ds,DSMatExtra[k],&G));
169:     PetscCall(BVMatProject(nep->V,nep->A[k],nep->V,G));
170:     PetscCall(DSRestoreMat(nep->ds,DSMatExtra[k],&G));
171:   }
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: /*@
176:    NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.

178:    Collective

180:    Input Parameters:
181: +  nep    - the nonlinear eigensolver context
182: .  lambda - scalar argument
183: .  x      - vector to be multiplied against
184: -  v      - workspace vector (used only in the case of split form)

186:    Output Parameters:
187: +  y   - result vector
188: .  A   - (optional) Function matrix, for callback interface only
189: -  B   - (unused) preconditioning matrix

191:    Note:
192:    If the nonlinear operator is represented in split form, the result
193:    y = T(lambda)*x is computed without building T(lambda) explicitly. In
194:    that case, parameters A and B are not used. Otherwise, the matrix
195:    T(lambda) is built and the effect is the same as a call to
196:    NEPComputeFunction() followed by a MatMult().

198:    Level: developer

200: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
201: @*/
202: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
203: {
204:   PetscInt       i;
205:   PetscScalar    alpha;

207:   PetscFunctionBegin;

216:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
217:     PetscCall(VecSet(y,0.0));
218:     for (i=0;i<nep->nt;i++) {
219:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
220:       PetscCall(MatMult(nep->A[i],x,v));
221:       PetscCall(VecAXPY(y,alpha,v));
222:     }
223:   } else {
224:     if (!A) A = nep->function;
225:     PetscCall(NEPComputeFunction(nep,lambda,A,A));
226:     PetscCall(MatMult(A,x,y));
227:   }
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

231: /*@
232:    NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.

234:    Collective

236:    Input Parameters:
237: +  nep    - the nonlinear eigensolver context
238: .  lambda - scalar argument
239: .  x      - vector to be multiplied against
240: -  v      - workspace vector (used only in the case of split form)

242:    Output Parameters:
243: +  y   - result vector
244: .  A   - (optional) Function matrix, for callback interface only
245: -  B   - (unused) preconditioning matrix

247:    Level: developer

249: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
250: @*/
251: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
252: {
253:   PetscInt       i;
254:   PetscScalar    alpha;
255:   Vec            w;

257:   PetscFunctionBegin;

266:   PetscCall(VecDuplicate(x,&w));
267:   PetscCall(VecCopy(x,w));
268:   PetscCall(VecConjugate(w));
269:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
270:     PetscCall(VecSet(y,0.0));
271:     for (i=0;i<nep->nt;i++) {
272:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
273:       PetscCall(MatMultTranspose(nep->A[i],w,v));
274:       PetscCall(VecAXPY(y,alpha,v));
275:     }
276:   } else {
277:     if (!A) A = nep->function;
278:     PetscCall(NEPComputeFunction(nep,lambda,A,A));
279:     PetscCall(MatMultTranspose(A,w,y));
280:   }
281:   PetscCall(VecDestroy(&w));
282:   PetscCall(VecConjugate(y));
283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }

286: /*@
287:    NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.

289:    Collective

291:    Input Parameters:
292: +  nep    - the nonlinear eigensolver context
293: .  lambda - scalar argument
294: .  x      - vector to be multiplied against
295: -  v      - workspace vector (used only in the case of split form)

297:    Output Parameters:
298: +  y   - result vector
299: -  A   - (optional) Jacobian matrix, for callback interface only

301:    Note:
302:    If the nonlinear operator is represented in split form, the result
303:    y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
304:    that case, parameter A is not used. Otherwise, the matrix
305:    T'(lambda) is built and the effect is the same as a call to
306:    NEPComputeJacobian() followed by a MatMult().

308:    Level: developer

310: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
311: @*/
312: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
313: {
314:   PetscInt       i;
315:   PetscScalar    alpha;

317:   PetscFunctionBegin;

325:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
326:     PetscCall(VecSet(y,0.0));
327:     for (i=0;i<nep->nt;i++) {
328:       PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
329:       PetscCall(MatMult(nep->A[i],x,v));
330:       PetscCall(VecAXPY(y,alpha,v));
331:     }
332:   } else {
333:     if (!A) A = nep->jacobian;
334:     PetscCall(NEPComputeJacobian(nep,lambda,A));
335:     PetscCall(MatMult(A,x,y));
336:   }
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: /*@
341:    NEPGetIterationNumber - Gets the current iteration number. If the
342:    call to NEPSolve() is complete, then it returns the number of iterations
343:    carried out by the solution method.

345:    Not Collective

347:    Input Parameter:
348: .  nep - the nonlinear eigensolver context

350:    Output Parameter:
351: .  its - number of iterations

353:    Note:
354:    During the i-th iteration this call returns i-1. If NEPSolve() is
355:    complete, then parameter "its" contains either the iteration number at
356:    which convergence was successfully reached, or failure was detected.
357:    Call NEPGetConvergedReason() to determine if the solver converged or
358:    failed and why.

360:    Level: intermediate

362: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
363: @*/
364: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
365: {
366:   PetscFunctionBegin;
368:   PetscAssertPointer(its,2);
369:   *its = nep->its;
370:   PetscFunctionReturn(PETSC_SUCCESS);
371: }

373: /*@
374:    NEPGetConverged - Gets the number of converged eigenpairs.

376:    Not Collective

378:    Input Parameter:
379: .  nep - the nonlinear eigensolver context

381:    Output Parameter:
382: .  nconv - number of converged eigenpairs

384:    Note:
385:    This function should be called after NEPSolve() has finished.

387:    Level: beginner

389: .seealso: NEPSetDimensions(), NEPSolve(), NEPGetEigenpair()
390: @*/
391: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
392: {
393:   PetscFunctionBegin;
395:   PetscAssertPointer(nconv,2);
396:   NEPCheckSolved(nep,1);
397:   *nconv = nep->nconv;
398:   PetscFunctionReturn(PETSC_SUCCESS);
399: }

401: /*@
402:    NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
403:    stopped.

405:    Not Collective

407:    Input Parameter:
408: .  nep - the nonlinear eigensolver context

410:    Output Parameter:
411: .  reason - negative value indicates diverged, positive value converged

413:    Options Database Key:
414: .  -nep_converged_reason - print the reason to a viewer

416:    Notes:
417:    Possible values for reason are
418: +  NEP_CONVERGED_TOL - converged up to tolerance
419: .  NEP_CONVERGED_USER - converged due to a user-defined condition
420: .  NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
421: .  NEP_DIVERGED_BREAKDOWN - generic breakdown in method
422: .  NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
423: -  NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
424:    unrestarted solver

426:    Can only be called after the call to NEPSolve() is complete.

428:    Level: intermediate

430: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason
431: @*/
432: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
433: {
434:   PetscFunctionBegin;
436:   PetscAssertPointer(reason,2);
437:   NEPCheckSolved(nep,1);
438:   *reason = nep->reason;
439:   PetscFunctionReturn(PETSC_SUCCESS);
440: }

442: /*@C
443:    NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
444:    NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.

446:    Collective

448:    Input Parameters:
449: +  nep - nonlinear eigensolver context
450: -  i   - index of the solution

452:    Output Parameters:
453: +  eigr - real part of eigenvalue
454: .  eigi - imaginary part of eigenvalue
455: .  Vr   - real part of eigenvector
456: -  Vi   - imaginary part of eigenvector

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

463:    If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
464:    configured with complex scalars the eigenvalue is stored
465:    directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
466:    set to zero). In any case, the user can pass NULL in Vr or Vi if one of
467:    them is not required.

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

473:    Level: beginner

475: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
476: @*/
477: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
478: {
479:   PetscInt       k;

481:   PetscFunctionBegin;
486:   NEPCheckSolved(nep,1);
487:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
488:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");

490:   PetscCall(NEPComputeVectors(nep));
491:   k = nep->perm[i];

493:   /* eigenvalue */
494: #if defined(PETSC_USE_COMPLEX)
495:   if (eigr) *eigr = nep->eigr[k];
496:   if (eigi) *eigi = 0;
497: #else
498:   if (eigr) *eigr = nep->eigr[k];
499:   if (eigi) *eigi = nep->eigi[k];
500: #endif

502:   /* eigenvector */
503:   PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
504:   PetscFunctionReturn(PETSC_SUCCESS);
505: }

507: /*@
508:    NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().

510:    Collective

512:    Input Parameters:
513: +  nep - eigensolver context
514: -  i   - index of the solution

516:    Output Parameters:
517: +  Wr   - real part of left eigenvector
518: -  Wi   - imaginary part of left eigenvector

520:    Notes:
521:    The caller must provide valid Vec objects, i.e., they must be created
522:    by the calling program with e.g. MatCreateVecs().

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

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

533:    Left eigenvectors are available only if the twosided flag was set, see
534:    NEPSetTwoSided().

536:    Level: intermediate

538: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
539: @*/
540: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
541: {
542:   PetscInt       k;

544:   PetscFunctionBegin;
549:   NEPCheckSolved(nep,1);
550:   PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
551:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
552:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
553:   PetscCall(NEPComputeVectors(nep));
554:   k = nep->perm[i];
555:   PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
556:   PetscFunctionReturn(PETSC_SUCCESS);
557: }

559: /*@
560:    NEPGetErrorEstimate - Returns the error estimate associated to the i-th
561:    computed eigenpair.

563:    Not Collective

565:    Input Parameters:
566: +  nep - nonlinear eigensolver context
567: -  i   - index of eigenpair

569:    Output Parameter:
570: .  errest - the error estimate

572:    Notes:
573:    This is the error estimate used internally by the eigensolver. The actual
574:    error bound can be computed with NEPComputeError().

576:    Level: advanced

578: .seealso: NEPComputeError()
579: @*/
580: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
581: {
582:   PetscFunctionBegin;
584:   PetscAssertPointer(errest,3);
585:   NEPCheckSolved(nep,1);
586:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
587:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
588:   *errest = nep->errest[nep->perm[i]];
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

592: /*
593:    NEPComputeResidualNorm_Private - Computes the norm of the residual vector
594:    associated with an eigenpair.

596:    Input Parameters:
597:      adj    - whether the adjoint T^* must be used instead of T
598:      lambda - eigenvalue
599:      x      - eigenvector
600:      w      - array of work vectors (two vectors in split form, one vector otherwise)
601: */
602: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
603: {
604:   Vec            y,z=NULL;

606:   PetscFunctionBegin;
607:   y = w[0];
608:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
609:   if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
610:   else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
611:   PetscCall(VecNorm(y,NORM_2,norm));
612:   PetscFunctionReturn(PETSC_SUCCESS);
613: }

615: /*@
616:    NEPComputeError - Computes the error (based on the residual norm) associated
617:    with the i-th computed eigenpair.

619:    Collective

621:    Input Parameters:
622: +  nep  - the nonlinear eigensolver context
623: .  i    - the solution index
624: -  type - the type of error to compute

626:    Output Parameter:
627: .  error - the error

629:    Notes:
630:    The error can be computed in various ways, all of them based on the residual
631:    norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
632:    eigenvector.

634:    Level: beginner

636: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
637: @*/
638: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
639: {
640:   Vec            xr,xi=NULL;
641:   PetscInt       j,nwork,issplit=0;
642:   PetscScalar    kr,ki,s;
643:   PetscReal      er,z=0.0,errorl,nrm;
644:   PetscBool      flg;

646:   PetscFunctionBegin;
650:   PetscAssertPointer(error,4);
651:   NEPCheckSolved(nep,1);

653:   /* allocate work vectors */
654: #if defined(PETSC_USE_COMPLEX)
655:   nwork = 2;
656: #else
657:   nwork = 3;
658: #endif
659:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
660:     issplit = 1;
661:     nwork++;  /* need an extra work vector for NEPComputeResidualNorm_Private */
662:   }
663:   PetscCall(NEPSetWorkVecs(nep,nwork));
664:   xr = nep->work[issplit+1];
665: #if !defined(PETSC_USE_COMPLEX)
666:   xi = nep->work[issplit+2];
667: #endif

669:   /* compute residual norms */
670:   PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
671: #if !defined(PETSC_USE_COMPLEX)
672:   PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
673: #endif
674:   PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
675:   PetscCall(VecNorm(xr,NORM_2,&er));

677:   /* if two-sided, compute left residual norm and take the maximum */
678:   if (nep->twosided) {
679:     PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
680:     PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
681:     *error = PetscMax(*error,errorl);
682:   }

684:   /* compute error */
685:   switch (type) {
686:     case NEP_ERROR_ABSOLUTE:
687:       break;
688:     case NEP_ERROR_RELATIVE:
689:       *error /= PetscAbsScalar(kr)*er;
690:       break;
691:     case NEP_ERROR_BACKWARD:
692:       if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
693:         PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
694:         PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
695:         PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
696:         PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
697:         *error /= nrm*er;
698:         break;
699:       }
700:       /* initialization of matrix norms */
701:       if (!nep->nrma[0]) {
702:         for (j=0;j<nep->nt;j++) {
703:           PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
704:           PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
705:           PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
706:         }
707:       }
708:       for (j=0;j<nep->nt;j++) {
709:         PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
710:         z = z + nep->nrma[j]*PetscAbsScalar(s);
711:       }
712:       *error /= z*er;
713:       break;
714:     default:
715:       SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
716:   }
717:   PetscFunctionReturn(PETSC_SUCCESS);
718: }

720: /*@
721:    NEPComputeFunction - Computes the function matrix T(lambda) that has been
722:    set with NEPSetFunction().

724:    Collective

726:    Input Parameters:
727: +  nep    - the NEP context
728: -  lambda - the scalar argument

730:    Output Parameters:
731: +  A   - Function matrix
732: -  B   - optional preconditioning matrix

734:    Notes:
735:    NEPComputeFunction() is typically used within nonlinear eigensolvers
736:    implementations, so most users would not generally call this routine
737:    themselves.

739:    Level: developer

741: .seealso: NEPSetFunction(), NEPGetFunction()
742: @*/
743: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
744: {
745:   PetscInt       i;
746:   PetscScalar    alpha;

748:   PetscFunctionBegin;
750:   NEPCheckProblem(nep,1);
751:   switch (nep->fui) {
752:   case NEP_USER_INTERFACE_CALLBACK:
753:     PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
754:     PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
755:     PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
756:     PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
757:     break;
758:   case NEP_USER_INTERFACE_SPLIT:
759:     PetscCall(MatZeroEntries(A));
760:     if (A != B) PetscCall(MatZeroEntries(B));
761:     for (i=0;i<nep->nt;i++) {
762:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
763:       PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
764:       if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
765:     }
766:     break;
767:   }
768:   PetscFunctionReturn(PETSC_SUCCESS);
769: }

771: /*@
772:    NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
773:    set with NEPSetJacobian().

775:    Collective

777:    Input Parameters:
778: +  nep    - the NEP context
779: -  lambda - the scalar argument

781:    Output Parameters:
782: .  A   - Jacobian matrix

784:    Notes:
785:    Most users should not need to explicitly call this routine, as it
786:    is used internally within the nonlinear eigensolvers.

788:    Level: developer

790: .seealso: NEPSetJacobian(), NEPGetJacobian()
791: @*/
792: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
793: {
794:   PetscInt       i;
795:   PetscScalar    alpha;

797:   PetscFunctionBegin;
799:   NEPCheckProblem(nep,1);
800:   switch (nep->fui) {
801:   case NEP_USER_INTERFACE_CALLBACK:
802:     PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
803:     PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
804:     PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
805:     PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
806:     break;
807:   case NEP_USER_INTERFACE_SPLIT:
808:     PetscCall(MatZeroEntries(A));
809:     for (i=0;i<nep->nt;i++) {
810:       PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
811:       PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
812:     }
813:     break;
814:   }
815:   PetscFunctionReturn(PETSC_SUCCESS);
816: }