1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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", arXiv:1910.11712, 2019.
17: */
19: #include <slepc/private/nepimpl.h> 20: #include <slepc/private/bvimpl.h> 21: #include <petscdraw.h>
23: static PetscBool cited = PETSC_FALSE;
24: static const char citation[] =
25: "@Article{slepc-nep,\n"
26: " author = \"C. Campos and J. E. Roman\",\n"
27: " title = \"{NEP}: a module for the parallel solution of nonlinear eigenvalue problems in {SLEPc}\",\n"
28: " journal = \"arXiv:1910.11712\",\n"
29: " url = \"https://arxiv.org/abs/1910.11712\",\n"
30: " year = \"2019\"\n"
31: "}\n";
33: PetscErrorCode NEPComputeVectors(NEP nep) 34: {
38: NEPCheckSolved(nep,1);
39: if (nep->state==NEP_STATE_SOLVED && nep->ops->computevectors) {
40: (*nep->ops->computevectors)(nep);
41: }
42: nep->state = NEP_STATE_EIGENVECTORS;
43: return(0);
44: }
46: /*@
47: NEPSolve - Solves the nonlinear eigensystem.
49: Collective on nep
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 binary - save the computed eigenvectors to the default binary viewer
57: . -nep_view_values - print 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
62: Level: beginner
64: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
65: @*/
66: PetscErrorCode NEPSolve(NEP nep) 67: {
69: PetscInt i;
73: if (nep->state>=NEP_STATE_SOLVED) return(0);
74: PetscCitationsRegister(citation,&cited);
75: PetscLogEventBegin(NEP_Solve,nep,0,0,0);
77: /* call setup */
78: NEPSetUp(nep);
79: nep->nconv = 0;
80: nep->its = 0;
81: for (i=0;i<nep->ncv;i++) {
82: nep->eigr[i] = 0.0;
83: nep->eigi[i] = 0.0;
84: nep->errest[i] = 0.0;
85: nep->perm[i] = i;
86: }
87: NEPViewFromOptions(nep,NULL,"-nep_view_pre");
88: RGViewFromOptions(nep->rg,NULL,"-rg_view");
90: /* call solver */
91: (*nep->ops->solve)(nep);
92: if (!nep->reason) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
93: nep->state = NEP_STATE_SOLVED;
95: /* Only the first nconv columns contain useful information */
96: BVSetActiveColumns(nep->V,0,nep->nconv);
97: if (nep->twosided) { BVSetActiveColumns(nep->W,0,nep->nconv); }
99: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
100: NEPComputeVectors(nep);
101: NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv);
102: nep->state = NEP_STATE_EIGENVECTORS;
103: }
105: /* sort eigenvalues according to nep->which parameter */
106: SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm);
107: PetscLogEventEnd(NEP_Solve,nep,0,0,0);
109: /* various viewers */
110: NEPViewFromOptions(nep,NULL,"-nep_view");
111: NEPConvergedReasonViewFromOptions(nep);
112: NEPErrorViewFromOptions(nep);
113: NEPValuesViewFromOptions(nep);
114: NEPVectorsViewFromOptions(nep);
116: /* Remove the initial subspace */
117: nep->nini = 0;
119: /* Reset resolvent information */
120: MatDestroy(&nep->resolvent);
121: return(0);
122: }
124: /*@
125: NEPProjectOperator - Computes the projection of the nonlinear operator.
127: Collective on nep
129: Input Parameters:
130: + nep - the nonlinear eigensolver context
131: . j0 - initial index
132: - j1 - final index
134: Notes:
135: This is available for split operator only.
137: The nonlinear operator T(lambda) is projected onto span(V), where V is
138: an orthonormal basis built internally by the solver. The projected
139: operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
140: computes all matrices Ei = V'*A_i*V, and stores them in the extra
141: matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
142: the previous ones are assumed to be available already.
144: Level: developer
146: .seealso: NEPSetSplitOperator()
147: @*/
148: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)149: {
151: PetscInt k;
152: Mat G;
158: NEPCheckProblem(nep,1);
159: NEPCheckSplit(nep,1);
160: BVSetActiveColumns(nep->V,j0,j1);
161: for (k=0;k<nep->nt;k++) {
162: DSGetMat(nep->ds,DSMatExtra[k],&G);
163: BVMatProject(nep->V,nep->A[k],nep->V,G);
164: DSRestoreMat(nep->ds,DSMatExtra[k],&G);
165: }
166: return(0);
167: }
169: /*@
170: NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
172: Collective on nep
174: Input Parameters:
175: + nep - the nonlinear eigensolver context
176: . lambda - scalar argument
177: . x - vector to be multiplied against
178: - v - workspace vector (used only in the case of split form)
180: Output Parameters:
181: + y - result vector
182: . A - Function matrix
183: - B - optional preconditioning matrix
185: Note:
186: If the nonlinear operator is represented in split form, the result
187: y = T(lambda)*x is computed without building T(lambda) explicitly. In
188: that case, parameters A and B are not used. Otherwise, the matrix
189: T(lambda) is built and the effect is the same as a call to
190: NEPComputeFunction() followed by a MatMult().
192: Level: developer
194: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
195: @*/
196: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)197: {
199: PetscInt i;
200: PetscScalar alpha;
211: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
212: VecSet(y,0.0);
213: for (i=0;i<nep->nt;i++) {
214: FNEvaluateFunction(nep->f[i],lambda,&alpha);
215: MatMult(nep->A[i],x,v);
216: VecAXPY(y,alpha,v);
217: }
218: } else {
219: NEPComputeFunction(nep,lambda,A,B);
220: MatMult(A,x,y);
221: }
222: return(0);
223: }
225: /*@
226: NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.
228: Collective on nep
230: Input Parameters:
231: + nep - the nonlinear eigensolver context
232: . lambda - scalar argument
233: . x - vector to be multiplied against
234: - v - workspace vector (used only in the case of split form)
236: Output Parameters:
237: + y - result vector
238: . A - Function matrix
239: - B - optional preconditioning matrix
241: Level: developer
243: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
244: @*/
245: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)246: {
248: PetscInt i;
249: PetscScalar alpha;
260: VecConjugate(x);
261: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
262: VecSet(y,0.0);
263: for (i=0;i<nep->nt;i++) {
264: FNEvaluateFunction(nep->f[i],lambda,&alpha);
265: MatMultTranspose(nep->A[i],x,v);
266: VecAXPY(y,alpha,v);
267: }
268: } else {
269: NEPComputeFunction(nep,lambda,A,B);
270: MatMultTranspose(A,x,y);
271: }
272: VecConjugate(x);
273: VecConjugate(y);
274: return(0);
275: }
277: /*@
278: NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
280: Collective on nep
282: Input Parameters:
283: + nep - the nonlinear eigensolver context
284: . lambda - scalar argument
285: . x - vector to be multiplied against
286: - v - workspace vector (used only in the case of split form)
288: Output Parameters:
289: + y - result vector
290: - A - Jacobian matrix
292: Note:
293: If the nonlinear operator is represented in split form, the result
294: y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
295: that case, parameter A is not used. Otherwise, the matrix
296: T'(lambda) is built and the effect is the same as a call to
297: NEPComputeJacobian() followed by a MatMult().
299: Level: developer
301: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
302: @*/
303: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)304: {
306: PetscInt i;
307: PetscScalar alpha;
317: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
318: VecSet(y,0.0);
319: for (i=0;i<nep->nt;i++) {
320: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
321: MatMult(nep->A[i],x,v);
322: VecAXPY(y,alpha,v);
323: }
324: } else {
325: NEPComputeJacobian(nep,lambda,A);
326: MatMult(A,x,y);
327: }
328: return(0);
329: }
331: /*@
332: NEPGetIterationNumber - Gets the current iteration number. If the
333: call to NEPSolve() is complete, then it returns the number of iterations
334: carried out by the solution method.
336: Not Collective
338: Input Parameter:
339: . nep - the nonlinear eigensolver context
341: Output Parameter:
342: . its - number of iterations
344: Note:
345: During the i-th iteration this call returns i-1. If NEPSolve() is
346: complete, then parameter "its" contains either the iteration number at
347: which convergence was successfully reached, or failure was detected.
348: Call NEPGetConvergedReason() to determine if the solver converged or
349: failed and why.
351: Level: intermediate
353: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
354: @*/
355: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)356: {
360: *its = nep->its;
361: return(0);
362: }
364: /*@
365: NEPGetConverged - Gets the number of converged eigenpairs.
367: Not Collective
369: Input Parameter:
370: . nep - the nonlinear eigensolver context
372: Output Parameter:
373: . nconv - number of converged eigenpairs
375: Note:
376: This function should be called after NEPSolve() has finished.
378: Level: beginner
380: .seealso: NEPSetDimensions(), NEPSolve()
381: @*/
382: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)383: {
387: NEPCheckSolved(nep,1);
388: *nconv = nep->nconv;
389: return(0);
390: }
392: /*@
393: NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
394: stopped.
396: Not Collective
398: Input Parameter:
399: . nep - the nonlinear eigensolver context
401: Output Parameter:
402: . reason - negative value indicates diverged, positive value converged
404: Notes:
406: Possible values for reason are
407: + NEP_CONVERGED_TOL - converged up to tolerance
408: . NEP_CONVERGED_USER - converged due to a user-defined condition
409: . NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
410: . NEP_DIVERGED_BREAKDOWN - generic breakdown in method
411: . NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
412: - NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
413: unrestarted solver
415: Can only be called after the call to NEPSolve() is complete.
417: Level: intermediate
419: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason420: @*/
421: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)422: {
426: NEPCheckSolved(nep,1);
427: *reason = nep->reason;
428: return(0);
429: }
431: /*@C
432: NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
433: NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
435: Logically Collective on nep
437: Input Parameters:
438: + nep - nonlinear eigensolver context
439: - i - index of the solution
441: Output Parameters:
442: + eigr - real part of eigenvalue
443: . eigi - imaginary part of eigenvalue
444: . Vr - real part of eigenvector
445: - Vi - imaginary part of eigenvector
447: Notes:
448: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
449: required. Otherwise, the caller must provide valid Vec objects, i.e.,
450: they must be created by the calling program with e.g. MatCreateVecs().
452: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
453: configured with complex scalars the eigenvalue is stored
454: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
455: set to zero). In any case, the user can pass NULL in Vr or Vi if one of
456: them is not required.
458: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
459: Eigenpairs are indexed according to the ordering criterion established
460: with NEPSetWhichEigenpairs().
462: Level: beginner
464: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
465: @*/
466: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)467: {
468: PetscInt k;
476: NEPCheckSolved(nep,1);
477: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
479: NEPComputeVectors(nep);
480: k = nep->perm[i];
482: /* eigenvalue */
483: #if defined(PETSC_USE_COMPLEX)
484: if (eigr) *eigr = nep->eigr[k];
485: if (eigi) *eigi = 0;
486: #else
487: if (eigr) *eigr = nep->eigr[k];
488: if (eigi) *eigi = nep->eigi[k];
489: #endif
491: /* eigenvector */
492: BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi);
493: return(0);
494: }
496: /*@
497: NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
499: Logically Collective on nep
501: Input Parameters:
502: + nep - eigensolver context
503: - i - index of the solution
505: Output Parameters:
506: + Wr - real part of left eigenvector
507: - Wi - imaginary part of left eigenvector
509: Notes:
510: The caller must provide valid Vec objects, i.e., they must be created
511: by the calling program with e.g. MatCreateVecs().
513: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
514: configured with complex scalars the eigenvector is stored directly in Wr
515: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
516: them is not required.
518: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
519: Eigensolutions are indexed according to the ordering criterion established
520: with NEPSetWhichEigenpairs().
522: Left eigenvectors are available only if the twosided flag was set, see
523: NEPSetTwoSided().
525: Level: intermediate
527: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
528: @*/
529: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)530: {
532: PetscInt k;
539: NEPCheckSolved(nep,1);
540: if (!nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
541: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
542: NEPComputeVectors(nep);
543: k = nep->perm[i];
544: BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi);
545: return(0);
546: }
548: /*@
549: NEPGetErrorEstimate - Returns the error estimate associated to the i-th
550: computed eigenpair.
552: Not Collective
554: Input Parameter:
555: + nep - nonlinear eigensolver context
556: - i - index of eigenpair
558: Output Parameter:
559: . errest - the error estimate
561: Notes:
562: This is the error estimate used internally by the eigensolver. The actual
563: error bound can be computed with NEPComputeRelativeError().
565: Level: advanced
567: .seealso: NEPComputeRelativeError()
568: @*/
569: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)570: {
574: NEPCheckSolved(nep,1);
575: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
576: *errest = nep->errest[nep->perm[i]];
577: return(0);
578: }
580: /*
581: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
582: associated with an eigenpair.
584: Input Parameters:
585: adj - whether the adjoint T^* must be used instead of T
586: lambda - eigenvalue
587: x - eigenvector
588: w - array of work vectors (two vectors in split form, one vector otherwise)
589: */
590: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)591: {
593: Vec y,z=NULL;
596: y = w[0];
597: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
598: if (adj) {
599: NEPApplyAdjoint(nep,lambda,x,z,y,nep->function,nep->function_pre);
600: } else {
601: NEPApplyFunction(nep,lambda,x,z,y,nep->function,nep->function_pre);
602: }
603: VecNorm(y,NORM_2,norm);
604: return(0);
605: }
607: /*@
608: NEPComputeError - Computes the error (based on the residual norm) associated
609: with the i-th computed eigenpair.
611: Collective on nep
613: Input Parameter:
614: + nep - the nonlinear eigensolver context
615: . i - the solution index
616: - type - the type of error to compute
618: Output Parameter:
619: . error - the error
621: Notes:
622: The error can be computed in various ways, all of them based on the residual
623: norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
624: eigenvector.
626: Level: beginner
628: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
629: @*/
630: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)631: {
633: Vec xr,xi=NULL;
634: PetscInt j,nwork,issplit=0;
635: PetscScalar kr,ki,s;
636: PetscReal er,z=0.0,errorl,nrm;
637: PetscBool flg;
644: NEPCheckSolved(nep,1);
646: /* allocate work vectors */
647: #if defined(PETSC_USE_COMPLEX)
648: nwork = 2;
649: #else
650: nwork = 3;
651: #endif
652: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
653: issplit = 1;
654: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
655: }
656: NEPSetWorkVecs(nep,nwork);
657: xr = nep->work[issplit+1];
658: #if !defined(PETSC_USE_COMPLEX)
659: xi = nep->work[issplit+2];
660: #endif
662: /* compute residual norms */
663: NEPGetEigenpair(nep,i,&kr,&ki,xr,xi);
664: #if !defined(PETSC_USE_COMPLEX)
665: if (ki) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented for complex eigenvalues with real scalars");
666: #endif
667: NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error);
668: VecNorm(xr,NORM_2,&er);
670: /* if two-sided, compute left residual norm and take the maximum */
671: if (nep->twosided) {
672: NEPGetLeftEigenvector(nep,i,xr,xi);
673: NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl);
674: *error = PetscMax(*error,errorl);
675: }
677: /* compute error */
678: switch (type) {
679: case NEP_ERROR_ABSOLUTE:
680: break;
681: case NEP_ERROR_RELATIVE:
682: *error /= PetscAbsScalar(kr)*er;
683: break;
684: case NEP_ERROR_BACKWARD:
685: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
686: NEPComputeFunction(nep,kr,nep->function,nep->function);
687: MatHasOperation(nep->function,MATOP_NORM,&flg);
688: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
689: MatNorm(nep->function,NORM_INFINITY,&nrm);
690: *error /= nrm*er;
691: break;
692: }
693: /* initialization of matrix norms */
694: if (!nep->nrma[0]) {
695: for (j=0;j<nep->nt;j++) {
696: MatHasOperation(nep->A[j],MATOP_NORM,&flg);
697: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
698: MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]);
699: }
700: }
701: for (j=0;j<nep->nt;j++) {
702: FNEvaluateFunction(nep->f[j],kr,&s);
703: z = z + nep->nrma[j]*PetscAbsScalar(s);
704: }
705: *error /= z*er;
706: break;
707: default:708: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
709: }
710: return(0);
711: }
713: /*@
714: NEPComputeFunction - Computes the function matrix T(lambda) that has been
715: set with NEPSetFunction().
717: Collective on nep
719: Input Parameters:
720: + nep - the NEP context
721: - lambda - the scalar argument
723: Output Parameters:
724: + A - Function matrix
725: - B - optional preconditioning matrix
727: Notes:
728: NEPComputeFunction() is typically used within nonlinear eigensolvers
729: implementations, so most users would not generally call this routine
730: themselves.
732: Level: developer
734: .seealso: NEPSetFunction(), NEPGetFunction()
735: @*/
736: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)737: {
739: PetscInt i;
740: PetscScalar alpha;
744: NEPCheckProblem(nep,1);
745: switch (nep->fui) {
746: case NEP_USER_INTERFACE_CALLBACK:
747: if (!nep->computefunction) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
748: PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0);
749: PetscStackPush("NEP user Function function");
750: (*nep->computefunction)(nep,lambda,A,B,nep->functionctx);
751: PetscStackPop;
752: PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0);
753: break;
754: case NEP_USER_INTERFACE_SPLIT:
755: MatZeroEntries(A);
756: for (i=0;i<nep->nt;i++) {
757: FNEvaluateFunction(nep->f[i],lambda,&alpha);
758: MatAXPY(A,alpha,nep->A[i],nep->mstr);
759: }
760: if (A != B) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented");
761: break;
762: }
763: return(0);
764: }
766: /*@
767: NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
768: set with NEPSetJacobian().
770: Collective on nep
772: Input Parameters:
773: + nep - the NEP context
774: - lambda - the scalar argument
776: Output Parameters:
777: . A - Jacobian matrix
779: Notes:
780: Most users should not need to explicitly call this routine, as it
781: is used internally within the nonlinear eigensolvers.
783: Level: developer
785: .seealso: NEPSetJacobian(), NEPGetJacobian()
786: @*/
787: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)788: {
790: PetscInt i;
791: PetscScalar alpha;
795: NEPCheckProblem(nep,1);
796: switch (nep->fui) {
797: case NEP_USER_INTERFACE_CALLBACK:
798: if (!nep->computejacobian) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
799: PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0);
800: PetscStackPush("NEP user Jacobian function");
801: (*nep->computejacobian)(nep,lambda,A,nep->jacobianctx);
802: PetscStackPop;
803: PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0);
804: break;
805: case NEP_USER_INTERFACE_SPLIT:
806: MatZeroEntries(A);
807: for (i=0;i<nep->nt;i++) {
808: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
809: MatAXPY(A,alpha,nep->A[i],nep->mstr);
810: }
811: break;
812: }
813: return(0);
814: }