Actual source code: nepsolve.c

slepc-3.17.1 2022-04-11
Report Typos and Errors
  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:            (to appear), arXiv:1910.11712, 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:   NEPCheckSolved(nep,1);
 40:   if (nep->state==NEP_STATE_SOLVED && nep->ops->computevectors) (*nep->ops->computevectors)(nep);
 41:   nep->state = NEP_STATE_EIGENVECTORS;
 42:   PetscFunctionReturn(0);
 43: }

 45: /*@
 46:    NEPSolve - Solves the nonlinear eigensystem.

 48:    Collective on nep

 50:    Input Parameter:
 51: .  nep - eigensolver context obtained from NEPCreate()

 53:    Options Database Keys:
 54: +  -nep_view - print information about the solver used
 55: .  -nep_view_vectors - view the computed eigenvectors
 56: .  -nep_view_values - view the computed eigenvalues
 57: .  -nep_converged_reason - print reason for convergence, and number of iterations
 58: .  -nep_error_absolute - print absolute errors of each eigenpair
 59: .  -nep_error_relative - print relative errors of each eigenpair
 60: -  -nep_error_backward - print backward errors of each eigenpair

 62:    Notes:
 63:    All the command-line options listed above admit an optional argument specifying
 64:    the viewer type and options. For instance, use '-nep_view_vectors binary:myvecs.bin'
 65:    to save the eigenvectors to a binary file, '-nep_view_values draw' to draw the computed
 66:    eigenvalues graphically, or '-nep_error_relative :myerr.m:ascii_matlab' to save
 67:    the errors in a file that can be executed in Matlab.

 69:    Level: beginner

 71: .seealso: NEPCreate(), NEPSetUp(), NEPDestroy(), NEPSetTolerances()
 72: @*/
 73: PetscErrorCode NEPSolve(NEP nep)
 74: {
 75:   PetscInt       i;

 78:   if (nep->state>=NEP_STATE_SOLVED) PetscFunctionReturn(0);
 79:   PetscCitationsRegister(citation,&cited);
 80:   PetscLogEventBegin(NEP_Solve,nep,0,0,0);

 82:   /* call setup */
 83:   NEPSetUp(nep);
 84:   nep->nconv = 0;
 85:   nep->its = 0;
 86:   for (i=0;i<nep->ncv;i++) {
 87:     nep->eigr[i]   = 0.0;
 88:     nep->eigi[i]   = 0.0;
 89:     nep->errest[i] = 0.0;
 90:     nep->perm[i]   = i;
 91:   }
 92:   NEPViewFromOptions(nep,NULL,"-nep_view_pre");
 93:   RGViewFromOptions(nep->rg,NULL,"-rg_view");

 95:   /* call solver */
 96:   (*nep->ops->solve)(nep);
 98:   nep->state = NEP_STATE_SOLVED;

100:   /* Only the first nconv columns contain useful information */
101:   BVSetActiveColumns(nep->V,0,nep->nconv);
102:   if (nep->twosided) BVSetActiveColumns(nep->W,0,nep->nconv);

104:   if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
105:     NEPComputeVectors(nep);
106:     NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv);
107:     nep->state = NEP_STATE_EIGENVECTORS;
108:   }

110:   /* sort eigenvalues according to nep->which parameter */
111:   SlepcSortEigenvalues(nep->sc,nep->nconv,nep->eigr,nep->eigi,nep->perm);
112:   PetscLogEventEnd(NEP_Solve,nep,0,0,0);

114:   /* various viewers */
115:   NEPViewFromOptions(nep,NULL,"-nep_view");
116:   NEPConvergedReasonViewFromOptions(nep);
117:   NEPErrorViewFromOptions(nep);
118:   NEPValuesViewFromOptions(nep);
119:   NEPVectorsViewFromOptions(nep);

121:   /* Remove the initial subspace */
122:   nep->nini = 0;

124:   /* Reset resolvent information */
125:   MatDestroy(&nep->resolvent);
126:   PetscFunctionReturn(0);
127: }

