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.
78: Level: beginner
80: .seealso: [](ch:nep), `NEPCreate()`, `NEPSetUp()`, `NEPDestroy()`, `NEPSetTolerances()`, `NEPGetConverged()`, `NEPGetConvergedReason()`
81: @*/
82: PetscErrorCode NEPSolve(NEP nep)
83: {
84: PetscInt i;
85: char str[16];
87: PetscFunctionBegin;
89: if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
90: PetscCall(PetscCitationsRegister(citation,&cited));
91: PetscCall(PetscLogEventBegin(NEP_Solve,nep,0,0,0));
93: /* call setup */
94: PetscCall(NEPSetUp(nep));
95: nep->nconv = 0;
96: nep->its = 0;
97: for (i=0;i<nep->ncv;i++) {
98: nep->eigr[i] = 0.0;
99: nep->eigi[i] = 0.0;
100: nep->errest[i] = 0.0;
101: nep->perm[i] = i;
102: }
103: PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view_pre"));
104: PetscCall(RGViewFromOptions(nep->rg,NULL,"-rg_view"));
106: /* call solver */
107: PetscUseTypeMethod(nep,solve);
108: PetscCheck(nep->reason,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
109: nep->state = NEP_STATE_SOLVED;
111: /* Only the first nconv columns contain useful information */
112: PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
113: if (nep->twosided) PetscCall(BVSetActiveColumns(nep->W,0,nep->nconv));
115: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
116: PetscCall(NEPComputeVectors(nep));
117: PetscCall(NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv));
118: nep->state = NEP_STATE_EIGENVECTORS;
119: }
121: /* sort eigenvalues according to nep->which parameter */
122: PetscCall(SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm));
123: PetscCall(PetscLogEventEnd(NEP_Solve,nep,0,0,0));
125: /* various viewers */
126: PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
127: PetscCall(NEPConvergedReasonViewFromOptions(nep));
128: PetscCall(NEPErrorViewFromOptions(nep));
129: PetscCall(NEPValuesViewFromOptions(nep));
130: PetscCall(NEPVectorsViewFromOptions(nep));
131: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
132: for (i=0;i<nep->nt;i++) {
133: PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_mat%" PetscInt_FMT,i));
134: PetscCall(MatViewFromOptions(nep->A[i],(PetscObject)nep,str));
135: PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_fn%" PetscInt_FMT,i));
136: PetscCall(FNViewFromOptions(nep->f[i],(PetscObject)nep,str));
137: }
138: }
140: /* Remove the initial subspace */
141: nep->nini = 0;
143: /* Reset resolvent information */
144: PetscCall(MatDestroy(&nep->resolvent));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*@
149: NEPProjectOperator - Computes the projection of the nonlinear operator.
151: Collective
153: Input Parameters:
154: + nep - the nonlinear eigensolver context
155: . j0 - initial index
156: - j1 - final index
158: Notes:
159: This is available for split operator only.
161: The nonlinear operator $T(\lambda)$ is projected onto $\operatorname{span}(V)$, where $V$ is
162: an orthonormal basis built internally by the solver. The projected
163: operator is equal to $\sum_i V^* A_i V f_i(\lambda)$, so this function
164: computes all matrices $E_i = V^* A_i V$, and stores them in the extra
165: matrices inside `DS`, see `DSMatType`. Only rows/columns in the range [`j0`,`j1`-1] are computed,
166: the previous ones are assumed to be available already.
168: Level: developer
170: .seealso: [](ch:nep), `NEPSetSplitOperator()`
171: @*/
172: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
173: {
174: PetscInt k;
175: Mat G;
177: PetscFunctionBegin;
181: NEPCheckProblem(nep,1);
182: NEPCheckSplit(nep,1);
183: PetscCall(BVSetActiveColumns(nep->V,j0,j1));
184: for (k=0;k<nep->nt;k++) {
185: PetscCall(DSGetMat(nep->ds,DSMatExtra[k],&G));
186: PetscCall(BVMatProject(nep->V,nep->A[k],nep->V,G));
187: PetscCall(DSRestoreMat(nep->ds,DSMatExtra[k],&G));
188: }
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: /*@
193: NEPApplyFunction - Applies the nonlinear function $T(\lambda)$ to a given vector.
195: Collective
197: Input Parameters:
198: + nep - the nonlinear eigensolver context
199: . lambda - scalar argument
200: . x - vector to be multiplied against
201: - v - workspace vector (used only in the case of split form)
203: Output Parameters:
204: + y - result vector
205: . A - (optional) Function matrix, for callback interface only
206: - B - (unused) preconditioning matrix
208: Note:
209: If the nonlinear operator is represented in split form, the result
210: $y = T(\lambda)x$ is computed without building $T(\lambda)$ explicitly. In
211: that case, parameters `A` and `B` are not used. Otherwise, the matrix
212: $T(\lambda)$ is built and the effect is the same as a call to
213: `NEPComputeFunction()` followed by a `MatMult()`.
215: Level: developer
217: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeFunction()`, `NEPApplyAdjoint()`
218: @*/
219: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
220: {
221: PetscInt i;
222: PetscScalar alpha;
224: PetscFunctionBegin;
233: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
234: PetscCall(VecSet(y,0.0));
235: for (i=0;i<nep->nt;i++) {
236: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
237: PetscCall(MatMult(nep->A[i],x,v));
238: PetscCall(VecAXPY(y,alpha,v));
239: }
240: } else {
241: if (!A) A = nep->function;
242: PetscCall(NEPComputeFunction(nep,lambda,A,A));
243: PetscCall(MatMult(A,x,y));
244: }
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: /*@
249: NEPApplyAdjoint - Applies the adjoint nonlinear function $T^*(\lambda)$ to a given vector.
251: Collective
253: Input Parameters:
254: + nep - the nonlinear eigensolver context
255: . lambda - scalar argument
256: . x - vector to be multiplied against
257: - v - workspace vector (used only in the case of split form)
259: Output Parameters:
260: + y - result vector
261: . A - (optional) Function matrix, for callback interface only
262: - B - (unused) preconditioning matrix
264: Note:
265: If the nonlinear operator is represented in split form, the result
266: $y = T^*(\lambda)x$ is computed without building $T(\lambda)$ explicitly. In
267: that case, parameters `A` and `B` are not used. Otherwise, the matrix
268: $T(\lambda)$ is built and the effect is the same as a call to
269: `NEPComputeFunction()` followed by a `MatMultHermitianTranspose()`.
271: Level: developer
273: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeFunction()`, `NEPApplyFunction()`
274: @*/
275: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
276: {
277: PetscInt i;
278: PetscScalar alpha;
279: Vec w;
281: PetscFunctionBegin;
290: PetscCall(VecDuplicate(x,&w));
291: PetscCall(VecCopy(x,w));
292: PetscCall(VecConjugate(w));
293: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
294: PetscCall(VecSet(y,0.0));
295: for (i=0;i<nep->nt;i++) {
296: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
297: PetscCall(MatMultHermitianTranspose(nep->A[i],w,v));
298: PetscCall(VecAXPY(y,alpha,v));
299: }
300: } else {
301: if (!A) A = nep->function;
302: PetscCall(NEPComputeFunction(nep,lambda,A,A));
303: PetscCall(MatMultHermitianTranspose(A,w,y));
304: }
305: PetscCall(VecDestroy(&w));
306: PetscCall(VecConjugate(y));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: /*@
311: NEPApplyJacobian - Applies the nonlinear Jacobian $T'(\lambda)$ to a given vector.
313: Collective
315: Input Parameters:
316: + nep - the nonlinear eigensolver context
317: . lambda - scalar argument
318: . x - vector to be multiplied against
319: - v - workspace vector (used only in the case of split form)
321: Output Parameters:
322: + y - result vector
323: - A - (optional) Jacobian matrix, for callback interface only
325: Note:
326: If the nonlinear operator is represented in split form, the result
327: $y = T'(\lambda)x$ is computed without building $T'(\lambda)$ explicitly. In
328: that case, parameter `A` is not used. Otherwise, the matrix
329: $T'(\lambda)$ is built and the effect is the same as a call to
330: `NEPComputeJacobian()` followed by a `MatMult()`.
332: Level: developer
334: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeJacobian()`
335: @*/
336: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
337: {
338: PetscInt i;
339: PetscScalar alpha;
341: PetscFunctionBegin;
349: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
350: PetscCall(VecSet(y,0.0));
351: for (i=0;i<nep->nt;i++) {
352: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
353: PetscCall(MatMult(nep->A[i],x,v));
354: PetscCall(VecAXPY(y,alpha,v));
355: }
356: } else {
357: if (!A) A = nep->jacobian;
358: PetscCall(NEPComputeJacobian(nep,lambda,A));
359: PetscCall(MatMult(A,x,y));
360: }
361: PetscFunctionReturn(PETSC_SUCCESS);
362: }
364: /*@
365: NEPGetIterationNumber - Gets the current iteration number. If the
366: call to `NEPSolve()` is complete, then it returns the number of iterations
367: carried out by the solution method.
369: Not Collective
371: Input Parameter:
372: . nep - the nonlinear eigensolver context
374: Output Parameter:
375: . its - number of iterations
377: Note:
378: During the $i$-th iteration this call returns $i-1$. If `NEPSolve()` is
379: complete, then parameter `its` contains either the iteration number at
380: which convergence was successfully reached, or failure was detected.
381: Call `NEPGetConvergedReason()` to determine if the solver converged or
382: failed and why.
384: Level: intermediate
386: .seealso: [](ch:nep), `NEPGetConvergedReason()`, `NEPSetTolerances()`
387: @*/
388: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
389: {
390: PetscFunctionBegin;
392: PetscAssertPointer(its,2);
393: *its = nep->its;
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: /*@
398: NEPGetConverged - Gets the number of converged eigenpairs.
400: Not Collective
402: Input Parameter:
403: . nep - the nonlinear eigensolver context
405: Output Parameter:
406: . nconv - number of converged eigenpairs
408: Note:
409: This function should be called after `NEPSolve()` has finished.
411: The value `nconv` may be different from the number of requested solutions
412: `nev`, but not larger than `ncv`, see `NEPSetDimensions()`.
414: Level: beginner
416: .seealso: [](ch:nep), `NEPSetDimensions()`, `NEPSolve()`, `NEPGetEigenpair()`
417: @*/
418: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
419: {
420: PetscFunctionBegin;
422: PetscAssertPointer(nconv,2);
423: NEPCheckSolved(nep,1);
424: *nconv = nep->nconv;
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: /*@
429: NEPGetConvergedReason - Gets the reason why the `NEPSolve()` iteration was
430: stopped.
432: Not Collective
434: Input Parameter:
435: . nep - the nonlinear eigensolver context
437: Output Parameter:
438: . reason - negative value indicates diverged, positive value converged, see
439: `NEPConvergedReason` for the possible values
441: Options Database Key:
442: . -nep_converged_reason - print reason for convergence/divergence, and number of iterations
444: Note:
445: If this routine is called before or doing the `NEPSolve()` the value of
446: `NEP_CONVERGED_ITERATING` is returned.
448: Level: intermediate
450: .seealso: [](ch:nep), `NEPSetTolerances()`, `NEPSolve()`, `NEPConvergedReason`
451: @*/
452: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
453: {
454: PetscFunctionBegin;
456: PetscAssertPointer(reason,2);
457: NEPCheckSolved(nep,1);
458: *reason = nep->reason;
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }
462: /*@
463: NEPGetEigenpair - Gets the `i`-th solution of the eigenproblem as computed by
464: `NEPSolve()`. The solution consists in both the eigenvalue and the eigenvector.
466: Collective
468: Input Parameters:
469: + nep - the nonlinear eigensolver context
470: - i - index of the solution
472: Output Parameters:
473: + eigr - real part of eigenvalue
474: . eigi - imaginary part of eigenvalue
475: . Vr - real part of eigenvector
476: - Vi - imaginary part of eigenvector
478: Notes:
479: It is allowed to pass `NULL` for `Vr` and `Vi`, if the eigenvector is not
480: required. Otherwise, the caller must provide valid `Vec` objects, i.e.,
481: they must be created by the calling program with e.g. `MatCreateVecs()`.
483: If the eigenvalue is real, then `eigi` and `Vi` are set to zero. If PETSc is
484: configured with complex scalars the eigenvalue is stored
485: directly in `eigr` (`eigi` is set to zero) and the eigenvector in `Vr` (`Vi`
486: is set to zero). In any case, the user can pass `NULL` in any argument that
487: is not required.
489: The index `i` should be a value between 0 and `nconv`-1 (see `NEPGetConverged()`).
490: Eigenpairs are indexed according to the ordering criterion established
491: with `NEPSetWhichEigenpairs()`.
493: The eigenvector is normalized to have unit norm.
495: Level: beginner
497: .seealso: [](ch:nep), `NEPSolve()`, `NEPGetConverged()`, `NEPSetWhichEigenpairs()`, `NEPGetLeftEigenvector()`
498: @*/
499: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
500: {
501: PetscInt k;
503: PetscFunctionBegin;
508: NEPCheckSolved(nep,1);
509: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
510: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
512: PetscCall(NEPComputeVectors(nep));
513: k = nep->perm[i];
515: /* eigenvalue */
516: #if defined(PETSC_USE_COMPLEX)
517: if (eigr) *eigr = nep->eigr[k];
518: if (eigi) *eigi = 0;
519: #else
520: if (eigr) *eigr = nep->eigr[k];
521: if (eigi) *eigi = nep->eigi[k];
522: #endif
524: /* eigenvector */
525: PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
526: PetscFunctionReturn(PETSC_SUCCESS);
527: }
529: /*@
530: NEPGetLeftEigenvector - Gets the `i`-th left eigenvector as computed by `NEPSolve()`.
532: Collective
534: Input Parameters:
535: + nep - the nonlinear eigensolver context
536: - i - index of the solution
538: Output Parameters:
539: + Wr - real part of left eigenvector
540: - Wi - imaginary part of left eigenvector
542: Notes:
543: The caller must provide valid `Vec` objects, i.e., they must be created
544: by the calling program with e.g. `MatCreateVecs()`.
546: If the corresponding eigenvalue is real, then `Wi` is set to zero. If PETSc is
547: configured with complex scalars the eigenvector is stored directly in `Wr`
548: (`Wi` is set to zero). In any case, the user can pass `NULL` in `Wr` or `Wi` if
549: one of them is not required.
551: The index `i` should be a value between 0 and `nconv`-1 (see `NEPGetConverged()`).
552: Eigensolutions are indexed according to the ordering criterion established
553: with `NEPSetWhichEigenpairs()`.
555: Left eigenvectors are available only if the twosided flag was set, see
556: `NEPSetTwoSided()`.
558: Level: intermediate
560: .seealso: [](ch:nep), `NEPGetEigenpair()`, `NEPGetConverged()`, `NEPSetWhichEigenpairs()`, `NEPSetTwoSided()`
561: @*/
562: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
563: {
564: PetscInt k;
566: PetscFunctionBegin;
571: NEPCheckSolved(nep,1);
572: PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
573: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
574: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
575: PetscCall(NEPComputeVectors(nep));
576: k = nep->perm[i];
577: PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: /*@
582: NEPGetErrorEstimate - Returns the error estimate associated to the `i`-th
583: computed eigenpair.
585: Not Collective
587: Input Parameters:
588: + nep - the nonlinear eigensolver context
589: - i - index of eigenpair
591: Output Parameter:
592: . errest - the error estimate
594: Note:
595: This is the error estimate used internally by the eigensolver. The actual
596: error bound can be computed with `NEPComputeError()`.
598: Level: advanced
600: .seealso: [](ch:nep), `NEPComputeError()`
601: @*/
602: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
603: {
604: PetscFunctionBegin;
606: PetscAssertPointer(errest,3);
607: NEPCheckSolved(nep,1);
608: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
609: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
610: *errest = nep->errest[nep->perm[i]];
611: PetscFunctionReturn(PETSC_SUCCESS);
612: }
614: /*
615: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
616: associated with an eigenpair.
618: Input Parameters:
619: adj - whether the adjoint T^* must be used instead of T
620: lambda - eigenvalue
621: x - eigenvector
622: w - array of work vectors (two vectors in split form, one vector otherwise)
623: */
624: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
625: {
626: Vec y,z=NULL;
628: PetscFunctionBegin;
629: y = w[0];
630: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
631: if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
632: else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
633: PetscCall(VecNorm(y,NORM_2,norm));
634: PetscFunctionReturn(PETSC_SUCCESS);
635: }
637: /*@
638: NEPComputeError - Computes the error (based on the residual norm) associated
639: with the `i`-th computed eigenpair.
641: Collective
643: Input Parameters:
644: + nep - the nonlinear eigensolver context
645: . i - the solution index
646: - type - the type of error to compute, see `NEPErrorType`
648: Output Parameter:
649: . error - the error
651: Notes:
652: The error can be computed in various ways, all of them based on the residual
653: norm computed as $\|T(\lambda)x\|_2$ where $(\lambda,x)$ is the approximate eigenpair.
655: If the computation of left eigenvectors was enabled with `NEPSetTwoSided()`,
656: then the error will be computed using the maximum of the value above and
657: the left residual norm $\|y^*T(\lambda)\|_2$, where $y$ is the approximate left
658: eigenvector.
660: Level: beginner
662: .seealso: [](ch:nep), `NEPErrorType`, `NEPSolve()`, `NEPGetErrorEstimate()`, `NEPSetTwoSided()`
663: @*/
664: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
665: {
666: Vec xr,xi=NULL;
667: PetscInt j,nwork,issplit=0;
668: PetscScalar kr,ki,s;
669: PetscReal er,z=0.0,errorl,nrm;
670: PetscBool flg;
672: PetscFunctionBegin;
676: PetscAssertPointer(error,4);
677: NEPCheckSolved(nep,1);
679: /* allocate work vectors */
680: #if defined(PETSC_USE_COMPLEX)
681: nwork = 2;
682: #else
683: nwork = 3;
684: #endif
685: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
686: issplit = 1;
687: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
688: }
689: PetscCall(NEPSetWorkVecs(nep,nwork));
690: xr = nep->work[issplit+1];
691: #if !defined(PETSC_USE_COMPLEX)
692: xi = nep->work[issplit+2];
693: #endif
695: /* compute residual norms */
696: PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
697: #if !defined(PETSC_USE_COMPLEX)
698: PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
699: #endif
700: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
701: PetscCall(VecNorm(xr,NORM_2,&er));
703: /* if two-sided, compute left residual norm and take the maximum */
704: if (nep->twosided) {
705: PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
706: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
707: *error = PetscMax(*error,errorl);
708: }
710: /* compute error */
711: switch (type) {
712: case NEP_ERROR_ABSOLUTE:
713: break;
714: case NEP_ERROR_RELATIVE:
715: *error /= PetscAbsScalar(kr)*er;
716: break;
717: case NEP_ERROR_BACKWARD:
718: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
719: PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
720: PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
721: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
722: PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
723: *error /= nrm*er;
724: break;
725: }
726: /* initialization of matrix norms */
727: if (!nep->nrma[0]) {
728: for (j=0;j<nep->nt;j++) {
729: PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
730: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
731: PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
732: }
733: }
734: for (j=0;j<nep->nt;j++) {
735: PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
736: z = z + nep->nrma[j]*PetscAbsScalar(s);
737: }
738: *error /= z*er;
739: break;
740: default:
741: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
742: }
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: /*@
747: NEPComputeFunction - Computes the function matrix $T(\lambda)$ that has been
748: set with `NEPSetFunction()`.
750: Collective
752: Input Parameters:
753: + nep - the nonlinear eigensolver context
754: - lambda - the scalar argument
756: Output Parameters:
757: + A - Function matrix
758: - B - optional preconditioning matrix
760: Note:
761: `NEPComputeFunction()` is typically used within nonlinear eigensolvers
762: implementations, so most users would not generally call this routine
763: themselves.
765: Level: developer
767: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetFunction()`, `NEPComputeJacobian()`
768: @*/
769: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
770: {
771: PetscInt i;
772: PetscScalar alpha;
774: PetscFunctionBegin;
776: NEPCheckProblem(nep,1);
777: switch (nep->fui) {
778: case NEP_USER_INTERFACE_CALLBACK:
779: PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
780: PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
781: PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
782: PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
783: break;
784: case NEP_USER_INTERFACE_SPLIT:
785: PetscCall(MatZeroEntries(A));
786: if (A != B) PetscCall(MatZeroEntries(B));
787: for (i=0;i<nep->nt;i++) {
788: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
789: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
790: if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
791: }
792: break;
793: }
794: PetscFunctionReturn(PETSC_SUCCESS);
795: }
797: /*@
798: NEPComputeJacobian - Computes the Jacobian matrix $T'(\lambda)$ that has been
799: set with `NEPSetJacobian()`.
801: Collective
803: Input Parameters:
804: + nep - the nonlinear eigensolver context
805: - lambda - the scalar argument
807: Output Parameter:
808: . A - Jacobian matrix
810: Notes:
811: Most users should not need to explicitly call this routine, as it
812: is used internally within the nonlinear eigensolvers.
814: Level: developer
816: .seealso: [](ch:nep), `NEPSetJacobian()`, `NEPGetJacobian()`, `NEPComputeFunction()`
817: @*/
818: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
819: {
820: PetscInt i;
821: PetscScalar alpha;
823: PetscFunctionBegin;
825: NEPCheckProblem(nep,1);
826: switch (nep->fui) {
827: case NEP_USER_INTERFACE_CALLBACK:
828: PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
829: PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
830: PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
831: PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
832: break;
833: case NEP_USER_INTERFACE_SPLIT:
834: PetscCall(MatZeroEntries(A));
835: for (i=0;i<nep->nt;i++) {
836: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
837: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
838: }
839: break;
840: }
841: PetscFunctionReturn(PETSC_SUCCESS);
842: }