Actual source code: nepbasic.c

slepc-3.10.1 2018-10-23
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2018, 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:    Basic NEP routines
 12: */

 14: #include <slepc/private/nepimpl.h>      /*I "slepcnep.h" I*/

 16: PetscFunctionList NEPList = 0;
 17: PetscBool         NEPRegisterAllCalled = PETSC_FALSE;
 18: PetscClassId      NEP_CLASSID = 0;
 19: PetscLogEvent     NEP_SetUp = 0,NEP_Solve = 0,NEP_Refine = 0,NEP_FunctionEval = 0,NEP_JacobianEval = 0,NEP_DerivativesEval = 0;

 21: /*@
 22:    NEPCreate - Creates the default NEP context.

 24:    Collective on MPI_Comm

 26:    Input Parameter:
 27: .  comm - MPI communicator

 29:    Output Parameter:
 30: .  nep - location to put the NEP context

 32:    Level: beginner

 34: .seealso: NEPSetUp(), NEPSolve(), NEPDestroy(), NEP
 35: @*/
 36: PetscErrorCode NEPCreate(MPI_Comm comm,NEP *outnep)
 37: {
 39:   NEP            nep;

 43:   *outnep = 0;
 44:   NEPInitializePackage();
 45:   SlepcHeaderCreate(nep,NEP_CLASSID,"NEP","Nonlinear Eigenvalue Problem","NEP",comm,NEPDestroy,NEPView);

 47:   nep->max_it          = 0;
 48:   nep->nev             = 1;
 49:   nep->ncv             = 0;
 50:   nep->mpd             = 0;
 51:   nep->nini            = 0;
 52:   nep->target          = 0.0;
 53:   nep->tol             = PETSC_DEFAULT;
 54:   nep->conv            = NEP_CONV_REL;
 55:   nep->stop            = NEP_STOP_BASIC;
 56:   nep->which           = (NEPWhich)0;
 57:   nep->problem_type    = (NEPProblemType)0;
 58:   nep->refine          = NEP_REFINE_NONE;
 59:   nep->npart           = 1;
 60:   nep->rtol            = PETSC_DEFAULT;
 61:   nep->rits            = PETSC_DEFAULT;
 62:   nep->scheme          = (NEPRefineScheme)0;
 63:   nep->trackall        = PETSC_FALSE;

 65:   nep->computefunction = NULL;
 66:   nep->computejacobian = NULL;
 67:   nep->functionctx     = NULL;
 68:   nep->jacobianctx     = NULL;
 69:   nep->computederivatives = NULL;
 70:   nep->derivativesctx  = NULL;
 71:   nep->converged       = NEPConvergedRelative;
 72:   nep->convergeduser   = NULL;
 73:   nep->convergeddestroy= NULL;
 74:   nep->stopping        = NEPStoppingBasic;
 75:   nep->stoppinguser    = NULL;
 76:   nep->stoppingdestroy = NULL;
 77:   nep->convergedctx    = NULL;
 78:   nep->stoppingctx     = NULL;
 79:   nep->numbermonitors  = 0;

 81:   nep->ds              = NULL;
 82:   nep->V               = NULL;
 83:   nep->rg              = NULL;
 84:   nep->function        = NULL;
 85:   nep->function_pre    = NULL;
 86:   nep->jacobian        = NULL;
 87:   nep->derivatives     = NULL;
 88:   nep->A               = NULL;
 89:   nep->f               = NULL;
 90:   nep->nt              = 0;
 91:   nep->mstr            = DIFFERENT_NONZERO_PATTERN;
 92:   nep->IS              = NULL;
 93:   nep->eigr            = NULL;
 94:   nep->eigi            = NULL;
 95:   nep->errest          = NULL;
 96:   nep->perm            = NULL;
 97:   nep->nwork           = 0;
 98:   nep->work            = NULL;
 99:   nep->data            = NULL;

101:   nep->state           = NEP_STATE_INITIAL;
102:   nep->nconv           = 0;
103:   nep->its             = 0;
104:   nep->n               = 0;
105:   nep->nloc            = 0;
106:   nep->nrma            = NULL;
107:   nep->fui             = (NEPUserInterface)0;
108:   nep->reason          = NEP_CONVERGED_ITERATING;

110:   PetscNewLog(nep,&nep->sc);
111:   *outnep = nep;
112:   return(0);
113: }