129: /*@
130:    NEPProjectOperator - Computes the projection of the nonlinear operator.

132:    Collective on nep

134:    Input Parameters:
135: +  nep - the nonlinear eigensolver context
136: .  j0  - initial index
137: -  j1  - final index

139:    Notes:
140:    This is available for split operator only.

142:    The nonlinear operator T(lambda) is projected onto span(V), where V is
143:    an orthonormal basis built internally by the solver. The projected
144:    operator is equal to sum_i V'*A_i*V*f_i(lambda), so this function
145:    computes all matrices Ei = V'*A_i*V, and stores them in the extra
146:    matrices inside DS. Only rows/columns in the range [j0,j1-1] are computed,
147:    the previous ones are assumed to be available already.

149:    Level: developer

151: .seealso: NEPSetSplitOperator()
152: @*/
153: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
154: {
155:   PetscInt       k;
156:   Mat            G;

161:   NEPCheckProblem(nep,1);
162:   NEPCheckSplit(nep,1);
163:   BVSetActiveColumns(nep->V,j0,j1);
164:   for (k=0;k<nep->nt;k++) {
165:     DSGetMat(nep->ds,DSMatExtra[k],&G);
166:     BVMatProject(nep->V,nep->A[k],nep->V,G);
167:     DSRestoreMat(nep->ds,DSMatExtra[k],&G);
168:   }
169:   PetscFunctionReturn(0);
170: }

172: /*@
173:    NEPApplyFunction - Applies the nonlinear function T(lambda) to a given vector.

175:    Collective on nep

177:    Input Parameters:
178: +  nep    - the nonlinear eigensolver context
179: .  lambda - scalar argument
180: .  x      - vector to be multiplied against
181: -  v      - workspace vector (used only in the case of split form)

183:    Output Parameters:
184: +  y   - result vector
185: .  A   - (optional) Function matrix, for callback interface only
186: -  B   - (unused) preconditioning matrix

188:    Note:
189:    If the nonlinear operator is represented in split form, the result
190:    y = T(lambda)*x is computed without building T(lambda) explicitly. In
191:    that case, parameters A and B are not used. Otherwise, the matrix
192:    T(lambda) is built and the effect is the same as a call to
193:    NEPComputeFunction() followed by a MatMult().

195:    Level: developer

197: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
198: @*/
199: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
200: {
201:   PetscInt       i;
202:   PetscScalar    alpha;


212:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
213:     VecSet(y,0.0);
214:     for (i=0;i<nep->nt;i++) {
215:       FNEvaluateFunction(nep->f[i],lambda,&alpha);
216:       MatMult(nep->A[i],x,v);
217:       VecAXPY(y,alpha,v);
218:     }
219:   } else {
220:     if (!A) A = nep->function;
221:     NEPComputeFunction(nep,lambda,A,A);
222:     MatMult(A,x,y);
223:   }
224:   PetscFunctionReturn(0);
225: }

227: /*@
228:    NEPApplyAdjoint - Applies the adjoint nonlinear function T(lambda)^* to a given vector.

230:    Collective on nep

232:    Input Parameters:
233: +  nep    - the nonlinear eigensolver context
234: .  lambda - scalar argument
235: .  x      - vector to be multiplied against
236: -  v      - workspace vector (used only in the case of split form)

238:    Output Parameters:
239: +  y   - result vector
240: .  A   - (optional) Function matrix, for callback interface only
241: -  B   - (unused) preconditioning matrix

243:    Level: developer

245: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
246: @*/
247: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
248: {
249:   PetscInt       i;
250:   PetscScalar    alpha;
251:   Vec            w;


261:   VecDuplicate(x,&w);
262:   VecCopy(x,w);
263:   VecConjugate(w);
264:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
265:     VecSet(y,0.0);
266:     for (i=0;i<nep->nt;i++) {
267:       FNEvaluateFunction(nep->f[i],lambda,&alpha);
268:       MatMultTranspose(nep->A[i],w,v);
269:       VecAXPY(y,alpha,v);
270:     }
271:   } else {
272:     if (!A) A = nep->function;
273:     NEPComputeFunction(nep,lambda,A,A);
274:     MatMultTranspose(A,w,y);
275:   }
276:   VecDestroy(&w);
277:   VecConjugate(y);
278:   PetscFunctionReturn(0);
279: }

