Actual source code: nepsolve.c

slepc-3.14.0 2020-09-30
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:   /* call solver */
 86:   (*nep->ops->solve)(nep);
 87:   if (!nep->reason) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_PLIB,"Internal error, solver returned without setting converged reason");
 88:   nep->state = NEP_STATE_SOLVED;

 90:   /* Only the first nconv columns contain useful information */
 91:   BVSetActiveColumns(nep->V,0,nep->nconv);
 92:   if (nep->twosided) { BVSetActiveColumns(nep->W,0,nep->nconv); }

 94:   if (nep->refine==NEP_REFINE_SIMPLE && nep->rits>0 && nep->nconv>0) {
 95:     NEPComputeVectors(nep);
 96:     NEPNewtonRefinementSimple(nep,&nep->rits,nep->rtol,nep->nconv);
 97:     nep->state = NEP_STATE_EIGENVECTORS;
 98:   }

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

104:   /* various viewers */
105:   NEPViewFromOptions(nep,NULL,"-nep_view");
106:   NEPConvergedReasonViewFromOptions(nep);
107:   NEPErrorViewFromOptions(nep);
108:   NEPValuesViewFromOptions(nep);
109:   NEPVectorsViewFromOptions(nep);

111:   /* Remove the initial subspace */
112:   nep->nini = 0;

114:   /* Reset resolvent information */
115:   MatDestroy(&nep->resolvent);
116:   return(0);
117: }

119: /*@
120:    NEPProjectOperator - Computes the projection of the nonlinear operator.

122:    Collective on nep

124:    Input Parameters:
125: +  nep - the nonlinear eigensolver context
126: .  j0  - initial index
127: -  j1  - final index

129:    Notes:
130:    This is available for split operator only.

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

139:    Level: developer

141: .seealso: NEPSetSplitOperator()
142: @*/
143: PetscErrorCode NEPProjectOperator(NEP nep,PetscInt j0,PetscInt j1)
144: {
146:   PetscInt       k;
147:   Mat            G;

153:   NEPCheckProblem(nep,1);
154:   NEPCheckSplit(nep,1);
155:   BVSetActiveColumns(nep->V,j0,j1);
156:   for (k=0;k<nep->nt;k++) {
157:     DSGetMat(nep->ds,DSMatExtra[k],&G);
158:     BVMatProject(nep->V,nep->A[k],nep->V,G);
159:     DSRestoreMat(nep->ds,DSMatExtra[k],&G);
160:   }
161:   return(0);
162: }

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

167:    Collective on nep

169:    Input Parameters:
170: +  nep    - the nonlinear eigensolver context
171: .  lambda - scalar argument
172: .  x      - vector to be multiplied against
173: -  v      - workspace vector (used only in the case of split form)

175:    Output Parameters:
176: +  y   - result vector
177: .  A   - Function matrix
178: -  B   - optional preconditioning matrix

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

187:    Level: developer

189: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyAdjoint()
190: @*/
191: PetscErrorCode NEPApplyFunction(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
192: {
194:   PetscInt       i;
195:   PetscScalar    alpha;


206:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
207:     VecSet(y,0.0);
208:     for (i=0;i<nep->nt;i++) {
209:       FNEvaluateFunction(nep->f[i],lambda,&alpha);
210:       MatMult(nep->A[i],x,v);
211:       VecAXPY(y,alpha,v);
212:     }
213:   } else {
214:     NEPComputeFunction(nep,lambda,A,B);
215:     MatMult(A,x,y);
216:   }
217:   return(0);
218: }

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

223:    Collective on nep

225:    Input Parameters:
226: +  nep    - the nonlinear eigensolver context
227: .  lambda - scalar argument
228: .  x      - vector to be multiplied against
229: -  v      - workspace vector (used only in the case of split form)

231:    Output Parameters:
232: +  y   - result vector
233: .  A   - Function matrix
234: -  B   - optional preconditioning matrix

236:    Level: developer

