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
12: */
14: #include <slepc/private/nepimpl.h> 15: #include <slepc/private/bvimpl.h> 16: #include <petscdraw.h>
18: static PetscBool cited = PETSC_FALSE;
19: static const char citation[] =
20: "@Article{slepc-nep,\n"
21: " author = \"C. Campos and J. E. Roman\",\n"
22: " title = \"{NEP}: a module for the parallel solution of nonlinear eigenvalue problems in {SLEPc}\",\n"
23: " journal = \"arXiv:1910.11712\",\n"
24: " url = \"https://arxiv.org/abs/1910.11712\",\n"
25: " year = \"2019\"\n"
26: "}\n";
28: PetscErrorCode NEPComputeVectors(NEP nep) 29: {
33: NEPCheckSolved(nep,1);
34: if (nep->state==NEP_STATE_SOLVED && nep->ops->computevectors) {
35: (*nep->ops->computevectors)(nep);
36: }
37: nep->state = NEP_STATE_EIGENVECTORS;
38: return(0);
39: }
41: /*@
42: NEPSolve - Solves the nonlinear eigensystem.
44: Collective on nep
46: Input Parameter:
47: . nep - eigensolver context obtained from NEPCreate()
49: Options Database Keys:
50: + -nep_view - print information about the solver used
51: . -nep_view_vectors binary - save the computed eigenvectors to the default binary viewer
52: . -nep_view_values - print computed eigenvalues
53: . -nep_converged_reason - print reason for convergence, and number of iterations
54: . -nep_error_absolute - print absolute errors of each eigenpair
55: - -nep_error_relative - print relative errors of each eigenpair
57: Level: beginner
59: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
60: @*/
61: PetscErrorCode NEPSolve(NEP nep) 62: {
64: PetscInt i;
68: if (nep->state>=NEP_STATE_SOLVED) return(0);
69: PetscCitationsRegister(citation,&cited);
70: PetscLogEventBegin(NEP_Solve,nep,0,0,0);
72: /* call setup */
73: NEPSetUp(nep);
74: nep->nconv = 0;
75: nep->its = 0;
76: for (i=0;i<nep->ncv;i++) {
77: nep->eigr[i] = 0.0;
78: nep->eigi[i] = 0.0;
79: nep->errest[i] = 0.0;
80: nep->perm[i] = i;
81: }
82: NEPViewFromOptions(nep,NULL,"-nep_view_pre");
83: RGViewFromOptions(nep->rg,NULL,"-rg_view");
85: /* call solver */
86: (*nep->ops->solve)(nep);
87: if (!nep->reason) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
88: nep->state = NEP_STATE_SOLVED;
90: /* Only the first nconv columns contain useful information */
91: BVSetActiveColumns(nep->V,0,nep->nconv);
92: if (nep->twosided) { BVSetActiveColumns(nep->W,0,nep->nconv); }
94: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
95: NEPComputeVectors(nep);
96: NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv);
97: nep->state = NEP_STATE_EIGENVECTORS;
98: }
100: /* sort eigenvalues according to nep->which parameter */
101: SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm);
102: PetscLogEventEnd(NEP_Solve,nep,0,0,0);
104: /* various viewers */
105: NEPViewFromOptions(nep,NULL,"-nep_view");
106: NEPConvergedReasonViewFromOptions(nep);
107: NEPErrorViewFromOptions(nep);
108: NEPValuesViewFromOptions(nep);
109: NEPVectorsViewFromOptions(nep);
111: /* Remove the initial subspace */
112: nep->nini = 0;
114: /* Reset resolvent information */
115: MatDestroy(&nep->resolvent);
116: return(0);
117: }
119: /*@
120: NEPProjectOperator - Computes the projection of the nonlinear operator.
122: Collective on nep
124: Input Parameters:
125: + nep - the nonlinear eigensolver context
126: . j0 - initial index
127: - j1 - final index
129: Notes:
130: This is available for split operator only.
132: The nonlinear operator T(lambda) is projected onto span(V), where V is
133: an orthonormal basis built internally by the solver. The projected
134: operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
135: computes all matrices Ei = V'*A_i*V, and stores them in the extra
136: matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
137: the previous ones are assumed to be available already.
139: Level: developer
141: .seealso: NEPSetSplitOperator()
142: @*/
143: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)144: {
146: PetscInt k;
147: Mat G;
153: NEPCheckProblem(nep,1);
154: NEPCheckSplit(nep,1);
155: BVSetActiveColumns(nep->V,j0,j1);
156: for (k=0;k<nep->nt;k++) {
157: DSGetMat(nep->ds,DSMatExtra[k],&G);
158: BVMatProject(nep->V,nep->A[k],nep->V,G);
159: DSRestoreMat(nep->ds,DSMatExtra[k],&G);
160: }
161: return(0);
162: }
164: /*@
165: NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
167: Collective on nep
169: Input Parameters:
170: + nep - the nonlinear eigensolver context
171: . lambda - scalar argument
172: . x - vector to be multiplied against
173: - v - workspace vector (used only in the case of split form)
175: Output Parameters:
176: + y - result vector
177: . A - Function matrix
178: - B - optional preconditioning matrix
180: Note:
181: If the nonlinear operator is represented in split form, the result
182: y = T(lambda)*x is computed without building T(lambda) explicitly. In
183: that case, parameters A and B are not used. Otherwise, the matrix
184: T(lambda) is built and the effect is the same as a call to
185: NEPComputeFunction() followed by a MatMult().
187: Level: developer
189: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
190: @*/
191: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)192: {
194: PetscInt i;
195: PetscScalar alpha;
206: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
207: VecSet(y,0.0);
208: for (i=0;i<nep->nt;i++) {
209: FNEvaluateFunction(nep->f[i],lambda,&alpha);
210: MatMult(nep->A[i],x,v);
211: VecAXPY(y,alpha,v);
212: }
213: } else {
214: NEPComputeFunction(nep,lambda,A,B);
215: MatMult(A,x,y);
216: }
217: return(0);
218: }
220: /*@
221: NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.
223: Collective on nep
225: Input Parameters:
226: + nep - the nonlinear eigensolver context
227: . lambda - scalar argument
228: . x - vector to be multiplied against
229: - v - workspace vector (used only in the case of split form)
231: Output Parameters:
232: + y - result vector
233: . A - Function matrix
234: - B - optional preconditioning matrix
236: Level: developer
238: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
239: @*/
240: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)241: {
243: PetscInt i;
244: PetscScalar alpha;
255: VecConjugate(x);
256: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
257: VecSet(y,0.0);
258: for (i=0;i<nep->nt;i++) {
259: FNEvaluateFunction(nep->f[i],lambda,&alpha);
260: MatMultTranspose(nep->A[i],x,v);
261: VecAXPY(y,alpha,v);
262: }
263: } else {
264: NEPComputeFunction(nep,lambda,A,B);
265: MatMultTranspose(A,x,y);
266: }
267: VecConjugate(x);
268: VecConjugate(y);
269: return(0);
270: }
272: /*@
273: NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
275: Collective on nep
277: Input Parameters:
278: + nep - the nonlinear eigensolver context
279: . lambda - scalar argument
280: . x - vector to be multiplied against
281: - v - workspace vector (used only in the case of split form)
283: Output Parameters:
284: + y - result vector
285: - A - Jacobian matrix
287: Note:
288: If the nonlinear operator is represented in split form, the result
289: y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
290: that case, parameter A is not used. Otherwise, the matrix
291: T'(lambda) is built and the effect is the same as a call to
292: NEPComputeJacobian() followed by a MatMult().
294: Level: developer
296: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
297: @*/
298: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)299: {
301: PetscInt i;
302: PetscScalar alpha;
312: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
313: VecSet(y,0.0);
314: for (i=0;i<nep->nt;i++) {
315: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
316: MatMult(nep->A[i],x,v);
317: VecAXPY(y,alpha,v);
318: }
319: } else {
320: NEPComputeJacobian(nep,lambda,A);
321: MatMult(A,x,y);
322: }
323: return(0);
324: }
326: /*@
327: NEPGetIterationNumber - Gets the current iteration number. If the
328: call to NEPSolve() is complete, then it returns the number of iterations
329: carried out by the solution method.
331: Not Collective
333: Input Parameter:
334: . nep - the nonlinear eigensolver context
336: Output Parameter:
337: . its - number of iterations
339: Note:
340: During the i-th iteration this call returns i-1. If NEPSolve() is
341: complete, then parameter "its" contains either the iteration number at
342: which convergence was successfully reached, or failure was detected.
343: Call NEPGetConvergedReason() to determine if the solver converged or
344: failed and why.
346: Level: intermediate
348: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
349: @*/
350: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)351: {
355: *its = nep->its;
356: return(0);
357: }
359: /*@
360: NEPGetConverged - Gets the number of converged eigenpairs.
362: Not Collective
364: Input Parameter:
365: . nep - the nonlinear eigensolver context
367: Output Parameter:
368: . nconv - number of converged eigenpairs
370: Note:
371: This function should be called after NEPSolve() has finished.
373: Level: beginner
375: .seealso: NEPSetDimensions(), NEPSolve()
376: @*/
377: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)378: {
382: NEPCheckSolved(nep,1);
383: *nconv = nep->nconv;
384: return(0);
385: }
387: /*@
388: NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
389: stopped.
391: Not Collective
393: Input Parameter:
394: . nep - the nonlinear eigensolver context
396: Output Parameter:
397: . reason - negative value indicates diverged, positive value converged
399: Notes:
401: Possible values for reason are
402: + NEP_CONVERGED_TOL - converged up to tolerance
403: . NEP_CONVERGED_USER - converged due to a user-defined condition
404: . NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
405: . NEP_DIVERGED_BREAKDOWN - generic breakdown in method
406: . NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
407: - NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
408: unrestarted solver
410: Can only be called after the call to NEPSolve() is complete.
412: Level: intermediate
414: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason415: @*/
416: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)417: {
421: NEPCheckSolved(nep,1);
422: *reason = nep->reason;
423: return(0);
424: }
426: /*@C
427: NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
428: NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
430: Logically Collective on nep
432: Input Parameters:
433: + nep - nonlinear eigensolver context
434: - i - index of the solution
436: Output Parameters:
437: + eigr - real part of eigenvalue
438: . eigi - imaginary part of eigenvalue
439: . Vr - real part of eigenvector
440: - Vi - imaginary part of eigenvector
442: Notes:
443: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
444: required. Otherwise, the caller must provide valid Vec objects, i.e.,
445: they must be created by the calling program with e.g. MatCreateVecs().
447: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
448: configured with complex scalars the eigenvalue is stored
449: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
450: set to zero). In any case, the user can pass NULL in Vr or Vi if one of
451: them is not required.
453: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
454: Eigenpairs are indexed according to the ordering criterion established
455: with NEPSetWhichEigenpairs().
457: Level: beginner
459: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
460: @*/
461: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)462: {
463: PetscInt k;
471: NEPCheckSolved(nep,1);
472: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
474: NEPComputeVectors(nep);
475: k = nep->perm[i];
477: /* eigenvalue */
478: #if defined(PETSC_USE_COMPLEX)
479: if (eigr) *eigr = nep->eigr[k];
480: if (eigi) *eigi = 0;
481: #else
482: if (eigr) *eigr = nep->eigr[k];
483: if (eigi) *eigi = nep->eigi[k];
484: #endif
486: /* eigenvector */
487: BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi);
488: return(0);
489: }
491: /*@
492: NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
494: Logically Collective on nep
496: Input Parameters:
497: + nep - eigensolver context
498: - i - index of the solution
500: Output Parameters:
501: + Wr - real part of left eigenvector
502: - Wi - imaginary part of left eigenvector
504: Notes:
505: The caller must provide valid Vec objects, i.e., they must be created
506: by the calling program with e.g. MatCreateVecs().
508: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
509: configured with complex scalars the eigenvector is stored directly in Wr
510: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
511: them is not required.
513: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
514: Eigensolutions are indexed according to the ordering criterion established
515: with NEPSetWhichEigenpairs().
517: Left eigenvectors are available only if the twosided flag was set, see
518: NEPSetTwoSided().
520: Level: intermediate
522: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
523: @*/
524: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)525: {
527: PetscInt k;
534: NEPCheckSolved(nep,1);
535: if (!nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
536: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
537: NEPComputeVectors(nep);
538: k = nep->perm[i];
539: BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi);
540: return(0);
541: }
543: /*@
544: NEPGetErrorEstimate - Returns the error estimate associated to the i-th
545: computed eigenpair.
547: Not Collective
549: Input Parameter:
550: + nep - nonlinear eigensolver context
551: - i - index of eigenpair
553: Output Parameter:
554: . errest - the error estimate
556: Notes:
557: This is the error estimate used internally by the eigensolver. The actual
558: error bound can be computed with NEPComputeRelativeError().
560: Level: advanced
562: .seealso: NEPComputeRelativeError()
563: @*/
564: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)565: {
569: NEPCheckSolved(nep,1);
570: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
571: *errest = nep->errest[nep->perm[i]];
572: return(0);
573: }
575: /*
576: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
577: associated with an eigenpair.
579: Input Parameters:
580: adj - whether the adjoint T^* must be used instead of T
581: lambda - eigenvalue
582: x - eigenvector
583: w - array of work vectors (two vectors in split form, one vector otherwise)
584: */
585: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)586: {
588: Vec y,z=NULL;
591: y = w[0];
592: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
593: if (adj) {
594: NEPApplyAdjoint(nep,lambda,x,z,y,nep->function,nep->function_pre);
595: } else {
596: NEPApplyFunction(nep,lambda,x,z,y,nep->function,nep->function_pre);
597: }
598: VecNorm(y,NORM_2,norm);
599: return(0);
600: }
602: /*@
603: NEPComputeError - Computes the error (based on the residual norm) associated
604: with the i-th computed eigenpair.
606: Collective on nep
608: Input Parameter:
609: + nep - the nonlinear eigensolver context
610: . i - the solution index
611: - type - the type of error to compute
613: Output Parameter:
614: . error - the error
616: Notes:
617: The error can be computed in various ways, all of them based on the residual
618: norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
619: eigenvector.
621: Level: beginner
623: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
624: @*/
625: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)626: {
628: Vec xr,xi=NULL;
629: PetscInt j,nwork,issplit=0;
630: PetscScalar kr,ki,s;
631: PetscReal er,z=0.0,errorl,nrm;
632: PetscBool flg;
639: NEPCheckSolved(nep,1);
641: /* allocate work vectors */
642: #if defined(PETSC_USE_COMPLEX)
643: nwork = 2;
644: #else
645: nwork = 3;
646: #endif
647: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
648: issplit = 1;
649: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
650: }
651: NEPSetWorkVecs(nep,nwork);
652: xr = nep->work[issplit+1];
653: #if !defined(PETSC_USE_COMPLEX)
654: xi = nep->work[issplit+2];
655: #endif
657: /* compute residual norms */
658: NEPGetEigenpair(nep,i,&kr,&ki,xr,xi);
659: #if !defined(PETSC_USE_COMPLEX)
660: if (ki) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented for complex eigenvalues with real scalars");
661: #endif
662: NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error);
663: VecNorm(xr,NORM_2,&er);
665: /* if two-sided, compute left residual norm and take the maximum */
666: if (nep->twosided) {
667: NEPGetLeftEigenvector(nep,i,xr,xi);
668: NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl);
669: *error = PetscMax(*error,errorl);
670: }
672: /* compute error */
673: switch (type) {
674: case NEP_ERROR_ABSOLUTE:
675: break;
676: case NEP_ERROR_RELATIVE:
677: *error /= PetscAbsScalar(kr)*er;
678: break;
679: case NEP_ERROR_BACKWARD:
680: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
681: NEPComputeFunction(nep,kr,nep->function,nep->function);
682: MatHasOperation(nep->function,MATOP_NORM,&flg);
683: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
684: MatNorm(nep->function,NORM_INFINITY,&nrm);
685: *error /= nrm*er;
686: break;
687: }
688: /* initialization of matrix norms */
689: if (!nep->nrma[0]) {
690: for (j=0;j<nep->nt;j++) {
691: MatHasOperation(nep->A[j],MATOP_NORM,&flg);
692: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
693: MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]);
694: }
695: }
696: for (j=0;j<nep->nt;j++) {
697: FNEvaluateFunction(nep->f[j],kr,&s);
698: z = z + nep->nrma[j]*PetscAbsScalar(s);
699: }
700: *error /= z*er;
701: break;
702: default:703: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
704: }
705: return(0);
706: }
708: /*@
709: NEPComputeFunction - Computes the function matrix T(lambda) that has been
710: set with NEPSetFunction().
712: Collective on nep
714: Input Parameters:
715: + nep - the NEP context
716: - lambda - the scalar argument
718: Output Parameters:
719: + A - Function matrix
720: - B - optional preconditioning matrix
722: Notes:
723: NEPComputeFunction() is typically used within nonlinear eigensolvers
724: implementations, so most users would not generally call this routine
725: themselves.
727: Level: developer
729: .seealso: NEPSetFunction(), NEPGetFunction()
730: @*/
731: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)732: {
734: PetscInt i;
735: PetscScalar alpha;
739: NEPCheckProblem(nep,1);
740: switch (nep->fui) {
741: case NEP_USER_INTERFACE_CALLBACK:
742: if (!nep->computefunction) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
743: PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0);
744: PetscStackPush("NEP user Function function");
745: (*nep->computefunction)(nep,lambda,A,B,nep->functionctx);
746: PetscStackPop;
747: PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0);
748: break;
749: case NEP_USER_INTERFACE_SPLIT:
750: MatZeroEntries(A);
751: for (i=0;i<nep->nt;i++) {
752: FNEvaluateFunction(nep->f[i],lambda,&alpha);
753: MatAXPY(A,alpha,nep->A[i],nep->mstr);
754: }
755: if (A != B) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented");
756: break;
757: }
758: return(0);
759: }
761: /*@
762: NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
763: set with NEPSetJacobian().
765: Collective on nep
767: Input Parameters:
768: + nep - the NEP context
769: - lambda - the scalar argument
771: Output Parameters:
772: . A - Jacobian matrix
774: Notes:
775: Most users should not need to explicitly call this routine, as it
776: is used internally within the nonlinear eigensolvers.
778: Level: developer
780: .seealso: NEPSetJacobian(), NEPGetJacobian()
781: @*/
782: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)783: {
785: PetscInt i;
786: PetscScalar alpha;
790: NEPCheckProblem(nep,1);
791: switch (nep->fui) {
792: case NEP_USER_INTERFACE_CALLBACK:
793: if (!nep->computejacobian) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
794: PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0);
795: PetscStackPush("NEP user Jacobian function");
796: (*nep->computejacobian)(nep,lambda,A,nep->jacobianctx);
797: PetscStackPop;
798: PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0);
799: break;
800: case NEP_USER_INTERFACE_SPLIT:
801: MatZeroEntries(A);
802: for (i=0;i<nep->nt;i++) {
803: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
804: MatAXPY(A,alpha,nep->A[i],nep->mstr);
805: }
806: break;
807: }
808: return(0);
809: }