115: /*@C
116:    NEPSetType - Selects the particular solver to be used in the NEP object.

118:    Logically Collective on NEP

120:    Input Parameters:
121: +  nep      - the nonlinear eigensolver context
122: -  type     - a known method

124:    Options Database Key:
125: .  -nep_type <method> - Sets the method; use -help for a list
126:     of available methods

128:    Notes:
129:    See "slepc/include/slepcnep.h" for available methods.

131:    Normally, it is best to use the NEPSetFromOptions() command and
132:    then set the NEP type from the options database rather than by using
133:    this routine.  Using the options database provides the user with
134:    maximum flexibility in evaluating the different available methods.
135:    The NEPSetType() routine is provided for those situations where it
136:    is necessary to set the iterative solver independently of the command
137:    line or options database.

139:    Level: intermediate

141: .seealso: NEPType
142: @*/
143: PetscErrorCode NEPSetType(NEP nep,NEPType type)
144: {
145:   PetscErrorCode ierr,(*r)(NEP);
146:   PetscBool      match;


152:   PetscObjectTypeCompare((PetscObject)nep,type,&match);
153:   if (match) return(0);

155:   PetscFunctionListFind(NEPList,type,&r);
156:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown NEP type given: %s",type);

158:   if (nep->ops->destroy) { (*nep->ops->destroy)(nep); }
159:   PetscMemzero(nep->ops,sizeof(struct _NEPOps));

161:   nep->state = NEP_STATE_INITIAL;
162:   PetscObjectChangeTypeName((PetscObject)nep,type);
163:   (*r)(nep);
164:   return(0);
165: }

167: /*@C
168:    NEPGetType - Gets the NEP type as a string from the NEP object.

170:    Not Collective

172:    Input Parameter:
173: .  nep - the eigensolver context

175:    Output Parameter:
176: .  name - name of NEP method

178:    Level: intermediate

180: .seealso: NEPSetType()
181: @*/
182: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
183: {
187:   *type = ((PetscObject)nep)->type_name;
188:   return(0);
189: }

191: /*@C
192:    NEPRegister - Adds a method to the nonlinear eigenproblem solver package.

194:    Not Collective

196:    Input Parameters:
197: +  name - name of a new user-defined solver
198: -  function - routine to create the solver context

200:    Notes:
201:    NEPRegister() may be called multiple times to add several user-defined solvers.

203:    Sample usage:
204: .vb
205:     NEPRegister("my_solver",MySolverCreate);
206: .ve

208:    Then, your solver can be chosen with the procedural interface via
209: $     NEPSetType(nep,"my_solver")
210:    or at runtime via the option
211: $     -nep_type my_solver

213:    Level: advanced

215: .seealso: NEPRegisterAll()
216: @*/
217: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
218: {

222:   NEPInitializePackage();
223:   PetscFunctionListAdd(&NEPList,name,function);
224:   return(0);
225: }

227: /*
228:    NEPReset_Problem - Destroys the problem matrices.
229: @*/
230: PetscErrorCode NEPReset_Problem(NEP nep)
231: {
233:   PetscInt       i;

237:   MatDestroy(&nep->function);
238:   MatDestroy(&nep->function_pre);
239:   MatDestroy(&nep->jacobian);
240:   MatDestroy(&nep->derivatives);
241:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
242:     MatDestroyMatrices(nep->nt,&nep->A);
243:     for (i=0;i<nep->nt;i++) {
244:       FNDestroy(&nep->f[i]);
245:     }
246:     PetscFree(nep->f);
247:     PetscFree(nep->nrma);
248:     nep->nt = 0;
249:   }
250:   return(0);
251: }
252: /*@
253:    NEPReset - Resets the NEP context to the initial state (prior to setup)
254:    and destroys any allocated Vecs and Mats.

256:    Collective on NEP

258:    Input Parameter:
259: .  nep - eigensolver context obtained from NEPCreate()

261:    Level: advanced

263: .seealso: NEPDestroy()
264: @*/
265: PetscErrorCode NEPReset(NEP nep)
266: {

271:   if (!nep) return(0);
272:   if (nep->ops->reset) { (nep->ops->reset)(nep); }
273:   if (nep->refineksp) { KSPReset(nep->refineksp); }
274:   NEPReset_Problem(nep);
275:   BVDestroy(&nep->V);
276:   VecDestroyVecs(nep->nwork,&nep->work);
277:   nep->nwork = 0;
278:   nep->state = NEP_STATE_INITIAL;
279:   return(0);
280: }

