Actual source code: nepsolve.c

slepc-main 2024-11-09
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_matk - view the split form matrix Ak (replace k by an integer from 0 to nt-1)
 57: .  -nep_view_fnk - view the split form function fk (replace k by an integer from 0 to nt-1)
 58: .  -nep_view_vectors - view the computed eigenvectors
 59: .  -nep_view_values - view the computed eigenvalues
 60: .  -nep_converged_reason - print reason for convergence, and number of iterations
 61: .  -nep_error_absolute - print absolute errors of each eigenpair
 62: .  -nep_error_relative - print relative errors of each eigenpair
 63: -  -nep_error_backward - print backward errors of each eigenpair

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

 72:    Level: beginner

 74: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
 75: @*/
 76: PetscErrorCode NEPSolve(NEP nep)
 77: {
 78:   PetscInt       i;
 79:   char           str[16];

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

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

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

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

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

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

119:   /* various viewers */
120:   PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
121:   PetscCall(NEPConvergedReasonViewFromOptions(nep));
122:   PetscCall(NEPErrorViewFromOptions(nep));
123:   PetscCall(NEPValuesViewFromOptions(nep));
124:   PetscCall(NEPVectorsViewFromOptions(nep));
125:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
126:     for (i=0;i<nep->nt;i++) {
127:       PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_mat%" PetscInt_FMT,i));
128:       PetscCall(MatViewFromOptions(nep->A[i],(PetscObject)nep,str));
129:       PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_fn%" PetscInt_FMT,i));
130:       PetscCall(FNViewFromOptions(nep->f[i],(PetscObject)nep,str));
131:     }
132:   }

134:   /* Remove the initial subspace */
135:   nep->nini = 0;

137:   /* Reset resolvent information */
138:   PetscCall(MatDestroy(&nep->resolvent));
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: /*@
143:    NEPProjectOperator - Computes the projection of the nonlinear operator.

145:    Collective

147:    Input Parameters:
148: +  nep - the nonlinear eigensolver context
149: .  j0  - initial index
150: -  j1  - final index

152:    Notes:
153:    This is available for split operator only.

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

162:    Level: developer

164: .seealso: NEPSetSplitOperator()
165: @*/
166: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
167: {
168:   PetscInt       k;
169:   Mat            G;

171:   PetscFunctionBegin;
175:   NEPCheckProblem(nep,1);
176:   NEPCheckSplit(nep,1);
177:   PetscCall(BVSetActiveColumns(nep->V,j0,j1));
178:   for (k=0;k<nep->nt;k++) {
179:     PetscCall(DSGetMat(nep->ds,DSMatExtra[k],&G));
180:     PetscCall(BVMatProject(nep->V,nep->A[k],nep->V,G));
181:     PetscCall(DSRestoreMat(nep->ds,DSMatExtra[k],&G));
182:   }
183:   PetscFunctionReturn(PETSC_SUCCESS);
184: }

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

189:    Collective

191:    Input Parameters:
192: +  nep    - the nonlinear eigensolver context
193: .  lambda - scalar argument
194: .  x      - vector to be multiplied against
195: -  v      - workspace vector (used only in the case of split form)

197:    Output Parameters:
198: +  y   - result vector
199: .  A   - (optional) Function matrix, for callback interface only
200: -  B   - (unused) preconditioning matrix

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

209:    Level: developer

211: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
212: @*/
213: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
214: {
215:   PetscInt       i;
216:   PetscScalar    alpha;

218:   PetscFunctionBegin;

227:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
228:     PetscCall(VecSet(y,0.0));
229:     for (i=0;i<nep->nt;i++) {
230:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
231:       PetscCall(MatMult(nep->A[i],x,v));
232:       PetscCall(VecAXPY(y,alpha,v));
233:     }
234:   } else {
235:     if (!A) A = nep->function;
236:     PetscCall(NEPComputeFunction(nep,lambda,A,A));
237:     PetscCall(MatMult(A,x,y));
238:   }
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

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

245:    Collective

247:    Input Parameters:
248: +  nep    - the nonlinear eigensolver context
249: .  lambda - scalar argument
250: .  x      - vector to be multiplied against
251: -  v      - workspace vector (used only in the case of split form)

