Actual source code: nepsolve.c

slepc-3.13.4 2020-09-02
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
 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(), NEPConvergedReason
411: @*/
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: }