Actual source code: nepsolve.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: NEP routines related to the solution process
13: References:
15: [1] C. Campos and J.E. Roman, "NEP: a module for the parallel solution
16: of nonlinear eigenvalue problems in SLEPc", ACM Trans. Math. Soft.
17: 47(3), 23:1--23:29, 2021.
18: */
20: #include <slepc/private/nepimpl.h>
21: #include <slepc/private/bvimpl.h>
22: #include <petscdraw.h>
24: static PetscBool cited = PETSC_FALSE;
25: static const char citation[] =
26: "@Article{slepc-nep,\n"
27: " author = \"C. Campos and J. E. Roman\",\n"
28: " title = \"{NEP}: a module for the parallel solution of nonlinear eigenvalue problems in {SLEPc}\",\n"
29: " journal = \"{ACM} Trans. Math. Software\",\n"
30: " volume = \"47\",\n"
31: " number = \"3\",\n"
32: " pages = \"23:1--23:29\",\n"
33: " year = \"2021\",\n"
34: " doi = \"10.1145/3447544\"\n"
35: "}\n";
37: PetscErrorCode NEPComputeVectors(NEP nep)
38: {
39: PetscFunctionBegin;
40: NEPCheckSolved(nep,1);
41: if (nep->state==NEP_STATE_SOLVED) PetscTryTypeMethod(nep,computevectors);
42: nep->state = NEP_STATE_EIGENVECTORS;
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*@
47: NEPSolve - Solves the nonlinear eigenproblem.
49: Collective
51: Input Parameter:
52: . nep - the nonlinear eigensolver context
54: Options Database Keys:
55: + -nep_view - print information about the solver once the solve is complete
56: . -nep_view_pre - print information about the solver before the solve starts
57: . -nep_view_matk - view the split form matrix $A_k$ (replace `k` by an integer from 0 to `nt`-1)
58: . -nep_view_fnk - view the split form function $f_k$ (replace `k` by an integer from 0 to `nt`-1)
59: . -nep_view_vectors - view the computed eigenvectors
60: . -nep_view_values - view the computed eigenvalues
61: . -nep_converged_reason - print reason for convergence/divergence, and number of iterations
62: . -nep_error_absolute - print absolute errors of each eigenpair
63: . -nep_error_relative - print relative errors of each eigenpair
64: - -nep_error_backward - print backward errors of each eigenpair
66: Notes:
67: `NEPSolve()` will return without generating an error regardless of whether
68: all requested solutions were computed or not. Call `NEPGetConverged()` to get the
69: actual number of computed solutions, and `NEPGetConvergedReason()` to determine if
70: the solver converged or failed and why.
72: All the command-line options listed above admit an optional argument specifying
73: the viewer type and options. For instance, use `-nep_view_vectors binary:myvecs.bin`
74: to save the eigenvectors to a binary file, `-nep_view_values draw` to draw the computed
75: eigenvalues graphically, or `-nep_error_relative :myerr.m:ascii_matlab` to save
76: the errors in a file that can be executed in Matlab.
77: See `PetscObjectViewFromOptions()` for more details.
79: Level: beginner
81: .seealso: [](ch:nep), `NEPCreate()`, `NEPSetUp()`, `NEPDestroy()`, `NEPSetTolerances()`, `NEPGetConverged()`, `NEPGetConvergedReason()`
82: @*/
83: PetscErrorCode NEPSolve(NEP nep)
84: {
85: PetscInt i;
86: char str[16];
88: PetscFunctionBegin;
90: if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
91: PetscCall(PetscCitationsRegister(citation,&cited));
92: PetscCall(PetscLogEventBegin(NEP_Solve,nep,0,0,0));
94: /* call setup */
95: PetscCall(NEPSetUp(nep));
96: nep->nconv = 0;
97: nep->its = 0;
98: for (i=0;i<nep->ncv;i++) {
99: nep->eigr[i] = 0.0;
100: nep->eigi[i] = 0.0;
101: nep->errest[i] = 0.0;
102: nep->perm[i] = i;
103: }
104: PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view_pre"));
105: PetscCall(RGViewFromOptions(nep->rg,NULL,"-rg_view"));
107: /* call solver */
108: PetscUseTypeMethod(nep,solve);
109: PetscCheck(nep->reason,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
110: nep->state = NEP_STATE_SOLVED;
112: /* Only the first nconv columns contain useful information */
113: PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
114: if (nep->twosided) PetscCall(BVSetActiveColumns(nep->W,0,nep->nconv));
116: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
117: PetscCall(NEPComputeVectors(nep));
118: PetscCall(NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv));
119: nep->state = NEP_STATE_EIGENVECTORS;
120: }
122: /* sort eigenvalues according to nep->which parameter */
123: PetscCall(SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm));
124: PetscCall(PetscLogEventEnd(NEP_Solve,nep,0,0,0));
126: /* various viewers */
127: PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
128: PetscCall(NEPConvergedReasonViewFromOptions(nep));
129: PetscCall(NEPErrorViewFromOptions(nep));
130: PetscCall(NEPValuesViewFromOptions(nep));
131: PetscCall(NEPVectorsViewFromOptions(nep));
132: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
133: for (i=0;i<nep->nt;i++) {
134: PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_mat%" PetscInt_FMT,i));
135: PetscCall(MatViewFromOptions(nep->A[i],(PetscObject)nep,str));
136: PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_fn%" PetscInt_FMT,i));
137: PetscCall(FNViewFromOptions(nep->f[i],(PetscObject)nep,str));
138: }
139: }
141: /* Remove the initial subspace */
142: nep->nini = 0;
144: /* Reset resolvent information */
145: PetscCall(MatDestroy(&nep->resolvent));
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: /*@
150: NEPProjectOperator - Computes the projection of the nonlinear operator.
152: Collective
154: Input Parameters:
155: + nep - the nonlinear eigensolver context
156: . j0 - initial index
157: - j1 - final index
159: Notes:
160: This is available for split operator only.
162: The nonlinear operator $T(\lambda)$ is projected onto $\operatorname{span}(V)$, where $V$ is
163: an orthonormal basis built internally by the solver. The projected
164: operator is equal to $\sum_i V^* A_i V f_i(\lambda)$, so this function
165: computes all matrices $E_i = V^* A_i V$, and stores them in the extra
166: matrices inside `DS`, see `DSMatType`. Only rows/columns in the range [`j0`,`j1`-1] are computed,
167: the previous ones are assumed to be available already.
169: Level: developer
171: .seealso: [](ch:nep), `NEPSetSplitOperator()`
172: @*/
173: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
174: {
175: PetscInt k;
176: Mat G;
178: PetscFunctionBegin;
182: NEPCheckProblem(nep,1);
183: NEPCheckSplit(nep,1);
184: PetscCall(BVSetActiveColumns(nep->V,j0,j1));
185: for (k=0;k<nep->nt;k++) {
186: PetscCall(DSGetMat(nep->ds,DSMatExtra[k],&G));
187: PetscCall(BVMatProject(nep->V,nep->A[k],nep->V,G));
188: PetscCall(DSRestoreMat(nep->ds,DSMatExtra[k],&G));
189: }
190: PetscFunctionReturn(PETSC_SUCCESS);
191: }
193: /*@
194: NEPApplyFunction - Applies the nonlinear function $T(\lambda)$ to a given vector.
196: Collective
198: Input Parameters:
199: + nep - the nonlinear eigensolver context
200: . lambda - scalar argument
201: . x - vector to be multiplied against
202: - v - workspace vector (used only in the case of split form)
204: Output Parameters:
205: + y - result vector
206: . A - (optional) Function matrix, for callback interface only
207: - B - (unused) preconditioning matrix
209: Note:
210: If the nonlinear operator is represented in split form, the result
211: $y = T(\lambda)x$ is computed without building $T(\lambda)$ explicitly. In
212: that case, parameters `A` and `B` are not used. Otherwise, the matrix
213: $T(\lambda)$ is built and the effect is the same as a call to
214: `NEPComputeFunction()` followed by a `MatMult()`.
216: Level: developer
218: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeFunction()`, `NEPApplyAdjoint()`
219: @*/
220: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
221: {
222: PetscInt i;
223: PetscScalar alpha;
225: PetscFunctionBegin;
234: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
235: PetscCall(VecSet(y,0.0));
236: for (i=0;i<nep->nt;i++) {
237: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
238: PetscCall(MatMult(nep->A[i],x,v));
239: PetscCall(VecAXPY(y,alpha,v));
240: }
241: } else {
242: if (!A) A = nep->function;
243: PetscCall(NEPComputeFunction(nep,lambda,A,A));
244: PetscCall(MatMult(A,x,y));
245: }
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: /*@
250: NEPApplyAdjoint - Applies the adjoint nonlinear function $T^*(\lambda)$ to a given vector.
252: Collective
254: Input Parameters:
255: + nep - the nonlinear eigensolver context
256: . lambda - scalar argument
257: . x - vector to be multiplied against
258: - v - workspace vector (used only in the case of split form)
260: Output Parameters:
261: + y - result vector
262: . A - (optional) Function matrix, for callback interface only
263: - B - (unused) preconditioning matrix
265: Note:
266: If the nonlinear operator is represented in split form, the result
267: $y = T^*(\lambda)x$ is computed without building $T(\lambda)$ explicitly. In
268: that case, parameters `A` and `B` are not used. Otherwise, the matrix
269: $T(\lambda)$ is built and the effect is the same as a call to
270: `NEPComputeFunction()` followed by a `MatMultHermitianTranspose()`.
272: Level: developer
274: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeFunction()`, `NEPApplyFunction()`
275: @*/
276: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
277: {
278: PetscInt i;
279: PetscScalar alpha;
280: Vec w;
282: PetscFunctionBegin;
291: PetscCall(VecDuplicate(x,&w));
292: PetscCall(VecCopy(x,w));
293: PetscCall(VecConjugate(w));
294: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
295: PetscCall(VecSet(y,0.0));
296: for (i=0;i<nep->nt;i++) {
297: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
298: PetscCall(MatMultHermitianTranspose(nep->A[i],w,v));
299: PetscCall(VecAXPY(y,alpha,v));
300: }
301: } else {
302: if (!A) A = nep->function;
303: PetscCall(NEPComputeFunction(nep,lambda,A,A));
304: PetscCall(MatMultHermitianTranspose(A,w,y));
305: }
306: PetscCall(VecDestroy(&w));
307: PetscCall(VecConjugate(y));
308: PetscFunctionReturn(PETSC_SUCCESS);
309: }
311: /*@
312: NEPApplyJacobian - Applies the nonlinear Jacobian $T'(\lambda)$ to a given vector.
314: Collective
316: Input Parameters:
317: + nep - the nonlinear eigensolver context
318: . lambda - scalar argument
319: . x - vector to be multiplied against
320: - v - workspace vector (used only in the case of split form)
322: Output Parameters:
323: + y - result vector
324: - A - (optional) Jacobian matrix, for callback interface only
326: Note:
327: If the nonlinear operator is represented in split form, the result
328: $y = T'(\lambda)x$ is computed without building $T'(\lambda)$ explicitly. In
329: that case, parameter `A` is not used. Otherwise, the matrix
330: $T'(\lambda)$ is built and the effect is the same as a call to
331: `NEPComputeJacobian()` followed by a `MatMult()`.
333: Level: developer
335: .seealso: [](ch:nep), `NEPSetSplitOperator()`, `NEPComputeJacobian()`
336: @*/
337: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
338: {
339: PetscInt i;
340: PetscScalar alpha;
342: PetscFunctionBegin;
350: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
351: PetscCall(VecSet(y,0.0));
352: for (i=0;i<nep->nt;i++) {
353: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
354: PetscCall(MatMult(nep->A[i],x,v));
355: PetscCall(VecAXPY(y,alpha,v));
356: }
357: } else {
358: if (!A) A = nep->jacobian;
359: PetscCall(NEPComputeJacobian(nep,lambda,A));
360: PetscCall(MatMult(A,x,y));
361: }
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@
366: NEPGetIterationNumber - Gets the current iteration number. If the
367: call to `NEPSolve()` is complete, then it returns the number of iterations
368: carried out by the solution method.
370: Not Collective
372: Input Parameter:
373: . nep - the nonlinear eigensolver context
375: Output Parameter:
376: . its - number of iterations
378: Note:
379: During the $i$-th iteration this call returns $i-1$. If `NEPSolve()` is
380: complete, then parameter `its` contains either the iteration number at
381: which convergence was successfully reached, or failure was detected.
382: Call `NEPGetConvergedReason()` to determine if the solver converged or
383: failed and why.
385: Level: intermediate
387: .seealso: [](ch:nep), `NEPGetConvergedReason()`, `NEPSetTolerances()`
388: @*/
389: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
390: {
391: PetscFunctionBegin;
393: PetscAssertPointer(its,2);
394: *its = nep->its;
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: /*@
399: NEPGetConverged - Gets the number of converged eigenpairs.
401: Not Collective
403: Input Parameter:
404: . nep - the nonlinear eigensolver context
406: Output Parameter:
407: . nconv - number of converged eigenpairs
409: Note:
410: This function should be called after `NEPSolve()` has finished.
412: The value `nconv` may be different from the number of requested solutions
413: `nev`, but not larger than `ncv`, see `NEPSetDimensions()`.
415: Level: beginner
417: .seealso: [](ch:nep), `NEPSetDimensions()`, `NEPSolve()`, `NEPGetEigenpair()`
418: @*/
419: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
420: {
421: PetscFunctionBegin;
423: PetscAssertPointer(nconv,2);
424: NEPCheckSolved(nep,1);
425: *nconv = nep->nconv;
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: /*@
430: NEPGetConvergedReason - Gets the reason why the `NEPSolve()` iteration was
431: stopped.
433: Not Collective
435: Input Parameter:
436: . nep - the nonlinear eigensolver context
438: Output Parameter:
439: . reason - negative value indicates diverged, positive value converged, see
440: `NEPConvergedReason` for the possible values
442: Options Database Key:
443: . -nep_converged_reason - print reason for convergence/divergence, and number of iterations
445: Note:
446: If this routine is called before or doing the `NEPSolve()` the value of
447: `NEP_CONVERGED_ITERATING` is returned.
449: Level: intermediate
451: .seealso: [](ch:nep), `NEPSetTolerances()`, `NEPSolve()`, `NEPConvergedReason`
452: @*/
453: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
454: {
455: PetscFunctionBegin;
457: PetscAssertPointer(reason,2);
458: NEPCheckSolved(nep,1);
459: *reason = nep->reason;
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: /*@
464: NEPGetEigenpair - Gets the `i`-th solution of the eigenproblem as computed by
465: `NEPSolve()`. The solution consists in both the eigenvalue and the eigenvector.
467: Collective
469: Input Parameters:
470: + nep - the nonlinear eigensolver context
471: - i - index of the solution
473: Output Parameters:
474: + eigr - real part of eigenvalue
475: . eigi - imaginary part of eigenvalue
476: . Vr - real part of eigenvector
477: - Vi - imaginary part of eigenvector
479: Notes:
480: It is allowed to pass `NULL` for `Vr` and `Vi`, if the eigenvector is not
481: required. Otherwise, the caller must provide valid `Vec` objects, i.e.,
482: they must be created by the calling program with e.g. `MatCreateVecs()`.
484: If the eigenvalue is real, then `eigi` and `Vi` are set to zero. If PETSc is
485: configured with complex scalars the eigenvalue is stored
486: directly in `eigr` (`eigi` is set to zero) and the eigenvector in `Vr` (`Vi`
487: is set to zero). In any case, the user can pass `NULL` in any argument that
488: is not required.
490: The index `i` should be a value between 0 and `nconv`-1 (see `NEPGetConverged()`).
491: Eigenpairs are indexed according to the ordering criterion established
492: with `NEPSetWhichEigenpairs()`.
494: The eigenvector is normalized to have unit norm.
496: Level: beginner
498: .seealso: [](ch:nep), `NEPSolve()`, `NEPGetConverged()`, `NEPSetWhichEigenpairs()`, `NEPGetLeftEigenvector()`
499: @*/
500: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
501: {
502: PetscInt k;
504: PetscFunctionBegin;
509: NEPCheckSolved(nep,1);
510: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
511: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
513: PetscCall(NEPComputeVectors(nep));
514: k = nep->perm[i];
516: /* eigenvalue */
517: #if defined(PETSC_USE_COMPLEX)
518: if (eigr) *eigr = nep->eigr[k];
519: if (eigi) *eigi = 0;
520: #else
521: if (eigr) *eigr = nep->eigr[k];
522: if (eigi) *eigi = nep->eigi[k];
523: #endif
525: /* eigenvector */
526: PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
527: PetscFunctionReturn(PETSC_SUCCESS);
528: }
530: /*@
531: NEPGetLeftEigenvector - Gets the `i`-th left eigenvector as computed by `NEPSolve()`.
533: Collective
535: Input Parameters:
536: + nep - the nonlinear eigensolver context
537: - i - index of the solution
539: Output Parameters:
540: + Wr - real part of left eigenvector
541: - Wi - imaginary part of left eigenvector
543: Notes:
544: The caller must provide valid `Vec` objects, i.e., they must be created
545: by the calling program with e.g. `MatCreateVecs()`.
547: If the corresponding eigenvalue is real, then `Wi` is set to zero. If PETSc is
548: configured with complex scalars the eigenvector is stored directly in `Wr`
549: (`Wi` is set to zero). In any case, the user can pass `NULL` in `Wr` or `Wi` if
550: one of them is not required.
552: The index `i` should be a value between 0 and `nconv`-1 (see `NEPGetConverged()`).
553: Eigensolutions are indexed according to the ordering criterion established
554: with `NEPSetWhichEigenpairs()`.
556: Left eigenvectors are available only if the twosided flag was set, see
557: `NEPSetTwoSided()`.
559: Level: intermediate
561: .seealso: [](ch:nep), `NEPGetEigenpair()`, `NEPGetConverged()`, `NEPSetWhichEigenpairs()`, `NEPSetTwoSided()`
562: @*/
563: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
564: {
565: PetscInt k;
567: PetscFunctionBegin;
572: NEPCheckSolved(nep,1);
573: PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
574: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
575: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
576: PetscCall(NEPComputeVectors(nep));
577: k = nep->perm[i];
578: PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /*@
583: NEPGetErrorEstimate - Returns the error estimate associated to the `i`-th
584: computed eigenpair.
586: Not Collective
588: Input Parameters:
589: + nep - the nonlinear eigensolver context
590: - i - index of eigenpair
592: Output Parameter:
593: . errest - the error estimate
595: Note:
596: This is the error estimate used internally by the eigensolver. The actual
597: error bound can be computed with `NEPComputeError()`.
599: Level: advanced
601: .seealso: [](ch:nep), `NEPComputeError()`
602: @*/
603: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
604: {
605: PetscFunctionBegin;
607: PetscAssertPointer(errest,3);
608: NEPCheckSolved(nep,1);
609: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
610: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
611: *errest = nep->errest[nep->perm[i]];
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*
616: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
617: associated with an eigenpair.
619: Input Parameters:
620: adj - whether the adjoint T^* must be used instead of T
621: lambda - eigenvalue
622: x - eigenvector
623: w - array of work vectors (two vectors in split form, one vector otherwise)
624: */
625: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
626: {
627: Vec y,z=NULL;
629: PetscFunctionBegin;
630: y = w[0];
631: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
632: if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
633: else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
634: PetscCall(VecNorm(y,NORM_2,norm));
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: /*@
639: NEPComputeError - Computes the error (based on the residual norm) associated
640: with the `i`-th computed eigenpair.
642: Collective
644: Input Parameters:
645: + nep - the nonlinear eigensolver context
646: . i - the solution index
647: - type - the type of error to compute, see `NEPErrorType`
649: Output Parameter:
650: . error - the error
652: Notes:
653: The error can be computed in various ways, all of them based on the residual
654: norm computed as $\|T(\lambda)x\|_2$ where $(\lambda,x)$ is the approximate eigenpair.
656: If the computation of left eigenvectors was enabled with `NEPSetTwoSided()`,
657: then the error will be computed using the maximum of the value above and
658: the left residual norm $\|y^*T(\lambda)\|_2$, where $y$ is the approximate left
659: eigenvector.
661: Level: beginner
663: .seealso: [](ch:nep), `NEPErrorType`, `NEPSolve()`, `NEPGetErrorEstimate()`, `NEPSetTwoSided()`
664: @*/
665: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
666: {
667: Vec xr,xi=NULL;
668: PetscInt j,nwork,issplit=0;
669: PetscScalar kr,ki,s;
670: PetscReal er,z=0.0,errorl,nrm;
671: PetscBool flg;
673: PetscFunctionBegin;
677: PetscAssertPointer(error,4);
678: NEPCheckSolved(nep,1);
680: /* allocate work vectors */
681: #if defined(PETSC_USE_COMPLEX)
682: nwork = 2;
683: #else
684: nwork = 3;
685: #endif
686: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
687: issplit = 1;
688: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
689: }
690: PetscCall(NEPSetWorkVecs(nep,nwork));
691: xr = nep->work[issplit+1];
692: #if !defined(PETSC_USE_COMPLEX)
693: xi = nep->work[issplit+2];
694: #endif
696: /* compute residual norms */
697: PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
698: #if !defined(PETSC_USE_COMPLEX)
699: PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
700: #endif
701: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
702: PetscCall(VecNorm(xr,NORM_2,&er));
704: /* if two-sided, compute left residual norm and take the maximum */
705: if (nep->twosided) {
706: PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
707: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
708: *error = PetscMax(*error,errorl);
709: }
711: /* compute error */
712: switch (type) {
713: case NEP_ERROR_ABSOLUTE:
714: break;
715: case NEP_ERROR_RELATIVE:
716: *error /= PetscAbsScalar(kr)*er;
717: break;
718: case NEP_ERROR_BACKWARD:
719: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
720: PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
721: PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
722: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
723: PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
724: *error /= nrm*er;
725: break;
726: }
727: /* initialization of matrix norms */
728: if (!nep->nrma[0]) {
729: for (j=0;j<nep->nt;j++) {
730: PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
731: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
732: PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
733: }
734: }
735: for (j=0;j<nep->nt;j++) {
736: PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
737: z = z + nep->nrma[j]*PetscAbsScalar(s);
738: }
739: *error /= z*er;
740: break;
741: default:
742: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
743: }
744: PetscFunctionReturn(PETSC_SUCCESS);
745: }
747: /*@
748: NEPComputeFunction - Computes the function matrix $T(\lambda)$ that has been
749: set with `NEPSetFunction()`.
751: Collective
753: Input Parameters:
754: + nep - the nonlinear eigensolver context
755: - lambda - the scalar argument
757: Output Parameters:
758: + A - Function matrix
759: - B - optional preconditioning matrix
761: Note:
762: `NEPComputeFunction()` is typically used within nonlinear eigensolvers
763: implementations, so most users would not generally call this routine
764: themselves.
766: Level: developer
768: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetFunction()`, `NEPComputeJacobian()`
769: @*/
770: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
771: {
772: PetscInt i;
773: PetscScalar alpha;
775: PetscFunctionBegin;
777: NEPCheckProblem(nep,1);
778: switch (nep->fui) {
779: case NEP_USER_INTERFACE_CALLBACK:
780: PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
781: PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
782: PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
783: PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
784: break;
785: case NEP_USER_INTERFACE_SPLIT:
786: PetscCall(MatZeroEntries(A));
787: if (A != B) PetscCall(MatZeroEntries(B));
788: for (i=0;i<nep->nt;i++) {
789: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
790: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
791: if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
792: }
793: break;
794: }
795: PetscFunctionReturn(PETSC_SUCCESS);
796: }
798: /*@
799: NEPComputeJacobian - Computes the Jacobian matrix $T'(\lambda)$ that has been
800: set with `NEPSetJacobian()`.
802: Collective
804: Input Parameters:
805: + nep - the nonlinear eigensolver context
806: - lambda - the scalar argument
808: Output Parameter:
809: . A - Jacobian matrix
811: Notes:
812: Most users should not need to explicitly call this routine, as it
813: is used internally within the nonlinear eigensolvers.
815: Level: developer
817: .seealso: [](ch:nep), `NEPSetJacobian()`, `NEPGetJacobian()`, `NEPComputeFunction()`
818: @*/
819: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
820: {
821: PetscInt i;
822: PetscScalar alpha;
824: PetscFunctionBegin;
826: NEPCheckProblem(nep,1);
827: switch (nep->fui) {
828: case NEP_USER_INTERFACE_CALLBACK:
829: PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
830: PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
831: PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
832: PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
833: break;
834: case NEP_USER_INTERFACE_SPLIT:
835: PetscCall(MatZeroEntries(A));
836: for (i=0;i<nep->nt;i++) {
837: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
838: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
839: }
840: break;
841: }
842: PetscFunctionReturn(PETSC_SUCCESS);
843: }