282: /*@
283:    NEPDestroy - Destroys the NEP context.

285:    Collective on NEP

287:    Input Parameter:
288: .  nep - eigensolver context obtained from NEPCreate()

290:    Level: beginner

292: .seealso: NEPCreate(), NEPSetUp(), NEPSolve()
293: @*/
294: PetscErrorCode NEPDestroy(NEP *nep)
295: {

299:   if (!*nep) return(0);
301:   if (--((PetscObject)(*nep))->refct > 0) { *nep = 0; return(0); }
302:   NEPReset(*nep);
303:   if ((*nep)->ops->destroy) { (*(*nep)->ops->destroy)(*nep); }
304:   if ((*nep)->eigr) {
305:     PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm);
306:   }
307:   RGDestroy(&(*nep)->rg);
308:   DSDestroy(&(*nep)->ds);
309:   KSPDestroy(&(*nep)->refineksp);
310:   PetscSubcommDestroy(&(*nep)->refinesubc);
311:   PetscFree((*nep)->sc);
312:   /* just in case the initial vectors have not been used */
313:   SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS);
314:   if ((*nep)->convergeddestroy) {
315:     (*(*nep)->convergeddestroy)((*nep)->convergedctx);
316:   }
317:   NEPMonitorCancel(*nep);
318:   PetscHeaderDestroy(nep);
319:   return(0);
320: }

322: /*@
323:    NEPSetBV - Associates a basis vectors object to the nonlinear eigensolver.

325:    Collective on NEP

327:    Input Parameters:
328: +  nep - eigensolver context obtained from NEPCreate()
329: -  bv  - the basis vectors object

331:    Note:
332:    Use NEPGetBV() to retrieve the basis vectors context (for example,
333:    to free it at the end of the computations).

335:    Level: advanced

337: .seealso: NEPGetBV()
338: @*/
339: PetscErrorCode NEPSetBV(NEP nep,BV bv)
340: {

347:   PetscObjectReference((PetscObject)bv);
348:   BVDestroy(&nep->V);
349:   nep->V = bv;
350:   PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->V);
351:   return(0);
352: }

354: /*@
355:    NEPGetBV - Obtain the basis vectors object associated to the nonlinear
356:    eigensolver object.

358:    Not Collective

360:    Input Parameters:
361: .  nep - eigensolver context obtained from NEPCreate()

363:    Output Parameter:
364: .  bv - basis vectors context

366:    Level: advanced

368: .seealso: NEPSetBV()
369: @*/
370: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
371: {

377:   if (!nep->V) {
378:     BVCreate(PetscObjectComm((PetscObject)nep),&nep->V);
379:     PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->V);
380:   }
381:   *bv = nep->V;
382:   return(0);
383: }

385: /*@
386:    NEPSetRG - Associates a region object to the nonlinear eigensolver.

388:    Collective on NEP

390:    Input Parameters:
391: +  nep - eigensolver context obtained from NEPCreate()
392: -  rg  - the region object

394:    Note:
395:    Use NEPGetRG() to retrieve the region context (for example,
396:    to free it at the end of the computations).

398:    Level: advanced

400: .seealso: NEPGetRG()
401: @*/
402: PetscErrorCode NEPSetRG(NEP nep,RG rg)
403: {

410:   PetscObjectReference((PetscObject)rg);
411:   RGDestroy(&nep->rg);
412:   nep->rg = rg;
413:   PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->rg);
414:   return(0);
415: }