253:    Output Parameters:
254: +  y   - result vector
255: .  A   - (optional) Function matrix, for callback interface only
256: -  B   - (unused) preconditioning matrix

258:    Level: developer

260: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
261: @*/
262: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
263: {
264:   PetscInt       i;
265:   PetscScalar    alpha;
266:   Vec            w;

268:   PetscFunctionBegin;

277:   PetscCall(VecDuplicate(x,&w));
278:   PetscCall(VecCopy(x,w));
279:   PetscCall(VecConjugate(w));
280:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
281:     PetscCall(VecSet(y,0.0));
282:     for (i=0;i<nep->nt;i++) {
283:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
284:       PetscCall(MatMultTranspose(nep->A[i],w,v));
285:       PetscCall(VecAXPY(y,alpha,v));
286:     }
287:   } else {
288:     if (!A) A = nep->function;
289:     PetscCall(NEPComputeFunction(nep,lambda,A,A));
290:     PetscCall(MatMultTranspose(A,w,y));
291:   }
292:   PetscCall(VecDestroy(&w));
293:   PetscCall(VecConjugate(y));
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

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

300:    Collective

302:    Input Parameters:
303: +  nep    - the nonlinear eigensolver context
304: .  lambda - scalar argument
305: .  x      - vector to be multiplied against
306: -  v      - workspace vector (used only in the case of split form)

308:    Output Parameters:
309: +  y   - result vector
310: -  A   - (optional) Jacobian matrix, for callback interface only

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

319:    Level: developer

321: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
322: @*/
323: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
324: {
325:   PetscInt       i;
326:   PetscScalar    alpha;

328:   PetscFunctionBegin;

336:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
337:     PetscCall(VecSet(y,0.0));
338:     for (i=0;i<nep->nt;i++) {
339:       PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
340:       PetscCall(MatMult(nep->A[i],x,v));
341:       PetscCall(VecAXPY(y,alpha,v));
342:     }
343:   } else {
344:     if (!A) A = nep->jacobian;
345:     PetscCall(NEPComputeJacobian(nep,lambda,A));
346:     PetscCall(MatMult(A,x,y));
347:   }
348:   PetscFunctionReturn(PETSC_SUCCESS);
349: }

351: /*@
352:    NEPGetIterationNumber - Gets the current iteration number. If the
353:    call to NEPSolve() is complete, then it returns the number of iterations
354:    carried out by the solution method.

356:    Not Collective

358:    Input Parameter:
359: .  nep - the nonlinear eigensolver context

361:    Output Parameter:
362: .  its - number of iterations

364:    Note:
365:    During the i-th iteration this call returns i-1. If NEPSolve() is
366:    complete, then parameter "its" contains either the iteration number at
367:    which convergence was successfully reached, or failure was detected.
368:    Call NEPGetConvergedReason() to determine if the solver converged or
369:    failed and why.

371:    Level: intermediate

373: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
374: @*/
375: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
376: {
377:   PetscFunctionBegin;
379:   PetscAssertPointer(its,2);
380:   *its = nep->its;
381:   PetscFunctionReturn(PETSC_SUCCESS);
382: }

