Actual source code: nepsolve.c
slepc-3.19.1 2023-05-15
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: NEP routines related to the solution process
13: References:
15: [1] C. Campos and J.E. Roman, "NEP: a module for the parallel solution
16: of nonlinear eigenvalue problems in SLEPc", ACM Trans. Math. Soft.
17: 47(3), 23:1--23:29, 2021.
18: */
20: #include <slepc/private/nepimpl.h>
21: #include <slepc/private/bvimpl.h>
22: #include <petscdraw.h>
24: static PetscBool cited = PETSC_FALSE;
25: static const char citation[] =
26: "@Article{slepc-nep,\n"
27: " author = \"C. Campos and J. E. Roman\",\n"
28: " title = \"{NEP}: a module for the parallel solution of nonlinear eigenvalue problems in {SLEPc}\",\n"
29: " journal = \"{ACM} Trans. Math. Software\",\n"
30: " volume = \"47\",\n"
31: " number = \"3\",\n"
32: " pages = \"23:1--23:29\",\n"
33: " year = \"2021\",\n"
34: " doi = \"10.1145/3447544\"\n"
35: "}\n";
37: PetscErrorCode NEPComputeVectors(NEP nep)
38: {
39: PetscFunctionBegin;
40: NEPCheckSolved(nep,1);
41: if (nep->state==NEP_STATE_SOLVED) PetscTryTypeMethod(nep,computevectors);
42: nep->state = NEP_STATE_EIGENVECTORS;
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*@
47: NEPSolve - Solves the nonlinear eigensystem.
49: Collective
51: Input Parameter:
52: . nep - eigensolver context obtained from NEPCreate()
54: Options Database Keys:
55: + -nep_view - print information about the solver used
56: . -nep_view_vectors - view the computed eigenvectors
57: . -nep_view_values - view the computed eigenvalues
58: . -nep_converged_reason - print reason for convergence, and number of iterations
59: . -nep_error_absolute - print absolute errors of each eigenpair
60: . -nep_error_relative - print relative errors of each eigenpair
61: - -nep_error_backward - print backward errors of each eigenpair
63: Notes:
64: All the command-line options listed above admit an optional argument specifying
65: the viewer type and options. For instance, use '-nep_view_vectors binary:myvecs.bin'
66: to save the eigenvectors to a binary file, '-nep_view_values draw' to draw the computed
67: eigenvalues graphically, or '-nep_error_relative :myerr.m:ascii_matlab' to save
68: the errors in a file that can be executed in Matlab.
70: Level: beginner
72: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
73: @*/
74: PetscErrorCode NEPSolve(NEP nep)
75: {
76: PetscInt i;
78: PetscFunctionBegin;
80: if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
81: PetscCall(PetscCitationsRegister(citation,&cited));
82: PetscCall(PetscLogEventBegin(NEP_Solve,nep,0,0,0));
84: /* call setup */
85: PetscCall(NEPSetUp(nep));
86: nep->nconv = 0;
87: nep->its = 0;
88: for (i=0;i<nep->ncv;i++) {
89: nep->eigr[i] = 0.0;
90: nep->eigi[i] = 0.0;
91: nep->errest[i] = 0.0;
92: nep->perm[i] = i;
93: }
94: PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view_pre"));
95: PetscCall(RGViewFromOptions(nep->rg,NULL,"-rg_view"));
97: /* call solver */
98: PetscUseTypeMethod(nep,solve);
99: PetscCheck(nep->reason,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
100: nep->state = NEP_STATE_SOLVED;
102: /* Only the first nconv columns contain useful information */
103: PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
104: if (nep->twosided) PetscCall(BVSetActiveColumns(nep->W,0,nep->nconv));
106: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
107: PetscCall(NEPComputeVectors(nep));
108: PetscCall(NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv));
109: nep->state = NEP_STATE_EIGENVECTORS;
110: }
112: /* sort eigenvalues according to nep->which parameter */
113: PetscCall(SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm));
114: PetscCall(PetscLogEventEnd(NEP_Solve,nep,0,0,0));
116: /* various viewers */
117: PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
118: PetscCall(NEPConvergedReasonViewFromOptions(nep));
119: PetscCall(NEPErrorViewFromOptions(nep));
120: PetscCall(NEPValuesViewFromOptions(nep));
121: PetscCall(NEPVectorsViewFromOptions(nep));
123: /* Remove the initial subspace */
124: nep->nini = 0;
126: /* Reset resolvent information */
127: PetscCall(MatDestroy(&nep->resolvent));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@
132: NEPProjectOperator - Computes the projection of the nonlinear operator.
134: Collective
136: Input Parameters:
137: + nep - the nonlinear eigensolver context
138: . j0 - initial index
139: - j1 - final index
141: Notes:
142: This is available for split operator only.
144: The nonlinear operator T(lambda) is projected onto span(V), where V is
145: an orthonormal basis built internally by the solver. The projected
146: operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
147: computes all matrices Ei = V'*A_i*V, and stores them in the extra
148: matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
149: the previous ones are assumed to be available already.
151: Level: developer
153: .seealso: NEPSetSplitOperator()
154: @*/
155: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
156: {
157: PetscInt k;
158: Mat G;
160: PetscFunctionBegin;
164: NEPCheckProblem(nep,1);
165: NEPCheckSplit(nep,1);
166: PetscCall(BVSetActiveColumns(nep->V,j0,j1));
167: for (k=0;k<nep->nt;k++) {
168: PetscCall(DSGetMat(nep->ds,DSMatExtra[k],&G));
169: PetscCall(BVMatProject(nep->V,nep->A[k],nep->V,G));
170: PetscCall(DSRestoreMat(nep->ds,DSMatExtra[k],&G));
171: }
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: /*@
176: NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
178: Collective
180: Input Parameters:
181: + nep - the nonlinear eigensolver context
182: . lambda - scalar argument
183: . x - vector to be multiplied against
184: - v - workspace vector (used only in the case of split form)
186: Output Parameters:
187: + y - result vector
188: . A - (optional) Function matrix, for callback interface only
189: - B - (unused) preconditioning matrix
191: Note:
192: If the nonlinear operator is represented in split form, the result
193: y = T(lambda)*x is computed without building T(lambda) explicitly. In
194: that case, parameters A and B are not used. Otherwise, the matrix
195: T(lambda) is built and the effect is the same as a call to
196: NEPComputeFunction() followed by a MatMult().
198: Level: developer
200: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
201: @*/
202: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
203: {
204: PetscInt i;
205: PetscScalar alpha;
207: PetscFunctionBegin;
216: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
217: PetscCall(VecSet(y,0.0));
218: for (i=0;i<nep->nt;i++) {
219: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
220: PetscCall(MatMult(nep->A[i],x,v));
221: PetscCall(VecAXPY(y,alpha,v));
222: }
223: } else {
224: if (!A) A = nep->function;
225: PetscCall(NEPComputeFunction(nep,lambda,A,A));
226: PetscCall(MatMult(A,x,y));
227: }
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*@
232: NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.
234: Collective
236: Input Parameters:
237: + nep - the nonlinear eigensolver context
238: . lambda - scalar argument
239: . x - vector to be multiplied against
240: - v - workspace vector (used only in the case of split form)
242: Output Parameters:
243: + y - result vector
244: . A - (optional) Function matrix, for callback interface only
245: - B - (unused) preconditioning matrix
247: Level: developer
249: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
250: @*/
251: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
252: {
253: PetscInt i;
254: PetscScalar alpha;
255: Vec w;
257: PetscFunctionBegin;
266: PetscCall(VecDuplicate(x,&w));
267: PetscCall(VecCopy(x,w));
268: PetscCall(VecConjugate(w));
269: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
270: PetscCall(VecSet(y,0.0));
271: for (i=0;i<nep->nt;i++) {
272: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
273: PetscCall(MatMultTranspose(nep->A[i],w,v));
274: PetscCall(VecAXPY(y,alpha,v));
275: }
276: } else {
277: if (!A) A = nep->function;
278: PetscCall(NEPComputeFunction(nep,lambda,A,A));
279: PetscCall(MatMultTranspose(A,w,y));
280: }
281: PetscCall(VecDestroy(&w));
282: PetscCall(VecConjugate(y));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@
287: NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
289: Collective
291: Input Parameters:
292: + nep - the nonlinear eigensolver context
293: . lambda - scalar argument
294: . x - vector to be multiplied against
295: - v - workspace vector (used only in the case of split form)
297: Output Parameters:
298: + y - result vector
299: - A - (optional) Jacobian matrix, for callback interface only
301: Note:
302: If the nonlinear operator is represented in split form, the result
303: y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
304: that case, parameter A is not used. Otherwise, the matrix
305: T'(lambda) is built and the effect is the same as a call to
306: NEPComputeJacobian() followed by a MatMult().
308: Level: developer
310: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
311: @*/
312: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
313: {
314: PetscInt i;
315: PetscScalar alpha;
317: PetscFunctionBegin;
325: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
326: PetscCall(VecSet(y,0.0));
327: for (i=0;i<nep->nt;i++) {
328: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
329: PetscCall(MatMult(nep->A[i],x,v));
330: PetscCall(VecAXPY(y,alpha,v));
331: }
332: } else {
333: if (!A) A = nep->jacobian;
334: PetscCall(NEPComputeJacobian(nep,lambda,A));
335: PetscCall(MatMult(A,x,y));
336: }
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }
340: /*@
341: NEPGetIterationNumber - Gets the current iteration number. If the
342: call to NEPSolve() is complete, then it returns the number of iterations
343: carried out by the solution method.
345: Not Collective
347: Input Parameter:
348: . nep - the nonlinear eigensolver context
350: Output Parameter:
351: . its - number of iterations
353: Note:
354: During the i-th iteration this call returns i-1. If NEPSolve() is
355: complete, then parameter "its" contains either the iteration number at
356: which convergence was successfully reached, or failure was detected.
357: Call NEPGetConvergedReason() to determine if the solver converged or
358: failed and why.
360: Level: intermediate
362: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
363: @*/
364: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
365: {
366: PetscFunctionBegin;
369: *its = nep->its;
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: /*@
374: NEPGetConverged - Gets the number of converged eigenpairs.
376: Not Collective
378: Input Parameter:
379: . nep - the nonlinear eigensolver context
381: Output Parameter:
382: . nconv - number of converged eigenpairs
384: Note:
385: This function should be called after NEPSolve() has finished.
387: Level: beginner
389: .seealso: NEPSetDimensions(), NEPSolve(), NEPGetEigenpair()
390: @*/
391: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
392: {
393: PetscFunctionBegin;
396: NEPCheckSolved(nep,1);
397: *nconv = nep->nconv;
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
401: /*@
402: NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
403: stopped.
405: Not Collective
407: Input Parameter:
408: . nep - the nonlinear eigensolver context
410: Output Parameter:
411: . reason - negative value indicates diverged, positive value converged
413: Options Database Key:
414: . -nep_converged_reason - print the reason to a viewer
416: Notes:
417: Possible values for reason are
418: + NEP_CONVERGED_TOL - converged up to tolerance
419: . NEP_CONVERGED_USER - converged due to a user-defined condition
420: . NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
421: . NEP_DIVERGED_BREAKDOWN - generic breakdown in method
422: . NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
423: - NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
424: unrestarted solver
426: Can only be called after the call to NEPSolve() is complete.
428: Level: intermediate
430: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason
431: @*/
432: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
433: {
434: PetscFunctionBegin;
437: NEPCheckSolved(nep,1);
438: *reason = nep->reason;
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*@C
443: NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
444: NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
446: Collective
448: Input Parameters:
449: + nep - nonlinear eigensolver context
450: - i - index of the solution
452: Output Parameters:
453: + eigr - real part of eigenvalue
454: . eigi - imaginary part of eigenvalue
455: . Vr - real part of eigenvector
456: - Vi - imaginary part of eigenvector
458: Notes:
459: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
460: required. Otherwise, the caller must provide valid Vec objects, i.e.,
461: they must be created by the calling program with e.g. MatCreateVecs().
463: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
464: configured with complex scalars the eigenvalue is stored
465: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
466: set to zero). In any case, the user can pass NULL in Vr or Vi if one of
467: them is not required.
469: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
470: Eigenpairs are indexed according to the ordering criterion established
471: with NEPSetWhichEigenpairs().
473: Level: beginner
475: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
476: @*/
477: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
478: {
479: PetscInt k;
481: PetscFunctionBegin;
486: NEPCheckSolved(nep,1);
487: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
488: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
490: PetscCall(NEPComputeVectors(nep));
491: k = nep->perm[i];
493: /* eigenvalue */
494: #if defined(PETSC_USE_COMPLEX)
495: if (eigr) *eigr = nep->eigr[k];
496: if (eigi) *eigi = 0;
497: #else
498: if (eigr) *eigr = nep->eigr[k];
499: if (eigi) *eigi = nep->eigi[k];
500: #endif
502: /* eigenvector */
503: PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: /*@
508: NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
510: Collective
512: Input Parameters:
513: + nep - eigensolver context
514: - i - index of the solution
516: Output Parameters:
517: + Wr - real part of left eigenvector
518: - Wi - imaginary part of left eigenvector
520: Notes:
521: The caller must provide valid Vec objects, i.e., they must be created
522: by the calling program with e.g. MatCreateVecs().
524: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
525: configured with complex scalars the eigenvector is stored directly in Wr
526: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
527: them is not required.
529: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
530: Eigensolutions are indexed according to the ordering criterion established
531: with NEPSetWhichEigenpairs().
533: Left eigenvectors are available only if the twosided flag was set, see
534: NEPSetTwoSided().
536: Level: intermediate
538: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
539: @*/
540: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
541: {
542: PetscInt k;
544: PetscFunctionBegin;
549: NEPCheckSolved(nep,1);
550: PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
551: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
552: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
553: PetscCall(NEPComputeVectors(nep));
554: k = nep->perm[i];
555: PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: /*@
560: NEPGetErrorEstimate - Returns the error estimate associated to the i-th
561: computed eigenpair.
563: Not Collective
565: Input Parameters:
566: + nep - nonlinear eigensolver context
567: - i - index of eigenpair
569: Output Parameter:
570: . errest - the error estimate
572: Notes:
573: This is the error estimate used internally by the eigensolver. The actual
574: error bound can be computed with NEPComputeError().
576: Level: advanced
578: .seealso: NEPComputeError()
579: @*/
580: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
581: {
582: PetscFunctionBegin;
585: NEPCheckSolved(nep,1);
586: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
587: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
588: *errest = nep->errest[nep->perm[i]];
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*
593: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
594: associated with an eigenpair.
596: Input Parameters:
597: adj - whether the adjoint T^* must be used instead of T
598: lambda - eigenvalue
599: x - eigenvector
600: w - array of work vectors (two vectors in split form, one vector otherwise)
601: */
602: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
603: {
604: Vec y,z=NULL;
606: PetscFunctionBegin;
607: y = w[0];
608: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
609: if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
610: else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
611: PetscCall(VecNorm(y,NORM_2,norm));
612: PetscFunctionReturn(PETSC_SUCCESS);
613: }
615: /*@
616: NEPComputeError - Computes the error (based on the residual norm) associated
617: with the i-th computed eigenpair.
619: Collective
621: Input Parameters:
622: + nep - the nonlinear eigensolver context
623: . i - the solution index
624: - type - the type of error to compute
626: Output Parameter:
627: . error - the error
629: Notes:
630: The error can be computed in various ways, all of them based on the residual
631: norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
632: eigenvector.
634: Level: beginner
636: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
637: @*/
638: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
639: {
640: Vec xr,xi=NULL;
641: PetscInt j,nwork,issplit=0;
642: PetscScalar kr,ki,s;
643: PetscReal er,z=0.0,errorl,nrm;
644: PetscBool flg;
646: PetscFunctionBegin;
651: NEPCheckSolved(nep,1);
653: /* allocate work vectors */
654: #if defined(PETSC_USE_COMPLEX)
655: nwork = 2;
656: #else
657: nwork = 3;
658: #endif
659: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
660: issplit = 1;
661: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
662: }
663: PetscCall(NEPSetWorkVecs(nep,nwork));
664: xr = nep->work[issplit+1];
665: #if !defined(PETSC_USE_COMPLEX)
666: xi = nep->work[issplit+2];
667: #endif
669: /* compute residual norms */
670: PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
671: #if !defined(PETSC_USE_COMPLEX)
672: PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
673: #endif
674: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
675: PetscCall(VecNorm(xr,NORM_2,&er));
677: /* if two-sided, compute left residual norm and take the maximum */
678: if (nep->twosided) {
679: PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
680: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
681: *error = PetscMax(*error,errorl);
682: }
684: /* compute error */
685: switch (type) {
686: case NEP_ERROR_ABSOLUTE:
687: break;
688: case NEP_ERROR_RELATIVE:
689: *error /= PetscAbsScalar(kr)*er;
690: break;
691: case NEP_ERROR_BACKWARD:
692: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
693: PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
694: PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
695: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
696: PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
697: *error /= nrm*er;
698: break;
699: }
700: /* initialization of matrix norms */
701: if (!nep->nrma[0]) {
702: for (j=0;j<nep->nt;j++) {
703: PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
704: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
705: PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
706: }
707: }
708: for (j=0;j<nep->nt;j++) {
709: PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
710: z = z + nep->nrma[j]*PetscAbsScalar(s);
711: }
712: *error /= z*er;
713: break;
714: default:
715: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
716: }
717: PetscFunctionReturn(PETSC_SUCCESS);
718: }
720: /*@
721: NEPComputeFunction - Computes the function matrix T(lambda) that has been
722: set with NEPSetFunction().
724: Collective
726: Input Parameters:
727: + nep - the NEP context
728: - lambda - the scalar argument
730: Output Parameters:
731: + A - Function matrix
732: - B - optional preconditioning matrix
734: Notes:
735: NEPComputeFunction() is typically used within nonlinear eigensolvers
736: implementations, so most users would not generally call this routine
737: themselves.
739: Level: developer
741: .seealso: NEPSetFunction(), NEPGetFunction()
742: @*/
743: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
744: {
745: PetscInt i;
746: PetscScalar alpha;
748: PetscFunctionBegin;
750: NEPCheckProblem(nep,1);
751: switch (nep->fui) {
752: case NEP_USER_INTERFACE_CALLBACK:
753: PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
754: PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
755: PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
756: PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
757: break;
758: case NEP_USER_INTERFACE_SPLIT:
759: PetscCall(MatZeroEntries(A));
760: if (A != B) PetscCall(MatZeroEntries(B));
761: for (i=0;i<nep->nt;i++) {
762: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
763: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
764: if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
765: }
766: break;
767: }
768: PetscFunctionReturn(PETSC_SUCCESS);
769: }
771: /*@
772: NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
773: set with NEPSetJacobian().
775: Collective
777: Input Parameters:
778: + nep - the NEP context
779: - lambda - the scalar argument
781: Output Parameters:
782: . A - Jacobian matrix
784: Notes:
785: Most users should not need to explicitly call this routine, as it
786: is used internally within the nonlinear eigensolvers.
788: Level: developer
790: .seealso: NEPSetJacobian(), NEPGetJacobian()
791: @*/
792: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
793: {
794: PetscInt i;
795: PetscScalar alpha;
797: PetscFunctionBegin;
799: NEPCheckProblem(nep,1);
800: switch (nep->fui) {
801: case NEP_USER_INTERFACE_CALLBACK:
802: PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
803: PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
804: PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
805: PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
806: break;
807: case NEP_USER_INTERFACE_SPLIT:
808: PetscCall(MatZeroEntries(A));
809: for (i=0;i<nep->nt;i++) {
810: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
811: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
812: }
813: break;
814: }
815: PetscFunctionReturn(PETSC_SUCCESS);
816: }