Actual source code: nepsolve.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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 NEP

109:    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 NEP

154:    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 NEP

210:    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 NEP

262:    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(), NEPConvergedReason
400: @*/
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 NEP

417:    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 NEP

481:    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 NEP

593:    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: }