384: /*@
385:    NEPGetConverged - Gets the number of converged eigenpairs.

387:    Not Collective

389:    Input Parameter:
390: .  nep - the nonlinear eigensolver context

392:    Output Parameter:
393: .  nconv - number of converged eigenpairs

395:    Note:
396:    This function should be called after NEPSolve() has finished.

398:    Level: beginner

400: .seealso: NEPSetDimensions(), NEPSolve(), NEPGetEigenpair()
401: @*/
402: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
403: {
404:   PetscFunctionBegin;
406:   PetscAssertPointer(nconv,2);
407:   NEPCheckSolved(nep,1);
408:   *nconv = nep->nconv;
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: /*@
413:    NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
414:    stopped.

416:    Not Collective

418:    Input Parameter:
419: .  nep - the nonlinear eigensolver context

421:    Output Parameter:
422: .  reason - negative value indicates diverged, positive value converged

424:    Options Database Key:
425: .  -nep_converged_reason - print the reason to a viewer

427:    Notes:
428:    Possible values for reason are
429: +  NEP_CONVERGED_TOL - converged up to tolerance
430: .  NEP_CONVERGED_USER - converged due to a user-defined condition
431: .  NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
432: .  NEP_DIVERGED_BREAKDOWN - generic breakdown in method
433: .  NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
434: -  NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
435:    unrestarted solver

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

439:    Level: intermediate

441: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason
442: @*/
443: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
444: {
445:   PetscFunctionBegin;
447:   PetscAssertPointer(reason,2);
448:   NEPCheckSolved(nep,1);
449:   *reason = nep->reason;
450:   PetscFunctionReturn(PETSC_SUCCESS);
451: }

453: /*@
454:    NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
455:    NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.

457:    Collective

459:    Input Parameters:
460: +  nep - nonlinear eigensolver context
461: -  i   - index of the solution

463:    Output Parameters:
464: +  eigr - real part of eigenvalue
465: .  eigi - imaginary part of eigenvalue
466: .  Vr   - real part of eigenvector
467: -  Vi   - imaginary part of eigenvector

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

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

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

484:    Level: beginner

486: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
487: @*/
488: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
489: {
490:   PetscInt       k;

492:   PetscFunctionBegin;
497:   NEPCheckSolved(nep,1);
498:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
499:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");

501:   PetscCall(NEPComputeVectors(nep));
502:   k = nep->perm[i];

504:   /* eigenvalue */
505: #if defined(PETSC_USE_COMPLEX)
506:   if (eigr) *eigr = nep->eigr[k];
507:   if (eigi) *eigi = 0;
508: #else
509:   if (eigr) *eigr = nep->eigr[k];
510:   if (eigi) *eigi = nep->eigi[k];
511: #endif

513:   /* eigenvector */
514:   PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
515:   PetscFunctionReturn(PETSC_SUCCESS);
516: }

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

521:    Collective

523:    Input Parameters:
524: +  nep - eigensolver context
525: -  i   - index of the solution

527:    Output Parameters:
528: +  Wr   - real part of left eigenvector
529: -  Wi   - imaginary part of left eigenvector

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

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

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

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

547:    Level: intermediate

549: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
550: @*/
551: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
552: {
553:   PetscInt       k;

555:   PetscFunctionBegin;
560:   NEPCheckSolved(nep,1);
561:   PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
562:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
563:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
564:   PetscCall(NEPComputeVectors(nep));
565:   k = nep->perm[i];
566:   PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
567:   PetscFunctionReturn(PETSC_SUCCESS);
568: }

570: /*@
571:    NEPGetErrorEstimate - Returns the error estimate associated to the i-th
572:    computed eigenpair.

574:    Not Collective

576:    Input Parameters:
577: +  nep - nonlinear eigensolver context
578: -  i   - index of eigenpair

580:    Output Parameter:
581: .  errest - the error estimate

583:    Notes:
584:    This is the error estimate used internally by the eigensolver. The actual
585:    error bound can be computed with NEPComputeError().

587:    Level: advanced

589: .seealso: NEPComputeError()
590: @*/
591: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
592: {
593:   PetscFunctionBegin;
595:   PetscAssertPointer(errest,3);
596:   NEPCheckSolved(nep,1);
597:   PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
598:   PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
599:   *errest = nep->errest[nep->perm[i]];
600:   PetscFunctionReturn(PETSC_SUCCESS);
601: }

603: /*
604:    NEPComputeResidualNorm_Private - Computes the norm of the residual vector
605:    associated with an eigenpair.

607:    Input Parameters:
608:      adj    - whether the adjoint T^* must be used instead of T
609:      lambda - eigenvalue
610:      x      - eigenvector
611:      w      - array of work vectors (two vectors in split form, one vector otherwise)
612: */
613: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
614: {
615:   Vec            y,z=NULL;

617:   PetscFunctionBegin;
618:   y = w[0];
619:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
620:   if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
621:   else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
622:   PetscCall(VecNorm(y,NORM_2,norm));
623:   PetscFunctionReturn(PETSC_SUCCESS);
624: }

626: /*@
627:    NEPComputeError - Computes the error (based on the residual norm) associated
628:    with the i-th computed eigenpair.

630:    Collective

632:    Input Parameters:
633: +  nep  - the nonlinear eigensolver context
634: .  i    - the solution index
635: -  type - the type of error to compute

637:    Output Parameter:
638: .  error - the error

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

645:    Level: beginner

647: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
648: @*/
649: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
650: {
651:   Vec            xr,xi=NULL;
652:   PetscInt       j,nwork,issplit=0;
653:   PetscScalar    kr,ki,s;
654:   PetscReal      er,z=0.0,errorl,nrm;
655:   PetscBool      flg;

657:   PetscFunctionBegin;
661:   PetscAssertPointer(error,4);
662:   NEPCheckSolved(nep,1);

664:   /* allocate work vectors */
665: #if defined(PETSC_USE_COMPLEX)
666:   nwork = 2;
667: #else
668:   nwork = 3;
669: #endif
670:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
671:     issplit = 1;
672:     nwork++;  /* need an extra work vector for NEPComputeResidualNorm_Private */
673:   }
674:   PetscCall(NEPSetWorkVecs(nep,nwork));
675:   xr = nep->work[issplit+1];
676: #if !defined(PETSC_USE_COMPLEX)
677:   xi = nep->work[issplit+2];
678: #endif

680:   /* compute residual norms */
681:   PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
682: #if !defined(PETSC_USE_COMPLEX)
683:   PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
684: #endif
685:   PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
686:   PetscCall(VecNorm(xr,NORM_2,&er));

688:   /* if two-sided, compute left residual norm and take the maximum */
689:   if (nep->twosided) {
690:     PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
691:     PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
692:     *error = PetscMax(*error,errorl);
693:   }

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

731: /*@
732:    NEPComputeFunction - Computes the function matrix T(lambda) that has been
733:    set with NEPSetFunction().

735:    Collective

737:    Input Parameters:
738: +  nep    - the NEP context
739: -  lambda - the scalar argument

741:    Output Parameters:
742: +  A   - Function matrix
743: -  B   - optional preconditioning matrix

745:    Notes:
746:    NEPComputeFunction() is typically used within nonlinear eigensolvers
747:    implementations, so most users would not generally call this routine
748:    themselves.

750:    Level: developer

752: .seealso: NEPSetFunction(), NEPGetFunction()
753: @*/
754: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
755: {
756:   PetscInt       i;
757:   PetscScalar    alpha;

759:   PetscFunctionBegin;
761:   NEPCheckProblem(nep,1);
762:   switch (nep->fui) {
763:   case NEP_USER_INTERFACE_CALLBACK:
764:     PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
765:     PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
766:     PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
767:     PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
768:     break;
769:   case NEP_USER_INTERFACE_SPLIT:
770:     PetscCall(MatZeroEntries(A));
771:     if (A != B) PetscCall(MatZeroEntries(B));
772:     for (i=0;i<nep->nt;i++) {
773:       PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
774:       PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
775:       if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
776:     }
777:     break;
778:   }
779:   PetscFunctionReturn(PETSC_SUCCESS);
780: }

782: /*@
783:    NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
784:    set with NEPSetJacobian().

786:    Collective

788:    Input Parameters:
789: +  nep    - the NEP context
790: -  lambda - the scalar argument

792:    Output Parameters:
793: .  A   - Jacobian matrix

795:    Notes:
796:    Most users should not need to explicitly call this routine, as it
797:    is used internally within the nonlinear eigensolvers.

799:    Level: developer

801: .seealso: NEPSetJacobian(), NEPGetJacobian()
802: @*/
803: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
804: {
805:   PetscInt       i;
806:   PetscScalar    alpha;

808:   PetscFunctionBegin;
810:   NEPCheckProblem(nep,1);
811:   switch (nep->fui) {
812:   case NEP_USER_INTERFACE_CALLBACK:
813:     PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
814:     PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
815:     PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
816:     PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
817:     break;
818:   case NEP_USER_INTERFACE_SPLIT:
819:     PetscCall(MatZeroEntries(A));
820:     for (i=0;i<nep->nt;i++) {
821:       PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
822:       PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
823:     }
824:     break;
825:   }
826:   PetscFunctionReturn(PETSC_SUCCESS);
827: }