1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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> /*I "slepcnep.h" I*/
15: #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
16: #include <petscdraw.h>
18: PetscErrorCode NEPComputeVectors(NEP nep) 19: {
23: NEPCheckSolved(nep,1);
24: if (nep->state==NEP_STATE_SOLVED && nep->ops->computevectors) {
25: (*nep->ops->computevectors)(nep);
26: }
27: nep->state = NEP_STATE_EIGENVECTORS;
28: return(0);
29: }
31: /*@
32: NEPSolve - Solves the nonlinear eigensystem.
34: Collective on NEP 36: Input Parameter:
37: . nep - eigensolver context obtained from NEPCreate()
39: Options Database Keys:
40: + -nep_view - print information about the solver used
41: . -nep_view_vectors binary - save the computed eigenvectors to the default binary viewer
42: . -nep_view_values - print computed eigenvalues
43: . -nep_converged_reason - print reason for convergence, and number of iterations
44: . -nep_error_absolute - print absolute errors of each eigenpair
45: - -nep_error_relative - print relative errors of each eigenpair
47: Level: beginner
49: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
50: @*/
51: PetscErrorCode NEPSolve(NEP nep) 52: {
54: PetscInt i;
58: if (nep->state>=NEP_STATE_SOLVED) return(0);
59: PetscLogEventBegin(NEP_Solve,nep,0,0,0);
61: /* call setup */
62: NEPSetUp(nep);
63: nep->nconv = 0;
64: nep->its = 0;
65: for (i=0;i<nep->ncv;i++) {
66: nep->eigr[i] = 0.0;
67: nep->eigi[i] = 0.0;
68: nep->errest[i] = 0.0;
69: nep->perm[i] = i;
70: }
71: NEPViewFromOptions(nep,NULL,"-nep_view_pre");
72: RGViewFromOptions(nep->rg,NULL,"-rg_view");
74: (*nep->ops->solve)(nep);
75: nep->state = NEP_STATE_SOLVED;
77: if (!nep->reason) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
79: if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
80: NEPComputeVectors(nep);
81: NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv);
82: nep->state = NEP_STATE_EIGENVECTORS;
83: }
85: /* sort eigenvalues according to nep->which parameter */
86: SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm);
87: PetscLogEventEnd(NEP_Solve,nep,0,0,0);
89: /* various viewers */
90: NEPViewFromOptions(nep,NULL,"-nep_view");
91: NEPReasonViewFromOptions(nep);
92: NEPErrorViewFromOptions(nep);
93: NEPValuesViewFromOptions(nep);
94: NEPVectorsViewFromOptions(nep);
96: /* Remove the initial subspace */
97: nep->nini = 0;
99: /* Reset resolvent information */
100: MatDestroy(&nep->resolvent);
101: return(0);
102: }
104: /*@
105: NEPProjectOperator - Computes the projection of the nonlinear operator.
107: Collective on NEP109: Input Parameters:
110: + nep - the nonlinear eigensolver context
111: . j0 - initial index
112: - j1 - final index
114: Notes:
115: This is available for split operator only.
117: The nonlinear operator T(lambda) is projected onto span(V), where V is
118: an orthonormal basis built internally by the solver. The projected
119: operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
120: computes all matrices Ei = V'*A_i*V, and stores them in the extra
121: matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
122: the previous ones are assumed to be available already.
124: Level: developer
126: .seealso: NEPSetSplitOperator()
127: @*/
128: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)129: {
131: PetscInt k;
132: Mat G;
138: NEPCheckProblem(nep,1);
139: NEPCheckSplit(nep,1);
140: BVSetActiveColumns(nep->V,j0,j1);
141: for (k=0;k<nep->nt;k++) {
142: DSGetMat(nep->ds,DSMatExtra[k],&G);
143: BVMatProject(nep->V,nep->A[k],nep->V,G);
144: DSRestoreMat(nep->ds,DSMatExtra[k],&G);
145: }
146: return(0);
147: }
149: /*@
150: NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.
152: Collective on NEP154: Input Parameters:
155: + nep - the nonlinear eigensolver context
156: . lambda - scalar argument
157: . x - vector to be multiplied against
158: - v - workspace vector (used only in the case of split form)
160: Output Parameters:
161: + y - result vector
162: . A - Function matrix
163: - B - optional preconditioning matrix
165: Note:
166: If the nonlinear operator is represented in split form, the result
167: y = T(lambda)*x is computed without building T(lambda) explicitly. In
168: that case, parameters A and B are not used. Otherwise, the matrix
169: T(lambda) is built and the effect is the same as a call to
170: NEPComputeFunction() followed by a MatMult().
172: Level: developer
174: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
175: @*/
176: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)177: {
179: PetscInt i;
180: PetscScalar alpha;
191: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
192: VecSet(y,0.0);
193: for (i=0;i<nep->nt;i++) {
194: FNEvaluateFunction(nep->f[i],lambda,&alpha);
195: MatMult(nep->A[i],x,v);
196: VecAXPY(y,alpha,v);
197: }
198: } else {
199: NEPComputeFunction(nep,lambda,A,B);
200: MatMult(A,x,y);
201: }
202: return(0);
203: }
205: /*@
206: NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.
208: Collective on NEP210: Input Parameters:
211: + nep - the nonlinear eigensolver context
212: . lambda - scalar argument
213: . x - vector to be multiplied against
214: - v - workspace vector (used only in the case of split form)
216: Output Parameters:
217: + y - result vector
218: . A - Function matrix
219: - B - optional preconditioning matrix
221: Level: developer
223: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
224: @*/
225: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)226: {
228: PetscInt i;
229: PetscScalar alpha;
240: VecConjugate(x);
241: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
242: VecSet(y,0.0);
243: for (i=0;i<nep->nt;i++) {
244: FNEvaluateFunction(nep->f[i],lambda,&alpha);
245: MatMultTranspose(nep->A[i],x,v);
246: VecAXPY(y,alpha,v);
247: }
248: } else {
249: NEPComputeFunction(nep,lambda,A,B);
250: MatMultTranspose(A,x,y);
251: }
252: VecConjugate(x);
253: VecConjugate(y);
254: return(0);
255: }
257: /*@
258: NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.
260: Collective on NEP262: Input Parameters:
263: + nep - the nonlinear eigensolver context
264: . lambda - scalar argument
265: . x - vector to be multiplied against
266: - v - workspace vector (used only in the case of split form)
268: Output Parameters:
269: + y - result vector
270: - A - Jacobian matrix
272: Note:
273: If the nonlinear operator is represented in split form, the result
274: y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
275: that case, parameter A is not used. Otherwise, the matrix
276: T'(lambda) is built and the effect is the same as a call to
277: NEPComputeJacobian() followed by a MatMult().
279: Level: developer
281: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
282: @*/
283: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)284: {
286: PetscInt i;
287: PetscScalar alpha;
297: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
298: VecSet(y,0.0);
299: for (i=0;i<nep->nt;i++) {
300: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
301: MatMult(nep->A[i],x,v);
302: VecAXPY(y,alpha,v);
303: }
304: } else {
305: NEPComputeJacobian(nep,lambda,A);
306: MatMult(A,x,y);
307: }
308: return(0);
309: }
311: /*@
312: NEPGetIterationNumber - Gets the current iteration number. If the
313: call to NEPSolve() is complete, then it returns the number of iterations
314: carried out by the solution method.
316: Not Collective
318: Input Parameter:
319: . nep - the nonlinear eigensolver context
321: Output Parameter:
322: . its - number of iterations
324: Note:
325: During the i-th iteration this call returns i-1. If NEPSolve() is
326: complete, then parameter "its" contains either the iteration number at
327: which convergence was successfully reached, or failure was detected.
328: Call NEPGetConvergedReason() to determine if the solver converged or
329: failed and why.
331: Level: intermediate
333: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
334: @*/
335: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)336: {
340: *its = nep->its;
341: return(0);
342: }
344: /*@
345: NEPGetConverged - Gets the number of converged eigenpairs.
347: Not Collective
349: Input Parameter:
350: . nep - the nonlinear eigensolver context
352: Output Parameter:
353: . nconv - number of converged eigenpairs
355: Note:
356: This function should be called after NEPSolve() has finished.
358: Level: beginner
360: .seealso: NEPSetDimensions(), NEPSolve()
361: @*/
362: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)363: {
367: NEPCheckSolved(nep,1);
368: *nconv = nep->nconv;
369: return(0);
370: }
372: /*@
373: NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
374: stopped.
376: Not Collective
378: Input Parameter:
379: . nep - the nonlinear eigensolver context
381: Output Parameter:
382: . reason - negative value indicates diverged, positive value converged
384: Notes:
386: Possible values for reason are
387: + NEP_CONVERGED_TOL - converged up to tolerance
388: . NEP_CONVERGED_USER - converged due to a user-defined condition
389: . NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
390: . NEP_DIVERGED_BREAKDOWN - generic breakdown in method
391: . NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
392: - NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
393: unrestarted solver
395: Can only be called after the call to NEPSolve() is complete.
397: Level: intermediate
399: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason400: @*/
401: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)402: {
406: NEPCheckSolved(nep,1);
407: *reason = nep->reason;
408: return(0);
409: }
411: /*@C
412: NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
413: NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.
415: Logically Collective on NEP417: Input Parameters:
418: + nep - nonlinear eigensolver context
419: - i - index of the solution
421: Output Parameters:
422: + eigr - real part of eigenvalue
423: . eigi - imaginary part of eigenvalue
424: . Vr - real part of eigenvector
425: - Vi - imaginary part of eigenvector
427: Notes:
428: It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
429: required. Otherwise, the caller must provide valid Vec objects, i.e.,
430: they must be created by the calling program with e.g. MatCreateVecs().
432: If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
433: configured with complex scalars the eigenvalue is stored
434: directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
435: set to zero). In any case, the user can pass NULL in Vr or Vi if one of
436: them is not required.
438: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
439: Eigenpairs are indexed according to the ordering criterion established
440: with NEPSetWhichEigenpairs().
442: Level: beginner
444: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
445: @*/
446: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)447: {
448: PetscInt k;
456: NEPCheckSolved(nep,1);
457: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
459: NEPComputeVectors(nep);
460: k = nep->perm[i];
462: /* eigenvalue */
463: #if defined(PETSC_USE_COMPLEX)
464: if (eigr) *eigr = nep->eigr[k];
465: if (eigi) *eigi = 0;
466: #else
467: if (eigr) *eigr = nep->eigr[k];
468: if (eigi) *eigi = nep->eigi[k];
469: #endif
471: /* eigenvector */
472: BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi);
473: return(0);
474: }
476: /*@
477: NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().
479: Logically Collective on NEP481: Input Parameters:
482: + nep - eigensolver context
483: - i - index of the solution
485: Output Parameters:
486: + Wr - real part of left eigenvector
487: - Wi - imaginary part of left eigenvector
489: Notes:
490: The caller must provide valid Vec objects, i.e., they must be created
491: by the calling program with e.g. MatCreateVecs().
493: If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
494: configured with complex scalars the eigenvector is stored directly in Wr
495: (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
496: them is not required.
498: The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
499: Eigensolutions are indexed according to the ordering criterion established
500: with NEPSetWhichEigenpairs().
502: Left eigenvectors are available only if the twosided flag was set, see
503: NEPSetTwoSided().
505: Level: intermediate
507: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
508: @*/
509: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)510: {
512: PetscInt k;
519: NEPCheckSolved(nep,1);
520: if (!nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
521: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
522: NEPComputeVectors(nep);
523: k = nep->perm[i];
524: BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi);
525: return(0);
526: }
528: /*@
529: NEPGetErrorEstimate - Returns the error estimate associated to the i-th
530: computed eigenpair.
532: Not Collective
534: Input Parameter:
535: + nep - nonlinear eigensolver context
536: - i - index of eigenpair
538: Output Parameter:
539: . errest - the error estimate
541: Notes:
542: This is the error estimate used internally by the eigensolver. The actual
543: error bound can be computed with NEPComputeRelativeError().
545: Level: advanced
547: .seealso: NEPComputeRelativeError()
548: @*/
549: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)550: {
554: NEPCheckSolved(nep,1);
555: if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
556: *errest = nep->errest[nep->perm[i]];
557: return(0);
558: }
560: /*
561: NEPComputeResidualNorm_Private - Computes the norm of the residual vector
562: associated with an eigenpair.
564: Input Parameters:
565: adj - whether the adjoint T^* must be used instead of T
566: lambda - eigenvalue
567: x - eigenvector
568: w - array of work vectors (two vectors in split form, one vector otherwise)
569: */
570: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)571: {
573: Vec y,z=NULL;
576: y = w[0];
577: if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
578: if (adj) {
579: NEPApplyAdjoint(nep,lambda,x,z,y,nep->function,nep->function_pre);
580: } else {
581: NEPApplyFunction(nep,lambda,x,z,y,nep->function,nep->function_pre);
582: }
583: VecNorm(y,NORM_2,norm);
584: return(0);
585: }
587: /*@
588: NEPComputeError - Computes the error (based on the residual norm) associated
589: with the i-th computed eigenpair.
591: Collective on NEP593: Input Parameter:
594: + nep - the nonlinear eigensolver context
595: . i - the solution index
596: - type - the type of error to compute
598: Output Parameter:
599: . error - the error
601: Notes:
602: The error can be computed in various ways, all of them based on the residual
603: norm computed as ||T(lambda)x||_2 where lambda is the eigenvalue and x is the
604: eigenvector.
606: Level: beginner
608: .seealso: NEPErrorType, NEPSolve(), NEPGetErrorEstimate()
609: @*/
610: PetscErrorCode NEPComputeError(NEP nep,PetscInt i,NEPErrorType type,PetscReal *error)611: {
613: Vec xr,xi=NULL;
614: PetscInt j,nwork,issplit=0;
615: PetscScalar kr,ki,s;
616: PetscReal er,z=0.0,errorl;
617: PetscBool flg;
624: NEPCheckSolved(nep,1);
626: /* allocate work vectors */
627: #if defined(PETSC_USE_COMPLEX)
628: nwork = 2;
629: #else
630: nwork = 3;
631: #endif
632: if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
633: issplit = 1;
634: nwork++; /* need an extra work vector for NEPComputeResidualNorm_Private */
635: }
636: NEPSetWorkVecs(nep,nwork);
637: xr = nep->work[issplit+1];
638: #if !defined(PETSC_USE_COMPLEX)
639: xi = nep->work[issplit+2];
640: #endif
642: /* compute residual norms */
643: NEPGetEigenpair(nep,i,&kr,&ki,xr,xi);
644: #if !defined(PETSC_USE_COMPLEX)
645: if (ki) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented for complex eigenvalues with real scalars");
646: #endif
647: NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error);
648: VecNorm(xr,NORM_2,&er);
650: /* if two-sided, compute left residual norm and take the maximum */
651: if (nep->twosided) {
652: NEPGetLeftEigenvector(nep,i,xr,xi);
653: NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl);
654: *error = PetscMax(*error,errorl);
655: }
657: /* compute error */
658: switch (type) {
659: case NEP_ERROR_ABSOLUTE:
660: break;
661: case NEP_ERROR_RELATIVE:
662: *error /= PetscAbsScalar(kr)*er;
663: break;
664: case NEP_ERROR_BACKWARD:
665: if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
666: *error = 0.0;
667: PetscInfo(nep,"Backward error only available in split form\n");
668: break;
669: }
670: /* initialization of matrix norms */
671: if (!nep->nrma[0]) {
672: for (j=0;j<nep->nt;j++) {
673: MatHasOperation(nep->A[j],MATOP_NORM,&flg);
674: if (!flg) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"The computation of backward errors requires a matrix norm operation");
675: MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]);
676: }
677: }
678: for (j=0;j<nep->nt;j++) {
679: FNEvaluateFunction(nep->f[j],kr,&s);
680: z = z + nep->nrma[j]*PetscAbsScalar(s);
681: }
682: *error /= z;
683: break;
684: default:685: SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
686: }
687: return(0);
688: }
690: /*@
691: NEPComputeFunction - Computes the function matrix T(lambda) that has been
692: set with NEPSetFunction().
694: Collective on NEP and Mat
696: Input Parameters:
697: + nep - the NEP context
698: - lambda - the scalar argument
700: Output Parameters:
701: + A - Function matrix
702: - B - optional preconditioning matrix
704: Notes:
705: NEPComputeFunction() is typically used within nonlinear eigensolvers
706: implementations, so most users would not generally call this routine
707: themselves.
709: Level: developer
711: .seealso: NEPSetFunction(), NEPGetFunction()
712: @*/
713: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)714: {
716: PetscInt i;
717: PetscScalar alpha;
721: NEPCheckProblem(nep,1);
722: switch (nep->fui) {
723: case NEP_USER_INTERFACE_CALLBACK:
724: if (!nep->computefunction) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
725: PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0);
726: PetscStackPush("NEP user Function function");
727: (*nep->computefunction)(nep,lambda,A,B,nep->functionctx);
728: PetscStackPop;
729: PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0);
730: break;
731: case NEP_USER_INTERFACE_SPLIT:
732: MatZeroEntries(A);
733: for (i=0;i<nep->nt;i++) {
734: FNEvaluateFunction(nep->f[i],lambda,&alpha);
735: MatAXPY(A,alpha,nep->A[i],nep->mstr);
736: }
737: if (A != B) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented");
738: break;
739: case NEP_USER_INTERFACE_DERIVATIVES:
740: PetscLogEventBegin(NEP_DerivativesEval,nep,A,B,0);
741: PetscStackPush("NEP user Derivatives function");
742: (*nep->computederivatives)(nep,lambda,0,A,nep->derivativesctx);
743: PetscStackPop;
744: PetscLogEventEnd(NEP_DerivativesEval,nep,A,B,0);
745: break;
746: }
747: return(0);
748: }
750: /*@
751: NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
752: set with NEPSetJacobian().
754: Collective on NEP and Mat
756: Input Parameters:
757: + nep - the NEP context
758: - lambda - the scalar argument
760: Output Parameters:
761: . A - Jacobian matrix
763: Notes:
764: Most users should not need to explicitly call this routine, as it
765: is used internally within the nonlinear eigensolvers.
767: Level: developer
769: .seealso: NEPSetJacobian(), NEPGetJacobian()
770: @*/
771: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)772: {
774: PetscInt i;
775: PetscScalar alpha;
779: NEPCheckProblem(nep,1);
780: switch (nep->fui) {
781: case NEP_USER_INTERFACE_CALLBACK:
782: if (!nep->computejacobian) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
783: PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0);
784: PetscStackPush("NEP user Jacobian function");
785: (*nep->computejacobian)(nep,lambda,A,nep->jacobianctx);
786: PetscStackPop;
787: PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0);
788: break;
789: case NEP_USER_INTERFACE_SPLIT:
790: MatZeroEntries(A);
791: for (i=0;i<nep->nt;i++) {
792: FNEvaluateDerivative(nep->f[i],lambda,&alpha);
793: MatAXPY(A,alpha,nep->A[i],nep->mstr);
794: }
795: break;
796: case NEP_USER_INTERFACE_DERIVATIVES:
797: PetscLogEventBegin(NEP_DerivativesEval,nep,A,0,0);
798: PetscStackPush("NEP user Derivatives function");
799: (*nep->computederivatives)(nep,lambda,1,A,nep->derivativesctx);
800: PetscStackPop;
801: PetscLogEventEnd(NEP_DerivativesEval,nep,A,0,0);
802: break;
803: }
804: return(0);
805: }