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: (*nep->ops->solve)(nep);
86: nep->state = NEP_STATE_SOLVED;
88: if (!nep->reason) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
90: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
91: NEPComputeVectors(nep);
92: NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv);
93: nep->state = NEP_STATE_EIGENVECTORS;
94: }
96: /* sort eigenvalues according to nep->which parameter */
97: SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm);
98: PetscLogEventEnd(NEP_Solve,nep,0,0,0);
100: /* various viewers */
101: NEPViewFromOptions(nep,NULL,"-nep_view");
102: NEPReasonViewFromOptions(nep);
103: NEPErrorViewFromOptions(nep);
104: NEPValuesViewFromOptions(nep);
105: NEPVectorsViewFromOptions(nep);
107: /* Remove the initial subspace */
108: nep->nini = 0;
110: /* Reset resolvent information */
111: MatDestroy(&nep->resolvent);
112: return(0);
113: }
115: /*@
116: NEPProjectOperator - Computes the projection of the nonlinear operator.
118: Collective on nep
120: Input Parameters:
121: + nep - the nonlinear eigensolver context
122: . j0 - initial index
123: - j1 - final index
125: Notes:
126: This is available for split operator only.
128: The nonlinear operator T(lambda) is projected onto span(V), where V is
129: an orthonormal basis built internally by the solver. The projected
130: operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
131: computes all matrices Ei = V'*A_i*V, and stores them in the extra
132: matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
133: the previous ones are assumed to be available already.
135: Level: developer
137: .seealso: NEPSetSplitOperator()
138: @*/
139: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)140: {
142: PetscInt k;
143: Mat G;
149: NEPCheckProblem(nep,1);
150: NEPCheckSplit(nep,1);
151: BVSetActiveColumns(nep->V,j0,j1);
152: for (k=0;k<nep->nt;k++) {
153: DSGetMat(nep->ds,DSMatExtra[k],&G);
154: BVMatProject(nep->V,nep->A[k],nep->V,G);
155: DSRestoreMat(nep->ds,DSMatExtra[k],&G);
156: }
157: return(0);
158: }
160: /*@
161: NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
163: Collective on nep
165: Input Parameters:
166: + nep - the nonlinear eigensolver context
167: . lambda - scalar argument
168: . x - vector to be multiplied against
169: - v - workspace vector (used only in the case of split form)
171: Output Parameters:
172: + y - result vector
173: . A - Function matrix
174: - B - optional preconditioning matrix
176: Note:
177: If the nonlinear operator is represented in split form, the result
178: y = T(lambda)*x is computed without building T(lambda) explicitly. In
179: that case, parameters A and B are not used. Otherwise, the matrix
180: T(lambda) is built and the effect is the same as a call to
181: NEPComputeFunction() followed by a MatMult().
183: Level: developer
185: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
186: @*/
187: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)188: {
190: PetscInt i;
191: PetscScalar alpha;
202: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
203: VecSet(y,0.0);
204: for (i=0;i<nep->nt;i++) {
205: FNEvaluateFunction(nep->f[i],lambda,&alpha);
206: MatMult(nep->A[i],x,v);
207: VecAXPY(y,alpha,v);
208: }
209: } else {
210: NEPComputeFunction(nep,lambda,A,B);
211: MatMult(A,x,y);
212: }
213: return(0);
214: }
216: /*@
217: NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.
219: Collective on nep
221: Input Parameters:
222: + nep - the nonlinear eigensolver context
223: . lambda - scalar argument
224: . x - vector to be multiplied against
225: - v - workspace vector (used only in the case of split form)
227: Output Parameters:
228: + y - result vector
229: . A - Function matrix
230: - B - optional preconditioning matrix
232: Level: developer
234: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
235: @*/
236: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)237: {
239: PetscInt i;
240: PetscScalar alpha;
251: VecConjugate(x);
252: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
253: VecSet(y,0.0);
254: for (i=0;i<nep->nt;i++) {
255: FNEvaluateFunction(nep->f[i],lambda,&alpha);
256: MatMultTranspose(nep->A[i],x,v);
257: VecAXPY(y,alpha,v);
258: }
259: } else {
260: NEPComputeFunction(nep,lambda,A,B);
261: MatMultTranspose(A,x,y);
262: }
263: VecConjugate(x);
264: VecConjugate(y);
265: return(0);
266: }
268: /*@
269: NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
271: Collective on nep
273: Input Parameters:
274: + nep - the nonlinear eigensolver context
275: . lambda - scalar argument
276: . x - vector to be multiplied against
277: - v - workspace vector (used only in the case of split form)
279: Output Parameters:
280: + y - result vector
281: - A - Jacobian matrix
283: Note:
284: If the nonlinear operator is represented in split form, the result
285: y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
286: that case, parameter A is not used. Otherwise, the matrix
287: T'(lambda) is built and the effect is the same as a call to
288: NEPComputeJacobian() followed by a MatMult().
290: Level: developer
292: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
293: @*/
294: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)295: {
297: PetscInt i;
298: PetscScalar alpha;
308: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
309: VecSet(y,0.0);
310: for (i=0;i<nep->nt;i++) {
311: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
312: MatMult(nep->A[i],x,v);
313: VecAXPY(y,alpha,v);
314: }
315: } else {
316: NEPComputeJacobian(nep,lambda,A);
317: MatMult(A,x,y);
318: }
319: return(0);
320: }
322: /*@
323: NEPGetIterationNumber - Gets the current iteration number. If the
324: call to NEPSolve() is complete, then it returns the number of iterations
325: carried out by the solution method.
327: Not Collective
329: Input Parameter:
330: . nep - the nonlinear eigensolver context
332: Output Parameter:
333: . its - number of iterations
335: Note:
336: During the i-th iteration this call returns i-1. If NEPSolve() is
337: complete, then parameter "its" contains either the iteration number at
338: which convergence was successfully reached, or failure was detected.
339: Call NEPGetConvergedReason() to determine if the solver converged or
340: failed and why.
342: Level: intermediate
344: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
345: @*/
346: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)347: {
351: *its = nep->its;
352: return(0);
353: }
355: /*@
356: NEPGetConverged - Gets the number of converged eigenpairs.
358: Not Collective
360: Input Parameter:
361: . nep - the nonlinear eigensolver context
363: Output Parameter:
364: . nconv - number of converged eigenpairs
366: Note:
367: This function should be called after NEPSolve() has finished.
369: Level: beginner
371: .seealso: NEPSetDimensions(), NEPSolve()
372: @*/
373: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)374: {
378: NEPCheckSolved(nep,1);
379: *nconv = nep->nconv;
380: return(0);
381: }
383: /*@
384: NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
385: stopped.
387: Not Collective
389: Input Parameter:
390: . nep - the nonlinear eigensolver context
392: Output Parameter:
393: . reason - negative value indicates diverged, positive value converged
395: Notes:
397: Possible values for reason are
398: + NEP_CONVERGED_TOL - converged up to tolerance
399: . NEP_CONVERGED_USER - converged due to a user-defined condition
400: . NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
401: . NEP_DIVERGED_BREAKDOWN - generic breakdown in method
402: . NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
403: - NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
404: unrestarted solver
406: Can only be called after the call to NEPSolve() is complete.
408: Level: intermediate
410: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason411: @*/
412: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)413: {
417: NEPCheckSolved(nep,1);
418: *reason = nep->reason;
419: return(0);
420: }
422: /*@C
423: NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
424: NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
426: Logically Collective on nep
428: Input Parameters:
429: + nep - nonlinear eigensolver context
430: - i - index of the solution
432: Output Parameters:
433: + eigr - real part of eigenvalue
434: . eigi - imaginary part of eigenvalue
435: . Vr - real part of eigenvector
436: - Vi - imaginary part of eigenvector
438: Notes:
439: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
440: required. Otherwise, the caller must provide valid Vec objects, i.e.,
441: they must be created by the calling program with e.g. MatCreateVecs().
443: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
444: configured with complex scalars the eigenvalue is stored
445: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
446: set to zero). In any case, the user can pass NULL in Vr or Vi if one of
447: them is not required.
449: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
450: Eigenpairs are indexed according to the ordering criterion established
451: with NEPSetWhichEigenpairs().
453: Level: beginner
455: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
456: @*/
457: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)458: {
459: PetscInt k;
467: NEPCheckSolved(nep,1);
468: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
470: NEPComputeVectors(nep);
471: k = nep->perm[i];
473: /* eigenvalue */
474: #if defined(PETSC_USE_COMPLEX)
475: if (eigr) *eigr = nep->eigr[k];
476: if (eigi) *eigi = 0;
477: #else
478: if (eigr) *eigr = nep->eigr[k];
479: if (eigi) *eigi = nep->eigi[k];
480: #endif
482: /* eigenvector */
483: BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi);
484: return(0);
485: }
487: /*@
488: NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
490: Logically Collective on nep
492: Input Parameters:
493: + nep - eigensolver context
494: - i - index of the solution
496: Output Parameters:
497: + Wr - real part of left eigenvector
498: - Wi - imaginary part of left eigenvector
500: Notes:
501: The caller must provide valid Vec objects, i.e., they must be created
502: by the calling program with e.g. MatCreateVecs().
504: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
505: configured with complex scalars the eigenvector is stored directly in Wr
506: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
507: them is not required.
509: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
510: Eigensolutions are indexed according to the ordering criterion established
511: with NEPSetWhichEigenpairs().
513: Left eigenvectors are available only if the twosided flag was set, see
514: NEPSetTwoSided().
516: Level: intermediate
518: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
519: @*/
520: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)521: {
523: PetscInt k;
530: NEPCheckSolved(nep,1);
531: if (!nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
532: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
533: NEPComputeVectors(nep);
534: k = nep->perm[i];
535: BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi);
536: return(0);
537: }
539: /*@
540: NEPGetErrorEstimate - Returns the error estimate associated to the i-th
541: computed eigenpair.
543: Not Collective
545: Input Parameter:
546: + nep - nonlinear eigensolver context
547: - i - index of eigenpair
549: Output Parameter:
550: . errest - the error estimate
552: Notes:
553: This is the error estimate used internally by the eigensolver. The actual
554: error bound can be computed with NEPComputeRelativeError().
556: Level: advanced
558: .seealso: NEPComputeRelativeError()
559: @*/
560: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)561: {
565: NEPCheckSolved(nep,1);
566: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
567: *errest = nep->errest[nep->perm[i]];
568: return(0);
569: }
571: /*
572: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
573: associated with an eigenpair.
575: Input Parameters:
576: adj - whether the adjoint T^* must be used instead of T
577: lambda - eigenvalue
578: x - eigenvector
579: w - array of work vectors (two vectors in split form, one vector otherwise)
580: */
581: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)582: {
584: Vec y,z=NULL;
587: y = w[0];
588: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
589: if (adj) {
590: NEPApplyAdjoint(nep,lambda,x,z,y,nep->function,nep->function_pre);
591: } else {
592: NEPApplyFunction(nep,lambda,x,z,y,nep->function,nep->function_pre);
593: }
594: VecNorm(y,NORM_2,norm);
595: return(0);
596: }
598: /*@
599: NEPComputeError - Computes the error (based on the residual norm) associated
600: with the i-th computed eigenpair.
602: Collective on nep
604: Input Parameter:
605: + nep - the nonlinear eigensolver context
606: . i - the solution index
607: - type - the type of error to compute
609: Output Parameter:
610: . error - the error
612: Notes:
613: The error can be computed in various ways, all of them based on the residual
614: norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
615: eigenvector.
617: Level: beginner
619: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
620: @*/
621: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)622: {
624: Vec xr,xi=NULL;
625: PetscInt j,nwork,issplit=0;
626: PetscScalar kr,ki,s;
627: PetscReal er,z=0.0,errorl,nrm;
628: PetscBool flg;
635: NEPCheckSolved(nep,1);
637: /* allocate work vectors */
638: #if defined(PETSC_USE_COMPLEX)
639: nwork = 2;
640: #else
641: nwork = 3;
642: #endif
643: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
644: issplit = 1;
645: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
646: }
647: NEPSetWorkVecs(nep,nwork);
648: xr = nep->work[issplit+1];
649: #if !defined(PETSC_USE_COMPLEX)
650: xi = nep->work[issplit+2];
651: #endif
653: /* compute residual norms */
654: NEPGetEigenpair(nep,i,&kr,&ki,xr,xi);
655: #if !defined(PETSC_USE_COMPLEX)
656: if (ki) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented for complex eigenvalues with real scalars");
657: #endif
658: NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error);
659: VecNorm(xr,NORM_2,&er);
661: /* if two-sided, compute left residual norm and take the maximum */
662: if (nep->twosided) {
663: NEPGetLeftEigenvector(nep,i,xr,xi);
664: NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl);
665: *error = PetscMax(*error,errorl);
666: }
668: /* compute error */
669: switch (type) {
670: case NEP_ERROR_ABSOLUTE:
671: break;
672: case NEP_ERROR_RELATIVE:
673: *error /= PetscAbsScalar(kr)*er;
674: break;
675: case NEP_ERROR_BACKWARD:
676: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
677: NEPComputeFunction(nep,kr,nep->function,nep->function);
678: MatHasOperation(nep->function,MATOP_NORM,&flg);
679: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
680: MatNorm(nep->function,NORM_INFINITY,&nrm);
681: *error /= nrm*er;
682: break;
683: }
684: /* initialization of matrix norms */
685: if (!nep->nrma[0]) {
686: for (j=0;j<nep->nt;j++) {
687: MatHasOperation(nep->A[j],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->A[j],NORM_INFINITY,&nep->nrma[j]);
690: }
691: }
692: for (j=0;j<nep->nt;j++) {
693: FNEvaluateFunction(nep->f[j],kr,&s);
694: z = z + nep->nrma[j]*PetscAbsScalar(s);
695: }
696: *error /= z*er;
697: break;
698: default:699: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
700: }
701: return(0);
702: }
704: /*@
705: NEPComputeFunction - Computes the function matrix T(lambda) that has been
706: set with NEPSetFunction().
708: Collective on nep
710: Input Parameters:
711: + nep - the NEP context
712: - lambda - the scalar argument
714: Output Parameters:
715: + A - Function matrix
716: - B - optional preconditioning matrix
718: Notes:
719: NEPComputeFunction() is typically used within nonlinear eigensolvers
720: implementations, so most users would not generally call this routine
721: themselves.
723: Level: developer
725: .seealso: NEPSetFunction(), NEPGetFunction()
726: @*/
727: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)728: {
730: PetscInt i;
731: PetscScalar alpha;
735: NEPCheckProblem(nep,1);
736: switch (nep->fui) {
737: case NEP_USER_INTERFACE_CALLBACK:
738: if (!nep->computefunction) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
739: PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0);
740: PetscStackPush("NEP user Function function");
741: (*nep->computefunction)(nep,lambda,A,B,nep->functionctx);
742: PetscStackPop;
743: PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0);
744: break;
745: case NEP_USER_INTERFACE_SPLIT:
746: MatZeroEntries(A);
747: for (i=0;i<nep->nt;i++) {
748: FNEvaluateFunction(nep->f[i],lambda,&alpha);
749: MatAXPY(A,alpha,nep->A[i],nep->mstr);
750: }
751: if (A != B) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented");
752: break;
753: }
754: return(0);
755: }
757: /*@
758: NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
759: set with NEPSetJacobian().
761: Collective on nep
763: Input Parameters:
764: + nep - the NEP context
765: - lambda - the scalar argument
767: Output Parameters:
768: . A - Jacobian matrix
770: Notes:
771: Most users should not need to explicitly call this routine, as it
772: is used internally within the nonlinear eigensolvers.
774: Level: developer
776: .seealso: NEPSetJacobian(), NEPGetJacobian()
777: @*/
778: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)779: {
781: PetscInt i;
782: PetscScalar alpha;
786: NEPCheckProblem(nep,1);
787: switch (nep->fui) {
788: case NEP_USER_INTERFACE_CALLBACK:
789: if (!nep->computejacobian) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
790: PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0);
791: PetscStackPush("NEP user Jacobian function");
792: (*nep->computejacobian)(nep,lambda,A,nep->jacobianctx);
793: PetscStackPop;
794: PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0);
795: break;
796: case NEP_USER_INTERFACE_SPLIT:
797: MatZeroEntries(A);
798: for (i=0;i<nep->nt;i++) {
799: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
800: MatAXPY(A,alpha,nep->A[i],nep->mstr);
801: }
802: break;
803: }
804: return(0);
805: }