Actual source code: nepsolve.c

slepc-3.14.1 2020-12-08
Report Typos and Errors
  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(), NEPConvergedReason
420: @*/
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: }