238: .seealso: NEPSetSplitOperator(), NEPComputeFunction(), NEPApplyFunction()
239: @*/
240: PetscErrorCode NEPApplyAdjoint(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A,Mat B)
241: {
243:   PetscInt       i;
244:   PetscScalar    alpha;


255:   VecConjugate(x);
256:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
257:     VecSet(y,0.0);
258:     for (i=0;i<nep->nt;i++) {
259:       FNEvaluateFunction(nep->f[i],lambda,&alpha);
260:       MatMultTranspose(nep->A[i],x,v);
261:       VecAXPY(y,alpha,v);
262:     }
263:   } else {
264:     NEPComputeFunction(nep,lambda,A,B);
265:     MatMultTranspose(A,x,y);
266:   }
267:   VecConjugate(x);
268:   VecConjugate(y);
269:   return(0);
270: }

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

275:    Collective on nep

277:    Input Parameters:
278: +  nep    - the nonlinear eigensolver context
279: .  lambda - scalar argument
280: .  x      - vector to be multiplied against
281: -  v      - workspace vector (used only in the case of split form)

283:    Output Parameters:
284: +  y   - result vector
285: -  A   - Jacobian matrix

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

294:    Level: developer

296: .seealso: NEPSetSplitOperator(), NEPComputeJacobian()
297: @*/
298: PetscErrorCode NEPApplyJacobian(NEP nep,PetscScalar lambda,Vec x,Vec v,Vec y,Mat A)
299: {
301:   PetscInt       i;
302:   PetscScalar    alpha;


312:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
313:     VecSet(y,0.0);
314:     for (i=0;i<nep->nt;i++) {
315:       FNEvaluateDerivative(nep->f[i],lambda,&alpha);
316:       MatMult(nep->A[i],x,v);
317:       VecAXPY(y,alpha,v);
318:     }
319:   } else {
320:     NEPComputeJacobian(nep,lambda,A);
321:     MatMult(A,x,y);
322:   }
323:   return(0);
324: }

326: /*@
327:    NEPGetIterationNumber - Gets the current iteration number. If the
328:    call to NEPSolve() is complete, then it returns the number of iterations
329:    carried out by the solution method.

331:    Not Collective

333:    Input Parameter:
334: .  nep - the nonlinear eigensolver context

336:    Output Parameter:
337: .  its - number of iterations

339:    Note:
340:    During the i-th iteration this call returns i-1. If NEPSolve() is
341:    complete, then parameter "its" contains either the iteration number at
342:    which convergence was successfully reached, or failure was detected.
343:    Call NEPGetConvergedReason() to determine if the solver converged or
344:    failed and why.

346:    Level: intermediate

348: .seealso: NEPGetConvergedReason(), NEPSetTolerances()
349: @*/
350: PetscErrorCode NEPGetIterationNumber(NEP nep,PetscInt *its)
351: {
355:   *its = nep->its;
356:   return(0);
357: }

359: /*@
360:    NEPGetConverged - Gets the number of converged eigenpairs.

362:    Not Collective

364:    Input Parameter:
365: .  nep - the nonlinear eigensolver context

367:    Output Parameter:
368: .  nconv - number of converged eigenpairs

370:    Note:
371:    This function should be called after NEPSolve() has finished.

373:    Level: beginner

375: .seealso: NEPSetDimensions(), NEPSolve()
376: @*/
377: PetscErrorCode NEPGetConverged(NEP nep,PetscInt *nconv)
378: {
382:   NEPCheckSolved(nep,1);
383:   *nconv = nep->nconv;
384:   return(0);
385: }