417: /*@
418:    NEPGetRG - Obtain the region object associated to the
419:    nonlinear eigensolver object.

421:    Not Collective

423:    Input Parameters:
424: .  nep - eigensolver context obtained from NEPCreate()

426:    Output Parameter:
427: .  rg - region context

429:    Level: advanced

431: .seealso: NEPSetRG()
432: @*/
433: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
434: {

440:   if (!nep->rg) {
441:     RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg);
442:     PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->rg);
443:   }
444:   *rg = nep->rg;
445:   return(0);
446: }

448: /*@
449:    NEPSetDS - Associates a direct solver object to the nonlinear eigensolver.

451:    Collective on NEP

453:    Input Parameters:
454: +  nep - eigensolver context obtained from NEPCreate()
455: -  ds  - the direct solver object

457:    Note:
458:    Use NEPGetDS() to retrieve the direct solver context (for example,
459:    to free it at the end of the computations).

461:    Level: advanced

463: .seealso: NEPGetDS()
464: @*/
465: PetscErrorCode NEPSetDS(NEP nep,DS ds)
466: {

473:   PetscObjectReference((PetscObject)ds);
474:   DSDestroy(&nep->ds);
475:   nep->ds = ds;
476:   PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->ds);
477:   return(0);
478: }

480: /*@
481:    NEPGetDS - Obtain the direct solver object associated to the
482:    nonlinear eigensolver object.

484:    Not Collective

486:    Input Parameters:
487: .  nep - eigensolver context obtained from NEPCreate()

489:    Output Parameter:
490: .  ds - direct solver context

492:    Level: advanced

494: .seealso: NEPSetDS()
495: @*/
496: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
497: {

503:   if (!nep->ds) {
504:     DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds);
505:     PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->ds);
506:   }
507:   *ds = nep->ds;
508:   return(0);
509: }

511: /*@
512:    NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
513:    object in the refinement phase.

515:    Not Collective

517:    Input Parameters:
518: .  nep - eigensolver context obtained from NEPCreate()

520:    Output Parameter:
521: .  ksp - ksp context

523:    Level: advanced

525: .seealso: NEPSetRefine()
526: @*/
527: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
528: {

534:   if (!nep->refineksp) {
535:     if (nep->npart>1) {
536:       /* Split in subcomunicators */
537:       PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc);
538:       PetscSubcommSetNumber(nep->refinesubc,nep->npart);
539:       PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
540:       PetscLogObjectMemory((PetscObject)nep,sizeof(PetscSubcomm));
541:     }
542:     KSPCreate((nep->npart==1)?PetscObjectComm((PetscObject)nep):PetscSubcommChild(nep->refinesubc),&nep->refineksp);
543:     PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->refineksp);
544:     KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix);
545:     KSPAppendOptionsPrefix(*ksp,"nep_refine_");
546:     KSPSetTolerances(nep->refineksp,SLEPC_DEFAULT_TOL,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
547:   }
548:   *ksp = nep->refineksp;
549:   return(0);
550: }

552: /*@
553:    NEPSetTarget - Sets the value of the target.

555:    Logically Collective on NEP

557:    Input Parameters:
558: +  nep    - eigensolver context
559: -  target - the value of the target

561:    Options Database Key:
562: .  -nep_target <scalar> - the value of the target

564:    Notes:
565:    The target is a scalar value used to determine the portion of the spectrum
566:    of interest. It is used in combination with NEPSetWhichEigenpairs().

568:    In the case of complex scalars, a complex value can be provided in the
569:    command line with [+/-][realnumber][+/-]realnumberi with no spaces, e.g.
570:    -nep_target 1.0+2.0i

572:    Level: intermediate

574: .seealso: NEPGetTarget(), NEPSetWhichEigenpairs()
575: @*/
576: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
577: {
581:   nep->target = target;
582:   return(0);
583: }

585: /*@
586:    NEPGetTarget - Gets the value of the target.

588:    Not Collective

590:    Input Parameter:
591: .  nep - eigensolver context

593:    Output Parameter:
594: .  target - the value of the target

596:    Note:
597:    If the target was not set by the user, then zero is returned.

599:    Level: intermediate

601: .seealso: NEPSetTarget()
602: @*/
603: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
604: {
608:   *target = nep->target;
609:   return(0);
610: }