281: /*@
282:    NEPApplyJacobian - Applies the nonlinear Jacobian T'(lambda) to a given vector.

284:    Collective on nep

286:    Input Parameters:
287: +  nep    - the nonlinear eigensolver context
288: .  lambda - scalar argument
289: .  x      - vector to be multiplied against
290: -  v      - workspace vector (used only in the case of split form)

292:    Output Parameters:
293: +  y   - result vector
294: -  A   - (optional) Jacobian matrix, for callback interface only

296:    Note:
297:    If the nonlinear operator is represented in split form, the result
298:    y = T'(lambda)*x is computed without building T'(lambda) explicitly. In
299:    that case, parameter A is not used. Otherwise, the matrix
300:    T'(lambda) is built and the effect is the same as a call to
301:    NEPComputeJacobian() followed by a MatMult().

303:    Level: developer

305: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
306: @*/
307: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
308: {
309:   PetscInt       i;
310:   PetscScalar    alpha;


319:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
320:     VecSet(y,0.0);
321:     for (i=0;i<nep->nt;i++) {
322:       FNEvaluateDerivative(nep->f[i],lambda,&alpha);
323:       MatMult(nep->A[i],x,v);
324:       VecAXPY(y,alpha,v);
325:     }
326:   } else {
327:     if (!A) A = nep->jacobian;
328:     NEPComputeJacobian(nep,lambda,A);
329:     MatMult(A,x,y);
330:   }
331:   PetscFunctionReturn(0);
332: }

334: /*@
335:    NEPGetIterationNumber - Gets the current iteration number. If the
336:    call to NEPSolve() is complete, then it returns the number of iterations
337:    carried out by the solution method.

339:    Not Collective

341:    Input Parameter:
342: .  nep - the nonlinear eigensolver context

344:    Output Parameter:
345: .  its - number of iterations

347:    Note:
348:    During the i-th iteration this call returns i-1. If NEPSolve() is
349:    complete, then parameter "its" contains either the iteration number at
350:    which convergence was successfully reached, or failure was detected.
351:    Call NEPGetConvergedReason() to determine if the solver converged or
352:    failed and why.

354:    Level: intermediate

356: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
357: @*/
358: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
359: {
362:   *its = nep->its;
363:   PetscFunctionReturn(0);
364: }

366: /*@
367:    NEPGetConverged - Gets the number of converged eigenpairs.

369:    Not Collective

371:    Input Parameter:
372: .  nep - the nonlinear eigensolver context

374:    Output Parameter:
375: .  nconv - number of converged eigenpairs

377:    Note:
378:    This function should be called after NEPSolve() has finished.

380:    Level: beginner

382: .seealso: NEPSetDimensions(), NEPSolve(), NEPGetEigenpair()
383: @*/
384: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
385: {
388:   NEPCheckSolved(nep,1);
389:   *nconv = nep->nconv;
390:   PetscFunctionReturn(0);
391: }

393: /*@
394:    NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
395:    stopped.

397:    Not Collective

399:    Input Parameter:
400: .  nep - the nonlinear eigensolver context

402:    Output Parameter:
403: .  reason - negative value indicates diverged, positive value converged

405:    Options Database Key:
406: .  -nep_converged_reason - print the reason to a viewer

408:    Notes:
409:    Possible values for reason are
410: +  NEP_CONVERGED_TOL - converged up to tolerance
411: .  NEP_CONVERGED_USER - converged due to a user-defined condition
412: .  NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
413: .  NEP_DIVERGED_BREAKDOWN - generic breakdown in method
414: .  NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
415: -  NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
416:    unrestarted solver

418:    Can only be called after the call to NEPSolve() is complete.

420:    Level: intermediate

422: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason
423: @*/
424: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
425: {
428:   NEPCheckSolved(nep,1);
429:   *reason = nep->reason;
430:   PetscFunctionReturn(0);
431: }

433: /*@C
434:    NEPGetEigenpair - Gets the i-th solution of the eigenproblem as computed by
435:    NEPSolve(). The solution consists in both the eigenvalue and the eigenvector.

437:    Logically Collective on nep

439:    Input Parameters:
440: +  nep - nonlinear eigensolver context
441: -  i   - index of the solution

443:    Output Parameters:
444: +  eigr - real part of eigenvalue
445: .  eigi - imaginary part of eigenvalue
446: .  Vr   - real part of eigenvector
447: -  Vi   - imaginary part of eigenvector

449:    Notes:
450:    It is allowed to pass NULL for Vr and Vi, if the eigenvector is not
451:    required. Otherwise, the caller must provide valid Vec objects, i.e.,
452:    they must be created by the calling program with e.g. MatCreateVecs().

454:    If the eigenvalue is real, then eigi and Vi are set to zero. If PETSc is
455:    configured with complex scalars the eigenvalue is stored
456:    directly in eigr (eigi is set to zero) and the eigenvector in Vr (Vi is
457:    set to zero). In any case, the user can pass NULL in Vr or Vi if one of
458:    them is not required.

460:    The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
461:    Eigenpairs are indexed according to the ordering criterion established
462:    with NEPSetWhichEigenpairs().

464:    Level: beginner

466: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
467: @*/
468: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
469: {
470:   PetscInt       k;

476:   NEPCheckSolved(nep,1);

480:   NEPComputeVectors(nep);
481:   k = nep->perm[i];

483:   /* eigenvalue */
484: #if defined(PETSC_USE_COMPLEX)
485:   if (eigr) *eigr = nep->eigr[k];
486:   if (eigi) *eigi = 0;
487: #else
488:   if (eigr) *eigr = nep->eigr[k];
489:   if (eigi) *eigi = nep->eigi[k];
490: #endif

492:   /* eigenvector */
493:   BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi);
494:   PetscFunctionReturn(0);
495: }