387: /*@
388:    NEPGetConvergedReason - Gets the reason why the NEPSolve() iteration was
389:    stopped.

391:    Not Collective

393:    Input Parameter:
394: .  nep - the nonlinear eigensolver context

396:    Output Parameter:
397: .  reason - negative value indicates diverged, positive value converged

399:    Notes:

401:    Possible values for reason are
402: +  NEP_CONVERGED_TOL - converged up to tolerance
403: .  NEP_CONVERGED_USER - converged due to a user-defined condition
404: .  NEP_DIVERGED_ITS - required more than max_it iterations to reach convergence
405: .  NEP_DIVERGED_BREAKDOWN - generic breakdown in method
406: .  NEP_DIVERGED_LINEAR_SOLVE - inner linear solve failed
407: -  NEP_DIVERGED_SUBSPACE_EXHAUSTED - run out of space for the basis in an
408:    unrestarted solver

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

412:    Level: intermediate

414: .seealso: NEPSetTolerances(), NEPSolve(), NEPConvergedReason
415: @*/
416: PetscErrorCode NEPGetConvergedReason(NEP nep,NEPConvergedReason *reason)
417: {
421:   NEPCheckSolved(nep,1);
422:   *reason = nep->reason;
423:   return(0);
424: }

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

430:    Logically Collective on nep

432:    Input Parameters:
433: +  nep - nonlinear eigensolver context
434: -  i   - index of the solution

436:    Output Parameters:
437: +  eigr - real part of eigenvalue
438: .  eigi - imaginary part of eigenvalue
439: .  Vr   - real part of eigenvector
440: -  Vi   - imaginary part of eigenvector

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

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

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

457:    Level: beginner

459: .seealso: NEPSolve(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPGetLeftEigenvector()
460: @*/
461: PetscErrorCode NEPGetEigenpair(NEP nep,PetscInt i,PetscScalar *eigr,PetscScalar *eigi,Vec Vr,Vec Vi)
462: {
463:   PetscInt       k;

471:   NEPCheckSolved(nep,1);
472:   if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");

474:   NEPComputeVectors(nep);
475:   k = nep->perm[i];

477:   /* eigenvalue */
478: #if defined(PETSC_USE_COMPLEX)
479:   if (eigr) *eigr = nep->eigr[k];
480:   if (eigi) *eigi = 0;
481: #else
482:   if (eigr) *eigr = nep->eigr[k];
483:   if (eigi) *eigi = nep->eigi[k];
484: #endif

486:   /* eigenvector */
487:   BV_GetEigenvector(nep->V,k,nep->eigi[k],Vr,Vi);
488:   return(0);
489: }

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

494:    Logically Collective on nep

496:    Input Parameters:
497: +  nep - eigensolver context
498: -  i   - index of the solution

500:    Output Parameters:
501: +  Wr   - real part of left eigenvector
502: -  Wi   - imaginary part of left eigenvector

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

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

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

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

520:    Level: intermediate

522: .seealso: NEPGetEigenpair(), NEPGetConverged(), NEPSetWhichEigenpairs(), NEPSetTwoSided()
523: @*/
524: PetscErrorCode NEPGetLeftEigenvector(NEP nep,PetscInt i,Vec Wr,Vec Wi)
525: {
527:   PetscInt       k;

534:   NEPCheckSolved(nep,1);
535:   if (!nep->twosided) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONGSTATE,"Must request left vectors with NEPSetTwoSided");
536:   if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
537:   NEPComputeVectors(nep);
538:   k = nep->perm[i];
539:   BV_GetEigenvector(nep->W,k,nep->eigi[k],Wr,Wi);
540:   return(0);
541: }

543: /*@
544:    NEPGetErrorEstimate - Returns the error estimate associated to the i-th
545:    computed eigenpair.

547:    Not Collective

549:    Input Parameter:
550: +  nep - nonlinear eigensolver context
551: -  i   - index of eigenpair

553:    Output Parameter:
554: .  errest - the error estimate

556:    Notes:
557:    This is the error estimate used internally by the eigensolver. The actual
558:    error bound can be computed with NEPComputeRelativeError().

560:    Level: advanced

562: .seealso: NEPComputeRelativeError()
563: @*/
564: PetscErrorCode NEPGetErrorEstimate(NEP nep,PetscInt i,PetscReal *errest)
565: {
569:   NEPCheckSolved(nep,1);
570:   if (i<0 || i>=nep->nconv) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Argument 2 out of range");
571:   *errest = nep->errest[nep->perm[i]];
572:   return(0);
573: }