612: /*@C
613:    NEPSetFunction - Sets the function to compute the nonlinear Function T(lambda)
614:    as well as the location to store the matrix.

616:    Logically Collective on NEP and Mat

618:    Input Parameters:
619: +  nep - the NEP context
620: .  A   - Function matrix
621: .  B   - preconditioner matrix (usually same as the Function)
622: .  fun - Function evaluation routine (if NULL then NEP retains any
623:          previously set value)
624: -  ctx - [optional] user-defined context for private data for the Function
625:          evaluation routine (may be NULL) (if NULL then NEP retains any
626:          previously set value)

628:    Calling Sequence of fun:
629: $   fun(NEP nep,PetscScalar lambda,Mat F,Mat P,void *ctx)

631: +  nep    - the NEP context
632: .  lambda - the scalar argument where T(.) must be evaluated
633: .  T      - matrix that will contain T(lambda)
634: .  P      - (optional) different matrix to build the preconditioner
635: -  ctx    - (optional) user-defined context, as set by NEPSetFunction()

637:    Level: beginner

639: .seealso: NEPGetFunction(), NEPSetJacobian()
640: @*/
641: PetscErrorCode NEPSetFunction(NEP nep,Mat A,Mat B,PetscErrorCode (*fun)(NEP,PetscScalar,Mat,Mat,void*),void *ctx)
642: {


652:   if (nep->state) { NEPReset(nep); }
653:   else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) { NEPReset_Problem(nep); }

655:   if (fun) nep->computefunction = fun;
656:   if (ctx) nep->functionctx     = ctx;
657:   if (A) {
658:     PetscObjectReference((PetscObject)A);
659:     MatDestroy(&nep->function);
660:     nep->function = A;
661:   }
662:   if (B) {
663:     PetscObjectReference((PetscObject)B);
664:     MatDestroy(&nep->function_pre);
665:     nep->function_pre = B;
666:   }
667:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
668:   nep->state = NEP_STATE_INITIAL;
669:   return(0);
670: }

672: /*@C
673:    NEPGetFunction - Returns the Function matrix and optionally the user
674:    provided context for evaluating the Function.

676:    Not Collective, but Mat object will be parallel if NEP object is

678:    Input Parameter:
679: .  nep - the nonlinear eigensolver context

681:    Output Parameters:
682: +  A   - location to stash Function matrix (or NULL)
683: .  B   - location to stash preconditioner matrix (or NULL)
684: .  fun - location to put Function function (or NULL)
685: -  ctx - location to stash Function context (or NULL)

687:    Level: advanced

689: .seealso: NEPSetFunction()
690: @*/
691: PetscErrorCode NEPGetFunction(NEP nep,Mat *A,Mat *B,PetscErrorCode (**fun)(NEP,PetscScalar,Mat,Mat,void*),void **ctx)
692: {
695:   NEPCheckCallback(nep,1);
696:   if (A)   *A   = nep->function;
697:   if (B)   *B   = nep->function_pre;
698:   if (fun) *fun = nep->computefunction;
699:   if (ctx) *ctx = nep->functionctx;
700:   return(0);
701: }

703: /*@C
704:    NEPSetJacobian - Sets the function to compute Jacobian T'(lambda) as well
705:    as the location to store the matrix.

707:    Logically Collective on NEP and Mat

709:    Input Parameters:
710: +  nep - the NEP context
711: .  A   - Jacobian matrix
712: .  jac - Jacobian evaluation routine (if NULL then NEP retains any
713:          previously set value)
714: -  ctx - [optional] user-defined context for private data for the Jacobian
715:          evaluation routine (may be NULL) (if NULL then NEP retains any
716:          previously set value)

718:    Calling Sequence of jac:
719: $   jac(NEP nep,PetscScalar lambda,Mat J,void *ctx)

721: +  nep    - the NEP context
722: .  lambda - the scalar argument where T'(.) must be evaluated
723: .  J      - matrix that will contain T'(lambda)
724: -  ctx    - (optional) user-defined context, as set by NEPSetJacobian()

726:    Level: beginner

728: .seealso: NEPSetFunction(), NEPGetJacobian()
729: @*/
730: PetscErrorCode NEPSetJacobian(NEP nep,Mat A,PetscErrorCode (*jac)(NEP,PetscScalar,Mat,void*),void *ctx)
731: {


739:   if (nep->state) { NEPReset(nep); }
740:   else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) { NEPReset_Problem(nep); }

