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: *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: }