497: /*@
498:    NEPGetLeftEigenvector - Gets the i-th left eigenvector as computed by NEPSolve().

500:    Logically Collective on nep

502:    Input Parameters:
503: +  nep - eigensolver context
504: -  i   - index of the solution

506:    Output Parameters:
507: +  Wr   - real part of left eigenvector
508: -  Wi   - imaginary part of left eigenvector

510:    Notes:
511:    The caller must provide valid Vec objects, i.e., they must be created
512:    by the calling program with e.g. MatCreateVecs().

514:    If the corresponding eigenvalue is real, then Wi is set to zero. If PETSc is
515:    configured with complex scalars the eigenvector is stored directly in Wr
516:    (Wi is set to zero). In any case, the user can pass NULL in Wr or Wi if one of
517:    them is not required.

519:    The index i should be a value between 0 and nconv-1 (see NEPGetConverged()).
520:    Eigensolutions are indexed according to the ordering criterion established
521:    with NEPSetWhichEigenpairs().

523:    Left eigenvectors are available only if the twosided flag was set, see
524:    NEPSetTwoSided().

526:    Level: intermediate

528: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
529: @*/
530: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
531: {
532:   PetscInt       k;

538:   NEPCheckSolved(nep,1);
542:   NEPComputeVectors(nep);
543:   k = nep->perm[i];
544:   BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi);
545:   PetscFunctionReturn(0);
546: }