575: /*
576:    NEPComputeResidualNorm_Private - Computes the norm of the residual vector
577:    associated with an eigenpair.

579:    Input Parameters:
580:      adj    - whether the adjoint T^* must be used instead of T
581:      lambda - eigenvalue
582:      x      - eigenvector
583:      w      - array of work vectors (two vectors in split form, one vector otherwise)
584: */
585: PetscErrorCode NEPComputeResidualNorm_Private(NEP nep,PetscBool adj,PetscScalar lambda,Vec x,Vec *w,PetscReal *norm)
586: {
588:   Vec            y,z=NULL;

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

639:   NEPCheckSolved(nep,1);

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

657:   /* compute residual norms */
658:   NEPGetEigenpair(nep,i,&kr,&ki,xr,xi);
659: #if !defined(PETSC_USE_COMPLEX)
660:   if (ki) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented for complex eigenvalues with real scalars");
661: #endif
662:   NEPComputeResidualNorm_Private(nep,PETSC_FALSE,kr,xr,nep->work,error);
663:   VecNorm(xr,NORM_2,&er);

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

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

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

712:    Collective on nep

714:    Input Parameters:
715: +  nep    - the NEP context
716: -  lambda - the scalar argument

718:    Output Parameters:
719: +  A   - Function matrix
720: -  B   - optional preconditioning matrix

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

727:    Level: developer

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

739:   NEPCheckProblem(nep,1);
740:   switch (nep->fui) {
741:   case NEP_USER_INTERFACE_CALLBACK:
742:     if (!nep->computefunction) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetFunction() first");
743:     PetscLogEventBegin(NEP_FunctionEval,nep,A,B,0);
744:     PetscStackPush("NEP user Function function");
745:     (*nep->computefunction)(nep,lambda,A,B,nep->functionctx);
746:     PetscStackPop;
747:     PetscLogEventEnd(NEP_FunctionEval,nep,A,B,0);
748:     break;
749:   case NEP_USER_INTERFACE_SPLIT:
750:     MatZeroEntries(A);
751:     for (i=0;i<nep->nt;i++) {
752:       FNEvaluateFunction(nep->f[i],lambda,&alpha);
753:       MatAXPY(A,alpha,nep->A[i],nep->mstr);
754:     }
755:     if (A != B) SETERRQ(PetscObjectComm((PetscObject)nep),1,"Not implemented");
756:     break;
757:   }
758:   return(0);
759: }

761: /*@
762:    NEPComputeJacobian - Computes the Jacobian matrix T'(lambda) that has been
763:    set with NEPSetJacobian().

765:    Collective on nep

767:    Input Parameters:
768: +  nep    - the NEP context
769: -  lambda - the scalar argument

771:    Output Parameters:
772: .  A   - Jacobian matrix

774:    Notes:
775:    Most users should not need to explicitly call this routine, as it
776:    is used internally within the nonlinear eigensolvers.

778:    Level: developer

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

790:   NEPCheckProblem(nep,1);
791:   switch (nep->fui) {
792:   case NEP_USER_INTERFACE_CALLBACK:
793:     if (!nep->computejacobian) SETERRQ(PetscObjectComm((PetscObject)nep),PETSC_ERR_USER,"Must call NEPSetJacobian() first");
794:     PetscLogEventBegin(NEP_JacobianEval,nep,A,0,0);
795:     PetscStackPush("NEP user Jacobian function");
796:     (*nep->computejacobian)(nep,lambda,A,nep->jacobianctx);
797:     PetscStackPop;
798:     PetscLogEventEnd(NEP_JacobianEval,nep,A,0,0);
799:     break;
800:   case NEP_USER_INTERFACE_SPLIT:
801:     MatZeroEntries(A);
802:     for (i=0;i<nep->nt;i++) {
803:       FNEvaluateDerivative(nep->f[i],lambda,&alpha);
804:       MatAXPY(A,alpha,nep->A[i],nep->mstr);
805:     }
806:     break;
807:   }
808:   return(0);
809: }