742:   if (jac) nep->computejacobian = jac;
743:   if (ctx) nep->jacobianctx     = ctx;
744:   if (A) {
745:     PetscObjectReference((PetscObject)A);
746:     MatDestroy(&nep->jacobian);
747:     nep->jacobian = A;
748:   }
749:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
750:   nep->state = NEP_STATE_INITIAL;
751:   return(0);
752: }

754: /*@C
755:    NEPGetJacobian - Returns the Jacobian matrix and optionally the user
756:    provided routine and context for evaluating the Jacobian.

758:    Not Collective, but Mat object will be parallel if NEP object is

760:    Input Parameter:
761: .  nep - the nonlinear eigensolver context

763:    Output Parameters:
764: +  A   - location to stash Jacobian matrix (or NULL)
765: .  jac - location to put Jacobian function (or NULL)
766: -  ctx - location to stash Jacobian context (or NULL)

768:    Level: advanced

770: .seealso: NEPSetJacobian()
771: @*/
772: PetscErrorCode NEPGetJacobian(NEP nep,Mat *A,PetscErrorCode (**jac)(NEP,PetscScalar,Mat,void*),void **ctx)
773: {
776:   NEPCheckCallback(nep,1);
777:   if (A)   *A   = nep->jacobian;
778:   if (jac) *jac = nep->computejacobian;
779:   if (ctx) *ctx = nep->jacobianctx;
780:   return(0);
781: }

783: /*@
784:    NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
785:    in split form.

787:    Collective on NEP, Mat and FN

789:    Input Parameters:
790: +  nep - the nonlinear eigensolver context
791: .  n   - number of terms in the split form
792: .  A   - array of matrices
793: .  f   - array of functions
794: -  str - structure flag for matrices

796:    Notes:
797:    The nonlinear operator is written as T(lambda) = sum_i A_i*f_i(lambda),
798:    for i=1,...,n. The derivative T'(lambda) can be obtained using the
799:    derivatives of f_i.

801:    The structure flag provides information about A_i's nonzero pattern
802:    (see MatStructure enum). If all matrices have the same pattern, then
803:    use SAME_NONZERO_PATTERN. If the patterns are different but contained
804:    in the pattern of the first one, then use SUBSET_NONZERO_PATTERN.
805:    Otherwise use DIFFERENT_NONZERO_PATTERN.

807:    This function must be called before NEPSetUp(). If it is called again
808:    after NEPSetUp() then the NEP object is reset.

810:    Level: beginner

812: .seealso: NEPGetSplitOperatorTerm(), NEPGetSplitOperatorInfo()
813:  @*/
814: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt n,Mat A[],FN f[],MatStructure str)
815: {
816:   PetscInt       i;

822:   if (n <= 0) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %D",n);

828:   for (i=0;i<n;i++) {
831:     PetscObjectReference((PetscObject)A[i]);
832:     PetscObjectReference((PetscObject)f[i]);
833:   }

835:   if (nep->state) { NEPReset(nep); }
836:   else { NEPReset_Problem(nep); }

838:   /* allocate space and copy matrices and functions */
839:   PetscMalloc1(n,&nep->A);
840:   PetscLogObjectMemory((PetscObject)nep,n*sizeof(Mat));
841:   for (i=0;i<n;i++) nep->A[i] = A[i];
842:   PetscMalloc1(n,&nep->f);
843:   PetscLogObjectMemory((PetscObject)nep,n*sizeof(FN));
844:   for (i=0;i<n;i++) nep->f[i] = f[i];
845:   PetscCalloc1(n,&nep->nrma);
846:   PetscLogObjectMemory((PetscObject)nep,n*sizeof(PetscReal));
847:   nep->nt    = n;
848:   nep->mstr  = str;
849:   nep->fui   = NEP_USER_INTERFACE_SPLIT;
850:   nep->state = NEP_STATE_INITIAL;
851:   return(0);
852: }

854: /*@
855:    NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
856:    the nonlinear operator in split form.

858:    Not collective, though parallel Mats and FNs are returned if the NEP is parallel

860:    Input Parameter:
861: +  nep - the nonlinear eigensolver context
862: -  k   - the index of the requested term (starting in 0)

864:    Output Parameters:
865: +  A - the matrix of the requested term
866: -  f - the function of the requested term

868:    Level: intermediate

870: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorInfo()
871: @*/
872: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
873: {
876:   NEPCheckSplit(nep,1);
877:   if (k<0 || k>=nep->nt) SETERRQ1(PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %D",nep->nt-1);
878:   if (A) *A = nep->A[k];
879:   if (f) *f = nep->f[k];
880:   return(0);
881: }

883: /*@
884:    NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
885:    the nonlinear operator, as well as the structure flag for matrices.

887:    Not collective

889:    Input Parameter:
890: .  nep - the nonlinear eigensolver context

892:    Output Parameters:
893: +  n   - the number of terms passed in NEPSetSplitOperator()
894: -  str - the matrix structure flag passed in NEPSetSplitOperator()

896:    Level: intermediate

898: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorTerm()
899: @*/
900: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
901: {
904:   NEPCheckSplit(nep,1);
905:   if (n)   *n = nep->nt;
906:   if (str) *str = nep->mstr;
907:   return(0);
908: }

910: /*@C
911:    NEPSetDerivatives - Sets the function to compute the k-th derivative T^(k)(lambda)
912:    for any value of k (including 0), as well as the location to store the matrix.

914:    Logically Collective on NEP and Mat

916:    Input Parameters:
917: +  nep - the NEP context
918: .  A   - the matrix to store the computed derivative
919: .  der - routing to evaluate the k-th derivative (if NULL then NEP retains any
920:          previously set value)
921: -  ctx - [optional] user-defined context for private data for the derivatives
922:          evaluation routine (may be NULL) (if NULL then NEP retains any
923:          previously set value)

925:    Level: beginner

927: .seealso: NEPSetFunction(), NEPGetDerivatives()
928: @*/
929: PetscErrorCode NEPSetDerivatives(NEP nep,Mat A,PetscErrorCode (*der)(NEP,PetscScalar,PetscInt,Mat,void*),void *ctx)
930: {


938:   if (nep->state) { NEPReset(nep); }
939:   else { NEPReset_Problem(nep); }

941:   if (der) nep->computederivatives = der;
942:   if (ctx) nep->derivativesctx     = ctx;
943:   if (A) {
944:     PetscObjectReference((PetscObject)A);
945:     MatDestroy(&nep->derivatives);
946:     nep->derivatives = A;
947:   }
948:   nep->fui   = NEP_USER_INTERFACE_DERIVATIVES;
949:   nep->state = NEP_STATE_INITIAL;
950:   return(0);
951: }

953: /*@C
954:    NEPGetDerivatives - Returns the derivatives matrix and optionally the user
955:    provided routine and context for evaluating the derivatives.

957:    Not Collective, but Mat object will be parallel if NEP object is

959:    Input Parameter:
960: .  nep - the nonlinear eigensolver context

962:    Output Parameters:
963: +  A   - location to stash the derivatives matrix (or NULL)
964: .  der - location to put derivatives function (or NULL)
965: -  ctx - location to stash derivatives context (or NULL)

967:    Level: advanced

969: .seealso: NEPSetDerivatives()
970: @*/
971: PetscErrorCode NEPGetDerivatives(NEP nep,Mat *A,PetscErrorCode (**der)(NEP,PetscScalar,PetscInt,Mat,void*),void **ctx)
972: {
975:   NEPCheckDerivatives(nep,1);
976:   if (A)   *A   = nep->derivatives;
977:   if (der) *der = nep->computederivatives;
978:   if (ctx) *ctx = nep->derivativesctx;
979:   return(0);
980: }