548: /*@
549:    NEPGetErrorEstimate - Returns the error estimate associated to the i-th
550:    computed eigenpair.

552:    Not Collective

554:    Input Parameters:
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 NEPComputeError().

565:    Level: advanced

567: .seealso: NEPComputeError()
568: @*/
569: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
570: {
573:   NEPCheckSolved(nep,1);
576:   *errest = nep->errest[nep->perm[i]];
577:   PetscFunctionReturn(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: {
592:   Vec            y,z=NULL;

594:   y = w[0];
595:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) z = w[1];
596:   if (adj) NEPApplyAdjoint(nep,lambda,x,z,y,NULL,NULL);
597:   else NEPApplyFunction(nep,lambda,x,z,y,NULL,NULL);
598:   VecNorm(y,NORM_2,norm);
599:   PetscFunctionReturn(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 Parameters:
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: {
627:   Vec            xr,xi=NULL;
628:   PetscInt       j,nwork,issplit=0;
629:   PetscScalar    kr,ki,s;
630:   PetscReal      er,z=0.0,errorl,nrm;
631:   PetscBool      flg;

637:   NEPCheckSolved(nep,1);

639:   /* allocate work vectors */
640: #if defined(PETSC_USE_COMPLEX)
641:   nwork = 2;
642: #else
643:   nwork = 3;
644: #endif
645:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
646:     issplit = 1;
647:     nwork++;  /* need an extra work vector for NEPComputeResidualNorm_Private */
648:   }
649:   NEPSetWorkVecs(nep,nwork);
650:   xr = nep->work[issplit+1];
651: #if !defined(PETSC_USE_COMPLEX)
652:   xi = nep->work[issplit+2];
653: #endif

655:   /* compute residual norms */
656:   NEPGetEigenpair(nep,i,&kr,&ki,xr,xi);
657: #if !defined(PETSC_USE_COMPLEX)
659: #endif
660:   NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error);
661:   VecNorm(xr,NORM_2,&er);

663:   /* if two-sided, compute left residual norm and take the maximum */
664:   if (nep->twosided) {
665:     NEPGetLeftEigenvector(nep,i,xr,xi);
666:     NEPComputeResidualNorm_Private(nep,PETSC_TRUE,kr,xr,nep->work,&errorl);
667:     *error = PetscMax(*error,errorl);
668:   }

670:   /* compute error */
671:   switch (type) {
672:     case NEP_ERROR_ABSOLUTE:
673:       break;
674:     case NEP_ERROR_RELATIVE:
675:       *error /= PetscAbsScalar(kr)*er;
676:       break;
677:     case NEP_ERROR_BACKWARD:
678:       if (nep->fui!=NEP_USER_INTERFACE_SPLIT) {
679:         NEPComputeFunction(nep,kr,nep->function,nep->function);
680:         MatHasOperation(nep->function,MATOP_NORM,&flg);
682:         MatNorm(nep->function,NORM_INFINITY,&nrm);
683:         *error /= nrm*er;
684:         break;
685:       }
686:       /* initialization of matrix norms */
687:       if (!nep->nrma[0]) {
688:         for (j=0;j<nep->nt;j++) {
689:           MatHasOperation(nep->A[j],MATOP_NORM,&flg);
691:           MatNorm(nep->A[j],NORM_INFINITY,&nep->nrma[j]);
692:         }
693:       }
694:       for (j=0;j<nep->nt;j++) {
695:         FNEvaluateFunction(nep->f[j],kr,&s);
696:         z = z + nep->nrma[j]*PetscAbsScalar(s);
697:       }
698:       *error /= z*er;
699:       break;
700:     default:
701:       SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
702:   }
703:   PetscFunctionReturn(0);
704: }

706: /*@
707:    NEPComputeFunction - Computes the function matrix T(lambda) that has been
708:    set with NEPSetFunction().

710:    Collective on nep

712:    Input Parameters:
713: +  nep    - the NEP context
714: -  lambda - the scalar argument

716:    Output Parameters:
717: +  A   - Function matrix
718: -  B   - optional preconditioning matrix

720:    Notes:
721:    NEPComputeFunction() is typically used within nonlinear eigensolvers
722:    implementations, so most users would not generally call this routine
723:    themselves.

725:    Level: developer

727: .seealso: NEPSetFunction(), NEPGetFunction()
728: @*/
729: PetscErrorCode NEPComputeFunction(NEP nep,PetscScalar lambda,Mat A,Mat B)
730: {
731:   PetscInt       i;
732:   PetscScalar    alpha;

735:   NEPCheckProblem(nep,1);
736:   switch (nep->fui) {
737:   case NEP_USER_INTERFACE_CALLBACK:
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:     if (A != B) MatZeroEntries(B);
748:     for (i=0;i<nep->nt;i++) {
749:       FNEvaluateFunction(nep->f[i],lambda,&alpha);
750:       MatAXPY(A,alpha,nep->A[i],nep->mstr);
751:       if (A != B) MatAXPY(B,alpha,nep->P[i],nep->mstrp);
752:     }
753:     break;
754:   }
755:   PetscFunctionReturn(0);
756: }

758: /*@
759:    NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
760:    set with NEPSetJacobian().

762:    Collective on nep

764:    Input Parameters:
765: +  nep    - the NEP context
766: -  lambda - the scalar argument

768:    Output Parameters:
769: .  A   - Jacobian matrix

771:    Notes:
772:    Most users should not need to explicitly call this routine, as it
773:    is used internally within the nonlinear eigensolvers.

775:    Level: developer

777: .seealso: NEPSetJacobian(), NEPGetJacobian()
778: @*/
779: PetscErrorCode NEPComputeJacobian(NEP nep,PetscScalar lambda,Mat A)
780: {
781:   PetscInt       i;
782:   PetscScalar    alpha;

785:   NEPCheckProblem(nep,1);
786:   switch (nep->fui) {
787:   case NEP_USER_INTERFACE_CALLBACK:
789:     PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0);
790:     PetscStackPush("NEP user Jacobian function");
791:     (*nep->computejacobian)(nep,lambda,A,nep->jacobianctx);
792:     PetscStackPop;
793:     PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0);
794:     break;
795:   case NEP_USER_INTERFACE_SPLIT:
796:     MatZeroEntries(A);
797:     for (i=0;i<nep->nt;i++) {
798:       FNEvaluateDerivative(nep->f[i],lambda,&alpha);
799:       MatAXPY(A,alpha,nep->A[i],nep->mstr);
800:     }
801:     break;
802:   }
803:   PetscFunctionReturn(0);
804: }