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 eigensystem.
49: Collective
51: Input Parameter:
52: . nep - the nonlinear eigensolver context
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: [](ch:nep), `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: [](ch:nep), `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: [](ch:nep), `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: [](ch:nep), `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: [](ch:nep), `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: [](ch:nep), `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: [](ch:nep), `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, see
423: `NEPConvergedReason` for the possible values
425: Options Database Key:
426: . -nep_converged_reason - print reason for convergence/divergence, and number of iterations
428: Note:
429: If this routine is called before or doing the `NEPSolve()` the value of
430: `NEP_CONVERGED_ITERATING` is returned.
432: Level: intermediate
434: .seealso: [](ch:nep), `NEPSetTolerances()`, `NEPSolve()`, `NEPConvergedReason`
435: @*/
436: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
437: {
438: PetscFunctionBegin;
440: PetscAssertPointer(reason,2);
441: NEPCheckSolved(nep,1);
442: *reason = nep->reason;
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }
446: /*@
447: NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
448: NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
450: Collective
452: Input Parameters:
453: + nep - the nonlinear eigensolver context
454: - i - index of the solution
456: Output Parameters:
457: + eigr - real part of eigenvalue
458: . eigi - imaginary part of eigenvalue
459: . Vr - real part of eigenvector
460: - Vi - imaginary part of eigenvector
462: Notes:
463: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
464: required. Otherwise, the caller must provide valid Vec objects, i.e.,
465: they must be created by the calling program with e.g. MatCreateVecs().
467: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
468: configured with complex scalars the eigenvalue is stored
469: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
470: set to zero). In any case, the user can pass NULL in Vr or Vi if one of
471: them is not required.
473: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
474: Eigenpairs are indexed according to the ordering criterion established
475: with NEPSetWhichEigenpairs().
477: Level: beginner
479: .seealso: [](ch:nep), `NEPSolve()`, `NEPGetConverged()`, `NEPSetWhichEigenpairs()`, `NEPGetLeftEigenvector()`
480: @*/
481: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
482: {
483: PetscInt k;
485: PetscFunctionBegin;
490: NEPCheckSolved(nep,1);
491: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
492: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
494: PetscCall(NEPComputeVectors(nep));
495: k = nep->perm[i];
497: /* eigenvalue */
498: #if defined(PETSC_USE_COMPLEX)
499: if (eigr) *eigr = nep->eigr[k];
500: if (eigi) *eigi = 0;
501: #else
502: if (eigr) *eigr = nep->eigr[k];
503: if (eigi) *eigi = nep->eigi[k];
504: #endif
506: /* eigenvector */
507: PetscCall(BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: /*@
512: NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
514: Collective
516: Input Parameters:
517: + nep - the nonlinear eigensolver context
518: - i - index of the solution
520: Output Parameters:
521: + Wr - real part of left eigenvector
522: - Wi - imaginary part of left eigenvector
524: Notes:
525: The caller must provide valid Vec objects, i.e., they must be created
526: by the calling program with e.g. MatCreateVecs().
528: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
529: configured with complex scalars the eigenvector is stored directly in Wr
530: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
531: them is not required.
533: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
534: Eigensolutions are indexed according to the ordering criterion established
535: with NEPSetWhichEigenpairs().
537: Left eigenvectors are available only if the twosided flag was set, see
538: NEPSetTwoSided().
540: Level: intermediate
542: .seealso: [](ch:nep), `NEPGetEigenpair()`, `NEPGetConverged()`, `NEPSetWhichEigenpairs()`, `NEPSetTwoSided()`
543: @*/
544: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
545: {
546: PetscInt k;
548: PetscFunctionBegin;
553: NEPCheckSolved(nep,1);
554: PetscCheck(nep->twosided,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
555: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
556: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
557: PetscCall(NEPComputeVectors(nep));
558: k = nep->perm[i];
559: PetscCall(BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi));
560: PetscFunctionReturn(PETSC_SUCCESS);
561: }
563: /*@
564: NEPGetErrorEstimate - Returns the error estimate associated to the i-th
565: computed eigenpair.
567: Not Collective
569: Input Parameters:
570: + nep - the nonlinear eigensolver context
571: - i - index of eigenpair
573: Output Parameter:
574: . errest - the error estimate
576: Notes:
577: This is the error estimate used internally by the eigensolver. The actual
578: error bound can be computed with NEPComputeError().
580: Level: advanced
582: .seealso: [](ch:nep), `NEPComputeError()`
583: @*/
584: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
585: {
586: PetscFunctionBegin;
588: PetscAssertPointer(errest,3);
589: NEPCheckSolved(nep,1);
590: PetscCheck(i>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
591: PetscCheck(i<nep->nconv,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see NEPGetConverged()");
592: *errest = nep->errest[nep->perm[i]];
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: /*
597: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
598: associated with an eigenpair.
600: Input Parameters:
601: adj - whether the adjoint T^* must be used instead of T
602: lambda - eigenvalue
603: x - eigenvector
604: w - array of work vectors (two vectors in split form, one vector otherwise)
605: */
606: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
607: {
608: Vec y,z=NULL;
610: PetscFunctionBegin;
611: y = w[0];
612: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
613: if (adj) PetscCall(NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL));
614: else PetscCall(NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL));
615: PetscCall(VecNorm(y,NORM_2,norm));
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*@
620: NEPComputeError - Computes the error (based on the residual norm) associated
621: with the i-th computed eigenpair.
623: Collective
625: Input Parameters:
626: + nep - the nonlinear eigensolver context
627: . i - the solution index
628: - type - the type of error to compute
630: Output Parameter:
631: . error - the error
633: Notes:
634: The error can be computed in various ways, all of them based on the residual
635: norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
636: eigenvector.
638: Level: beginner
640: .seealso: [](ch:nep), `NEPErrorType`, `NEPSolve()`, `NEPGetErrorEstimate()`
641: @*/
642: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)
643: {
644: Vec xr,xi=NULL;
645: PetscInt j,nwork,issplit=0;
646: PetscScalar kr,ki,s;
647: PetscReal er,z=0.0,errorl,nrm;
648: PetscBool flg;
650: PetscFunctionBegin;
654: PetscAssertPointer(error,4);
655: NEPCheckSolved(nep,1);
657: /* allocate work vectors */
658: #if defined(PETSC_USE_COMPLEX)
659: nwork = 2;
660: #else
661: nwork = 3;
662: #endif
663: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
664: issplit = 1;
665: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
666: }
667: PetscCall(NEPSetWorkVecs(nep,nwork));
668: xr = nep->work[issplit+1];
669: #if !defined(PETSC_USE_COMPLEX)
670: xi = nep->work[issplit+2];
671: #endif
673: /* compute residual norms */
674: PetscCall(NEPGetEigenpair(nep,i,&kr,&ki,xr,xi));
675: #if !defined(PETSC_USE_COMPLEX)
676: PetscCheck(ki==0.0,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"Not implemented for complex eigenvalues with real scalars");
677: #endif
678: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error));
679: PetscCall(VecNorm(xr,NORM_2,&er));
681: /* if two-sided, compute left residual norm and take the maximum */
682: if (nep->twosided) {
683: PetscCall(NEPGetLeftEigenvector(nep,i,xr,xi));
684: PetscCall(NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl));
685: *error = PetscMax(*error,errorl);
686: }
688: /* compute error */
689: switch (type) {
690: case NEP_ERROR_ABSOLUTE:
691: break;
692: case NEP_ERROR_RELATIVE:
693: *error /= PetscAbsScalar(kr)*er;
694: break;
695: case NEP_ERROR_BACKWARD:
696: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
697: PetscCall(NEPComputeFunction(nep,kr,nep->function,nep->function));
698: PetscCall(MatHasOperation(nep->function,MATOP_NORM,&flg));
699: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
700: PetscCall(MatNorm(nep->function,NORM_INFINITY,&nrm));
701: *error /= nrm*er;
702: break;
703: }
704: /* initialization of matrix norms */
705: if (!nep->nrma[0]) {
706: for (j=0;j<nep->nt;j++) {
707: PetscCall(MatHasOperation(nep->A[j],MATOP_NORM,&flg));
708: PetscCheck(flg,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
709: PetscCall(MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]));
710: }
711: }
712: for (j=0;j<nep->nt;j++) {
713: PetscCall(FNEvaluateFunction(nep->f[j],kr,&s));
714: z = z + nep->nrma[j]*PetscAbsScalar(s);
715: }
716: *error /= z*er;
717: break;
718: default:
719: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
720: }
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /*@
725: NEPComputeFunction - Computes the function matrix T(lambda) that has been
726: set with NEPSetFunction().
728: Collective
730: Input Parameters:
731: + nep - the nonlinear eigensolver context
732: - lambda - the scalar argument
734: Output Parameters:
735: + A - Function matrix
736: - B - optional preconditioning matrix
738: Notes:
739: NEPComputeFunction() is typically used within nonlinear eigensolvers
740: implementations, so most users would not generally call this routine
741: themselves.
743: Level: developer
745: .seealso: [](ch:nep), `NEPSetFunction()`, `NEPGetFunction()`
746: @*/
747: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
748: {
749: PetscInt i;
750: PetscScalar alpha;
752: PetscFunctionBegin;
754: NEPCheckProblem(nep,1);
755: switch (nep->fui) {
756: case NEP_USER_INTERFACE_CALLBACK:
757: PetscCheck(nep->computefunction,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
758: PetscCall(PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0));
759: PetscCallBack("NEP user Function function",(*nep->computefunction)(nep,lambda,A,B,nep->functionctx));
760: PetscCall(PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0));
761: break;
762: case NEP_USER_INTERFACE_SPLIT:
763: PetscCall(MatZeroEntries(A));
764: if (A != B) PetscCall(MatZeroEntries(B));
765: for (i=0;i<nep->nt;i++) {
766: PetscCall(FNEvaluateFunction(nep->f[i],lambda,&alpha));
767: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
768: if (A != B) PetscCall(MatAXPY(B,alpha,nep->P[i],nep->mstrp));
769: }
770: break;
771: }
772: PetscFunctionReturn(PETSC_SUCCESS);
773: }
775: /*@
776: NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
777: set with NEPSetJacobian().
779: Collective
781: Input Parameters:
782: + nep - the nonlinear eigensolver context
783: - lambda - the scalar argument
785: Output Parameter:
786: . A - Jacobian matrix
788: Notes:
789: Most users should not need to explicitly call this routine, as it
790: is used internally within the nonlinear eigensolvers.
792: Level: developer
794: .seealso: [](ch:nep), `NEPSetJacobian()`, `NEPGetJacobian()`
795: @*/
796: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
797: {
798: PetscInt i;
799: PetscScalar alpha;
801: PetscFunctionBegin;
803: NEPCheckProblem(nep,1);
804: switch (nep->fui) {
805: case NEP_USER_INTERFACE_CALLBACK:
806: PetscCheck(nep->computejacobian,PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
807: PetscCall(PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0));
808: PetscCallBack("NEP user Jacobian function",(*nep->computejacobian)(nep,lambda,A,nep->jacobianctx));
809: PetscCall(PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0));
810: break;
811: case NEP_USER_INTERFACE_SPLIT:
812: PetscCall(MatZeroEntries(A));
813: for (i=0;i<nep->nt;i++) {
814: PetscCall(FNEvaluateDerivative(nep->f[i],lambda,&alpha));
815: PetscCall(MatAXPY(A,alpha,nep->A[i],nep->mstr));
816: }
817: break;
818: }
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }