Actual source code: nepsolve.c
slepc-3.23.1 2025-05-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
10: /*
11: NEP routines related to the solution process
13: References:
15: [1] C. Campos and J.E. Roman, "NEP: a module for the parallel solution
16: of nonlinear eigenvalue problems in SLEPc", ACM Trans. Math. Soft.
17: 47(3), 23:1--23:29, 2021.
18: */
20: #include <slepc/private/nepimpl.h>
21: #include <slepc/private/bvimpl.h>
22: #include <petscdraw.h>
24: static PetscBool cited = PETSC_FALSE;
25: static const char citation[] =
26: "@Article{slepc-nep,\n"
27: " author = \"C. Campos and J. E. Roman\",\n"
28: " title = \"{NEP}: a module for the parallel solution of nonlinear eigenvalue problems in {SLEPc}\",\n"
29: " journal = \"{ACM} Trans. Math. Software\",\n"
30: " volume = \"47\",\n"
31: " number = \"3\",\n"
32: " pages = \"23:1--23:29\",\n"
33: " year = \"2021\",\n"
34: " doi = \"10.1145/3447544\"\n"
35: "}\n";
37: PetscErrorCode NEPComputeVectors(NEP nep)
38: {
39: PetscFunctionBegin;
40: NEPCheckSolved(nep,1);
41: if (nep->state==NEP_STATE_SOLVED) PetscTryTypeMethod(nep,computevectors);
42: nep->state = NEP_STATE_EIGENVECTORS;
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: /*@
47: NEPSolve - Solves the nonlinear eigensystem.
49: Collective
51: Input Parameter:
52: . nep - eigensolver context obtained from NEPCreate()
54: Options Database Keys:
55: + -nep_view - print information about the solver used
56: . -nep_view_matk - view the split form matrix Ak (replace k by an integer from 0 to nt-1)
57: . -nep_view_fnk - view the split form function fk (replace k by an integer from 0 to nt-1)
58: . -nep_view_vectors - view the computed eigenvectors
59: . -nep_view_values - view the computed eigenvalues
60: . -nep_converged_reason - print reason for convergence, and number of iterations
61: . -nep_error_absolute - print absolute errors of each eigenpair
62: . -nep_error_relative - print relative errors of each eigenpair
63: - -nep_error_backward - print backward errors of each eigenpair
65: Notes:
66: All the command-line options listed above admit an optional argument specifying
67: the viewer type and options. For instance, use '-nep_view_vectors binary:myvecs.bin'
68: to save the eigenvectors to a binary file, '-nep_view_values draw' to draw the computed
69: eigenvalues graphically, or '-nep_error_relative :myerr.m:ascii_matlab' to save
70: the errors in a file that can be executed in Matlab.
72: Level: beginner
74: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
75: @*/
76: PetscErrorCode NEPSolve(NEP nep)
77: {
78: PetscInt i;
79: char str[16];
81: PetscFunctionBegin;
83: if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
84: PetscCall(PetscCitationsRegister(citation,&cited));
85: PetscCall(PetscLogEventBegin(NEP_Solve,nep,0,0,0));
87: /* call setup */
88: PetscCall(NEPSetUp(nep));
89: nep->nconv = 0;
90: nep->its = 0;
91: for (i=0;i<nep->ncv;i++) {
92: nep->eigr[i] = 0.0;
93: nep->eigi[i] = 0.0;
94: nep->errest[i] = 0.0;
95: nep->perm[i] = i;
96: }
97: PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view_pre"));
98: PetscCall(RGViewFromOptions(nep->rg,NULL,"-rg_view"));
100: /* call solver */
101: PetscUseTypeMethod(nep,solve);
102: PetscCheck(nep->reason,PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
103: nep->state = NEP_STATE_SOLVED;
105: /* Only the first nconv columns contain useful information */
106: PetscCall(BVSetActiveColumns(nep->V,0,nep->nconv));
107: if (nep->twosided) PetscCall(BVSetActiveColumns(nep->W,0,nep->nconv));
109: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
110: PetscCall(NEPComputeVectors(nep));
111: PetscCall(NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv));
112: nep->state = NEP_STATE_EIGENVECTORS;
113: }
115: /* sort eigenvalues according to nep->which parameter */
116: PetscCall(SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm));
117: PetscCall(PetscLogEventEnd(NEP_Solve,nep,0,0,0));
119: /* various viewers */
120: PetscCall(NEPViewFromOptions(nep,NULL,"-nep_view"));
121: PetscCall(NEPConvergedReasonViewFromOptions(nep));
122: PetscCall(NEPErrorViewFromOptions(nep));
123: PetscCall(NEPValuesViewFromOptions(nep));
124: PetscCall(NEPVectorsViewFromOptions(nep));
125: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
126: for (i=0;i<nep->nt;i++) {
127: PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_mat%" PetscInt_FMT,i));
128: PetscCall(MatViewFromOptions(nep->A[i],(PetscObject)nep,str));
129: PetscCall(PetscSNPrintf(str,sizeof(str),"-nep_view_fn%" PetscInt_FMT,i));
130: PetscCall(FNViewFromOptions(nep->f[i],(PetscObject)nep,str));
131: }
132: }
134: /* Remove the initial subspace */
135: nep->nini = 0;
137: /* Reset resolvent information */
138: PetscCall(MatDestroy(&nep->resolvent));
139: PetscFunctionReturn(PETSC_SUCCESS);
140: }
142: /*@
143: NEPProjectOperator - Computes the projection of the nonlinear operator.
145: Collective
147: Input Parameters:
148: + nep - the nonlinear eigensolver context
149: . j0 - initial index
150: - j1 - final index
152: Notes:
153: This is available for split operator only.
155: The nonlinear operator T(lambda) is projected onto span(V), where V is
156: an orthonormal basis built internally by the solver. The projected
157: operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
158: computes all matrices Ei = V'*A_i*V, and stores them in the extra
159: matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
160: the previous ones are assumed to be available already.
162: Level: developer
164: .seealso: NEPSetSplitOperator()
165: @*/
166: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
167: {
168: PetscInt k;
169: Mat G;
171: PetscFunctionBegin;
175: NEPCheckProblem(nep,1);
176: NEPCheckSplit(nep,1);
177: PetscCall(BVSetActiveColumns(nep->V,j0,j1));
178: for (k=0;k<nep->nt;k++) {
179: PetscCall(DSGetMat(nep->ds,DSMatExtra[k],&G));
180: PetscCall(BVMatProject(nep->V,nep->A[k],nep->V,G));
181: PetscCall(DSRestoreMat(nep->ds,DSMatExtra[k],&G));
182: }
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@
187: NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
189: Collective
191: Input Parameters:
192: + nep - the nonlinear eigensolver context
193: . lambda - scalar argument
194: . x - vector to be multiplied against
195: - v - workspace vector (used only in the case of split form)
197: Output Parameters:
198: + y - result vector
199: . A - (optional) Function matrix, for callback interface only
200: - B - (unused) preconditioning matrix
202: Note:
203: If the nonlinear operator is represented in split form, the result
204: y = T(lambda)*x is computed without building T(lambda) explicitly. In
205: that case, parameters A and B are not used. Otherwise, the matrix
206: T(lambda) is built and the effect is the same as a call to
207: NEPComputeFunction() followed by a MatMult().
209: Level: developer
211: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
212: @*/
213: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
214: {
215: PetscInt i;
216: PetscScalar alpha;
218: PetscFunctionBegin;
227: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
228: PetscCall(VecSet(y,0.0));
229: for (i=0;i<nep->nt;i++) {
230: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
231: PetscCall(MatMult(nep->A[i],x,v));
232: PetscCall(VecAXPY(y,alpha,v));
233: }
234: } else {
235: if (!A) A = nep->function;
236: PetscCall(NEPComputeFunction(nep,lambda,A,A));
237: PetscCall(MatMult(A,x,y));
238: }
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: /*@
243: NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.
245: Collective
247: Input Parameters:
248: + nep - the nonlinear eigensolver context
249: . lambda - scalar argument
250: . x - vector to be multiplied against
251: - v - workspace vector (used only in the case of split form)
253: Output Parameters:
254: + y - result vector
255: . A - (optional) Function matrix, for callback interface only
256: - B - (unused) preconditioning matrix
258: Level: developer
260: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
261: @*/
262: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
263: {
264: PetscInt i;
265: PetscScalar alpha;
266: Vec w;
268: PetscFunctionBegin;
277: PetscCall(VecDuplicate(x,&w));
278: PetscCall(VecCopy(x,w));
279: PetscCall(VecConjugate(w));
280: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
281: PetscCall(VecSet(y,0.0));
282: for (i=0;i<nep->nt;i++) {
283: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
284: PetscCall(MatMultTranspose(nep->A[i],w,v));
285: PetscCall(VecAXPY(y,alpha,v));
286: }
287: } else {
288: if (!A) A = nep->function;
289: PetscCall(NEPComputeFunction(nep,lambda,A,A));
290: PetscCall(MatMultTranspose(A,w,y));
291: }
292: PetscCall(VecDestroy(&w));
293: PetscCall(VecConjugate(y));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*@
298: NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
300: Collective
302: Input Parameters:
303: + nep - the nonlinear eigensolver context
304: . lambda - scalar argument
305: . x - vector to be multiplied against
306: - v - workspace vector (used only in the case of split form)
308: Output Parameters:
309: + y - result vector
310: - A - (optional) Jacobian matrix, for callback interface only
312: Note:
313: If the nonlinear operator is represented in split form, the result
314: y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
315: that case, parameter A is not used. Otherwise, the matrix
316: T'(lambda) is built and the effect is the same as a call to
317: NEPComputeJacobian() followed by a MatMult().
319: Level: developer
321: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
322: @*/
323: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
324: {
325: PetscInt i;
326: PetscScalar alpha;
328: PetscFunctionBegin;
336: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
337: PetscCall(VecSet(y,0.0));
338: for (i=0;i<nep->nt;i++) {
339: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
340: PetscCall(MatMult(nep->A[i],x,v));
341: PetscCall(VecAXPY(y,alpha,v));
342: }
343: } else {
344: if (!A) A = nep->jacobian;
345: PetscCall(NEPComputeJacobian(nep,lambda,A));
346: PetscCall(MatMult(A,x,y));
347: }
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: /*@
352: NEPGetIterationNumber - Gets the current iteration number. If the
353: call to NEPSolve() is complete, then it returns the number of iterations
354: carried out by the solution method.
356: Not Collective
358: Input Parameter:
359: . nep - the nonlinear eigensolver context
361: Output Parameter:
362: . its - number of iterations
364: Note:
365: During the i-th iteration this call returns i-1. If NEPSolve() is
366: complete, then parameter "its" contains either the iteration number at
367: which convergence was successfully reached, or failure was detected.
368: Call NEPGetConvergedReason() to determine if the solver converged or
369: failed and why.
371: Level: intermediate
373: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
374: @*/
375: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
376: {
377: PetscFunctionBegin;
379: PetscAssertPointer(its,2);
380: *its = nep->its;
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*@
385: NEPGetConverged - Gets the number of converged eigenpairs.
387: Not Collective
389: Input Parameter:
390: . nep - the nonlinear eigensolver context
392: Output Parameter:
393: . nconv - number of converged eigenpairs
395: Note:
396: This function should be called after NEPSolve() has finished.
398: Level: beginner
400: .seealso: NEPSetDimensions(), NEPSolve(), NEPGetEigenpair()
401: @*/
402: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
403: {
404: PetscFunctionBegin;
406: PetscAssertPointer(nconv,2);
407: NEPCheckSolved(nep,1);
408: *nconv = nep->nconv;
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /*@
413: NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
414: stopped.
416: Not Collective
418: Input Parameter:
419: . nep - the nonlinear eigensolver context
421: Output Parameter:
422: . reason - negative value indicates diverged, positive value converged
424: Options Database Key:
425: . -nep_converged_reason - print the reason to a viewer
427: Notes:
428: Possible values for reason are
429: + NEP_CONVERGED_TOL - converged up to tolerance
430: . NEP_CONVERGED_USER - converged due to a user-defined condition
431: . NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
432: . NEP_DIVERGED_BREAKDOWN - generic breakdown in method
433: . NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
434: - NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
435: unrestarted solver
437: Can only be called after the call to NEPSolve() is complete.
439: Level: intermediate
441: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason
442: @*/
443: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
444: {
445: PetscFunctionBegin;
447: PetscAssertPointer(reason,2);
448: NEPCheckSolved(nep,1);
449: *reason = nep->reason;
450: PetscFunctionReturn(PETSC_SUCCESS);
451: }
453: /*@
454: NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
455: NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
457: Collective
459: Input Parameters:
460: + nep - nonlinear eigensolver context
461: - i - index of the solution
463: Output Parameters:
464: + eigr - real part of eigenvalue
465: . eigi - imaginary part of eigenvalue
466: . Vr - real part of eigenvector
467: - Vi - imaginary part of eigenvector
469: Notes:
470: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
471: required. Otherwise, the caller must provide valid Vec objects, i.e.,
472: they must be created by the calling program with e.g. MatCreateVecs().
474: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
475: configured with complex scalars the eigenvalue is stored
476: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
477: set to zero). In any case, the user can pass NULL in Vr or Vi if one of
478: them is not required.
480: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
481: Eigenpairs are indexed according to the ordering criterion established
482: with NEPSetWhichEigenpairs().
484: Level: beginner
486: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
487: @*/
488: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
489: {
490: PetscInt k;
492: PetscFunctionBegin;
497: NEPCheckSolved(nep,1);
498: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
499: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
501: PetscCall(NEPComputeVectors(nep));
502: k = nep->perm[i];
504: /* eigenvalue */
505: #if defined(PETSC_USE_COMPLEX)
506: if (eigr) *eigr = nep->eigr[k];
507: if (eigi) *eigi = 0;
508: #else
509: if (eigr) *eigr = nep->eigr[k];
510: if (eigi) *eigi = nep->eigi[k];
511: #endif
513: /* eigenvector */
514: PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: /*@
519: NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
521: Collective
523: Input Parameters:
524: + nep - eigensolver context
525: - i - index of the solution
527: Output Parameters:
528: + Wr - real part of left eigenvector
529: - Wi - imaginary part of left eigenvector
531: Notes:
532: The caller must provide valid Vec objects, i.e., they must be created
533: by the calling program with e.g. MatCreateVecs().
535: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
536: configured with complex scalars the eigenvector is stored directly in Wr
537: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
538: them is not required.
540: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
541: Eigensolutions are indexed according to the ordering criterion established
542: with NEPSetWhichEigenpairs().
544: Left eigenvectors are available only if the twosided flag was set, see
545: NEPSetTwoSided().
547: Level: intermediate
549: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
550: @*/
551: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
552: {
553: PetscInt k;
555: PetscFunctionBegin;
560: NEPCheckSolved(nep,1);
561: PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
562: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
563: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
564: PetscCall(NEPComputeVectors(nep));
565: k = nep->perm[i];
566: PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*@
571: NEPGetErrorEstimate - Returns the error estimate associated to the i-th
572: computed eigenpair.
574: Not Collective
576: Input Parameters:
577: + nep - nonlinear eigensolver context
578: - i - index of eigenpair
580: Output Parameter:
581: . errest - the error estimate
583: Notes:
584: This is the error estimate used internally by the eigensolver. The actual
585: error bound can be computed with NEPComputeError().
587: Level: advanced
589: .seealso: NEPComputeError()
590: @*/
591: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
592: {
593: PetscFunctionBegin;
595: PetscAssertPointer(errest,3);
596: NEPCheckSolved(nep,1);
597: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
598: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
599: *errest = nep->errest[nep->perm[i]];
600: PetscFunctionReturn(PETSC_SUCCESS);
601: }
603: /*
604: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
605: associated with an eigenpair.
607: Input Parameters:
608: adj - whether the adjoint T^* must be used instead of T
609: lambda - eigenvalue
610: x - eigenvector
611: w - array of work vectors (two vectors in split form, one vector otherwise)
612: */
613: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
614: {
615: Vec y,z=NULL;
617: PetscFunctionBegin;
618: y = w[0];
619: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
620: if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
621: else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
622: PetscCall(VecNorm(y,NORM_2,norm));
623: PetscFunctionReturn(PETSC_SUCCESS);
624: }
626: /*@
627: NEPComputeError - Computes the error (based on the residual norm) associated
628: with the i-th computed eigenpair.
630: Collective
632: Input Parameters:
633: + nep - the nonlinear eigensolver context
634: . i - the solution index
635: - type - the type of error to compute
637: Output Parameter:
638: . error - the error
640: Notes:
641: The error can be computed in various ways, all of them based on the residual
642: norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
643: eigenvector.
645: Level: beginner
647: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
648: @*/
649: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
650: {
651: Vec xr,xi=NULL;
652: PetscInt j,nwork,issplit=0;
653: PetscScalar kr,ki,s;
654: PetscReal er,z=0.0,errorl,nrm;
655: PetscBool flg;
657: PetscFunctionBegin;
661: PetscAssertPointer(error,4);
662: NEPCheckSolved(nep,1);
664: /* allocate work vectors */
665: #if defined(PETSC_USE_COMPLEX)
666: nwork = 2;
667: #else
668: nwork = 3;
669: #endif
670: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
671: issplit = 1;
672: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
673: }
674: PetscCall(NEPSetWorkVecs(nep,nwork));
675: xr = nep->work[issplit+1];
676: #if !defined(PETSC_USE_COMPLEX)
677: xi = nep->work[issplit+2];
678: #endif
680: /* compute residual norms */
681: PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
682: #if !defined(PETSC_USE_COMPLEX)
683: PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
684: #endif
685: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
686: PetscCall(VecNorm(xr,NORM_2,&er));
688: /* if two-sided, compute left residual norm and take the maximum */
689: if (nep->twosided) {
690: PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
691: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
692: *error = PetscMax(*error,errorl);
693: }
695: /* compute error */
696: switch (type) {
697: case NEP_ERROR_ABSOLUTE:
698: break;
699: case NEP_ERROR_RELATIVE:
700: *error /= PetscAbsScalar(kr)*er;
701: break;
702: case NEP_ERROR_BACKWARD:
703: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
704: PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
705: PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
706: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
707: PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
708: *error /= nrm*er;
709: break;
710: }
711: /* initialization of matrix norms */
712: if (!nep->nrma[0]) {
713: for (j=0;j<nep->nt;j++) {
714: PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
715: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
716: PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
717: }
718: }
719: for (j=0;j<nep->nt;j++) {
720: PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
721: z = z + nep->nrma[j]*PetscAbsScalar(s);
722: }
723: *error /= z*er;
724: break;
725: default:
726: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
727: }
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: /*@
732: NEPComputeFunction - Computes the function matrix T(lambda) that has been
733: set with NEPSetFunction().
735: Collective
737: Input Parameters:
738: + nep - the NEP context
739: - lambda - the scalar argument
741: Output Parameters:
742: + A - Function matrix
743: - B - optional preconditioning matrix
745: Notes:
746: NEPComputeFunction() is typically used within nonlinear eigensolvers
747: implementations, so most users would not generally call this routine
748: themselves.
750: Level: developer
752: .seealso: NEPSetFunction(), NEPGetFunction()
753: @*/
754: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
755: {
756: PetscInt i;
757: PetscScalar alpha;
759: PetscFunctionBegin;
761: NEPCheckProblem(nep,1);
762: switch (nep->fui) {
763: case NEP_USER_INTERFACE_CALLBACK:
764: PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
765: PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
766: PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
767: PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
768: break;
769: case NEP_USER_INTERFACE_SPLIT:
770: PetscCall(MatZeroEntries(A));
771: if (A != B) PetscCall(MatZeroEntries(B));
772: for (i=0;i<nep->nt;i++) {
773: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
774: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
775: if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
776: }
777: break;
778: }
779: PetscFunctionReturn(PETSC_SUCCESS);
780: }
782: /*@
783: NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
784: set with NEPSetJacobian().
786: Collective
788: Input Parameters:
789: + nep - the NEP context
790: - lambda - the scalar argument
792: Output Parameters:
793: . A - Jacobian matrix
795: Notes:
796: Most users should not need to explicitly call this routine, as it
797: is used internally within the nonlinear eigensolvers.
799: Level: developer
801: .seealso: NEPSetJacobian(), NEPGetJacobian()
802: @*/
803: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
804: {
805: PetscInt i;
806: PetscScalar alpha;
808: PetscFunctionBegin;
810: NEPCheckProblem(nep,1);
811: switch (nep->fui) {
812: case NEP_USER_INTERFACE_CALLBACK:
813: PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
814: PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
815: PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
816: PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
817: break;
818: case NEP_USER_INTERFACE_SPLIT:
819: PetscCall(MatZeroEntries(A));
820: for (i=0;i<nep->nt;i++) {
821: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
822: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
823: }
824: break;
825: }
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }