Actual source code: nepbasic.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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>

 16: /* Logging support */
 17: PetscClassId      NEP_CLASSID = 0;
 18: PetscLogEvent     NEP_SetUp = 0,NEP_Solve = 0,NEP_Refine = 0,NEP_FunctionEval = 0,NEP_JacobianEval = 0,NEP_Resolvent = 0,NEP_CISS_SVD = 0;

 20: /* List of registered NEP routines */
 21: PetscFunctionList NEPList = 0;
 22: PetscBool         NEPRegisterAllCalled = PETSC_FALSE;

 24: /* List of registered NEP monitors */
 25: PetscFunctionList NEPMonitorList              = NULL;
 26: PetscFunctionList NEPMonitorCreateList        = NULL;
 27: PetscFunctionList NEPMonitorDestroyList       = NULL;
 28: PetscBool         NEPMonitorRegisterAllCalled = PETSC_FALSE;

 30: /*@
 31:    NEPCreate - Creates the default NEP context.

 33:    Collective

 35:    Input Parameter:
 36: .  comm - MPI communicator

 38:    Output Parameter:
 39: .  outnep - location to put the NEP context

 41:    Level: beginner

 43: .seealso: NEPSetUp(), NEPSolve(), NEPDestroy(), NEP
 44: @*/
 45: PetscErrorCode NEPCreate(MPI_Comm comm,NEP *outnep)
 46: {
 47:   NEP            nep;

 50:   *outnep = 0;
 51:   NEPInitializePackage();
 52:   SlepcHeaderCreate(nep,NEP_CLASSID,"NEP","Nonlinear Eigenvalue Problem","NEP",comm,NEPDestroy,NEPView);

 54:   nep->max_it          = PETSC_DEFAULT;
 55:   nep->nev             = 1;
 56:   nep->ncv             = PETSC_DEFAULT;
 57:   nep->mpd             = PETSC_DEFAULT;
 58:   nep->nini            = 0;
 59:   nep->target          = 0.0;
 60:   nep->tol             = PETSC_DEFAULT;
 61:   nep->conv            = NEP_CONV_REL;
 62:   nep->stop            = NEP_STOP_BASIC;
 63:   nep->which           = (NEPWhich)0;
 64:   nep->problem_type    = (NEPProblemType)0;
 65:   nep->refine          = NEP_REFINE_NONE;
 66:   nep->npart           = 1;
 67:   nep->rtol            = PETSC_DEFAULT;
 68:   nep->rits            = PETSC_DEFAULT;
 69:   nep->scheme          = (NEPRefineScheme)0;
 70:   nep->trackall        = PETSC_FALSE;
 71:   nep->twosided        = PETSC_FALSE;

 73:   nep->computefunction = NULL;
 74:   nep->computejacobian = NULL;
 75:   nep->functionctx     = NULL;
 76:   nep->jacobianctx     = NULL;
 77:   nep->converged       = NEPConvergedRelative;
 78:   nep->convergeduser   = NULL;
 79:   nep->convergeddestroy= NULL;
 80:   nep->stopping        = NEPStoppingBasic;
 81:   nep->stoppinguser    = NULL;
 82:   nep->stoppingdestroy = NULL;
 83:   nep->convergedctx    = NULL;
 84:   nep->stoppingctx     = NULL;
 85:   nep->numbermonitors  = 0;

 87:   nep->ds              = NULL;
 88:   nep->V               = NULL;
 89:   nep->W               = NULL;
 90:   nep->rg              = NULL;
 91:   nep->function        = NULL;
 92:   nep->function_pre    = NULL;
 93:   nep->jacobian        = NULL;
 94:   nep->A               = NULL;
 95:   nep->f               = NULL;
 96:   nep->nt              = 0;
 97:   nep->mstr            = UNKNOWN_NONZERO_PATTERN;
 98:   nep->P               = NULL;
 99:   nep->mstrp           = UNKNOWN_NONZERO_PATTERN;
100:   nep->IS              = NULL;
101:   nep->eigr            = NULL;
102:   nep->eigi            = NULL;
103:   nep->errest          = NULL;
104:   nep->perm            = NULL;
105:   nep->nwork           = 0;
106:   nep->work            = NULL;
107:   nep->data            = NULL;

109:   nep->state           = NEP_STATE_INITIAL;
110:   nep->nconv           = 0;
111:   nep->its             = 0;
112:   nep->n               = 0;
113:   nep->nloc            = 0;
114:   nep->nrma            = NULL;
115:   nep->fui             = (NEPUserInterface)0;
116:   nep->useds           = PETSC_FALSE;
117:   nep->resolvent       = NULL;
118:   nep->reason          = NEP_CONVERGED_ITERATING;

120:   PetscNewLog(nep,&nep->sc);
121:   *outnep = nep;
122:   PetscFunctionReturn(0);
123: }

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

128:    Logically Collective on nep

130:    Input Parameters:
131: +  nep      - the nonlinear eigensolver context
132: -  type     - a known method

134:    Options Database Key:
135: .  -nep_type <method> - Sets the method; use -help for a list
136:     of available methods

138:    Notes:
139:    See "slepc/include/slepcnep.h" for available methods.

141:    Normally, it is best to use the NEPSetFromOptions() command and
142:    then set the NEP type from the options database rather than by using
143:    this routine.  Using the options database provides the user with
144:    maximum flexibility in evaluating the different available methods.
145:    The NEPSetType() routine is provided for those situations where it
146:    is necessary to set the iterative solver independently of the command
147:    line or options database.

149:    Level: intermediate

151: .seealso: NEPType
152: @*/
153: PetscErrorCode NEPSetType(NEP nep,NEPType type)
154: {
155:   PetscErrorCode (*r)(NEP);
156:   PetscBool      match;


161:   PetscObjectTypeCompare((PetscObject)nep,type,&match);
162:   if (match) PetscFunctionReturn(0);

164:   PetscFunctionListFind(NEPList,type,&r);

167:   if (nep->ops->destroy) (*nep->ops->destroy)(nep);
168:   PetscMemzero(nep->ops,sizeof(struct _NEPOps));

170:   nep->state = NEP_STATE_INITIAL;
171:   PetscObjectChangeTypeName((PetscObject)nep,type);
172:   (*r)(nep);
173:   PetscFunctionReturn(0);
174: }

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

179:    Not Collective

181:    Input Parameter:
182: .  nep - the eigensolver context

184:    Output Parameter:
185: .  type - name of NEP method

187:    Level: intermediate

189: .seealso: NEPSetType()
190: @*/
191: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
192: {
195:   *type = ((PetscObject)nep)->type_name;
196:   PetscFunctionReturn(0);
197: }

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

202:    Not Collective

204:    Input Parameters:
205: +  name - name of a new user-defined solver
206: -  function - routine to create the solver context

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

211:    Sample usage:
212: .vb
213:     NEPRegister("my_solver",MySolverCreate);
214: .ve

216:    Then, your solver can be chosen with the procedural interface via
217: $     NEPSetType(nep,"my_solver")
218:    or at runtime via the option
219: $     -nep_type my_solver

221:    Level: advanced

223: .seealso: NEPRegisterAll()
224: @*/
225: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
226: {
227:   NEPInitializePackage();
228:   PetscFunctionListAdd(&NEPList,name,function);
229:   PetscFunctionReturn(0);
230: }

232: /*@C
233:    NEPMonitorRegister - Adds NEP monitor routine.

235:    Not Collective

237:    Input Parameters:
238: +  name    - name of a new monitor routine
239: .  vtype   - a PetscViewerType for the output
240: .  format  - a PetscViewerFormat for the output
241: .  monitor - monitor routine
242: .  create  - creation routine, or NULL
243: -  destroy - destruction routine, or NULL

245:    Notes:
246:    NEPMonitorRegister() may be called multiple times to add several user-defined monitors.

248:    Sample usage:
249: .vb
250:    NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
251: .ve

253:    Then, your monitor can be chosen with the procedural interface via
254: $      NEPMonitorSetFromOptions(nep,"-nep_monitor_my_monitor","my_monitor",NULL)
255:    or at runtime via the option
256: $      -nep_monitor_my_monitor

258:    Level: advanced

260: .seealso: NEPMonitorRegisterAll()
261: @*/
262: PetscErrorCode NEPMonitorRegister(const char name[],PetscViewerType vtype,PetscViewerFormat format,PetscErrorCode (*monitor)(NEP,PetscInt,PetscInt,PetscScalar*,PetscScalar*,PetscReal*,PetscInt,PetscViewerAndFormat*),PetscErrorCode (*create)(PetscViewer,PetscViewerFormat,void*,PetscViewerAndFormat**),PetscErrorCode (*destroy)(PetscViewerAndFormat**))
263: {
264:   char           key[PETSC_MAX_PATH_LEN];

266:   NEPInitializePackage();
267:   SlepcMonitorMakeKey_Internal(name,vtype,format,key);
268:   PetscFunctionListAdd(&NEPMonitorList,key,monitor);
269:   if (create)  PetscFunctionListAdd(&NEPMonitorCreateList,key,create);
270:   if (destroy) PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy);
271:   PetscFunctionReturn(0);
272: }

274: /*
275:    NEPReset_Problem - Destroys the problem matrices.
276: */
277: PetscErrorCode NEPReset_Problem(NEP nep)
278: {
279:   PetscInt       i;

282:   MatDestroy(&nep->function);
283:   MatDestroy(&nep->function_pre);
284:   MatDestroy(&nep->jacobian);
285:   if (nep->fui==NEP_USER_INTERFACE_SPLIT) {
286:     MatDestroyMatrices(nep->nt,&nep->A);
287:     for (i=0;i<nep->nt;i++) FNDestroy(&nep->f[i]);
288:     PetscFree(nep->f);
289:     PetscFree(nep->nrma);
290:     if (nep->P) MatDestroyMatrices(nep->nt,&nep->P);
291:     nep->nt = 0;
292:   }
293:   PetscFunctionReturn(0);
294: }
295: /*@
296:    NEPReset - Resets the NEP context to the initial state (prior to setup)
297:    and destroys any allocated Vecs and Mats.

299:    Collective on nep

301:    Input Parameter:
302: .  nep - eigensolver context obtained from NEPCreate()

304:    Level: advanced

306: .seealso: NEPDestroy()
307: @*/
308: PetscErrorCode NEPReset(NEP nep)
309: {
311:   if (!nep) PetscFunctionReturn(0);
312:   if (nep->ops->reset) (nep->ops->reset)(nep);
313:   if (nep->refineksp) KSPReset(nep->refineksp);
314:   NEPReset_Problem(nep);
315:   BVDestroy(&nep->V);
316:   BVDestroy(&nep->W);
317:   VecDestroyVecs(nep->nwork,&nep->work);
318:   MatDestroy(&nep->resolvent);
319:   nep->nwork = 0;
320:   nep->state = NEP_STATE_INITIAL;
321:   PetscFunctionReturn(0);
322: }

324: /*@C
325:    NEPDestroy - Destroys the NEP context.

327:    Collective on nep

329:    Input Parameter:
330: .  nep - eigensolver context obtained from NEPCreate()

332:    Level: beginner

334: .seealso: NEPCreate(), NEPSetUp(), NEPSolve()
335: @*/
336: PetscErrorCode NEPDestroy(NEP *nep)
337: {
338:   if (!*nep) PetscFunctionReturn(0);
340:   if (--((PetscObject)(*nep))->refct > 0) { *nep = 0; PetscFunctionReturn(0); }
341:   NEPReset(*nep);
342:   if ((*nep)->ops->destroy) (*(*nep)->ops->destroy)(*nep);
343:   if ((*nep)->eigr) PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm);
344:   RGDestroy(&(*nep)->rg);
345:   DSDestroy(&(*nep)->ds);
346:   KSPDestroy(&(*nep)->refineksp);
347:   PetscSubcommDestroy(&(*nep)->refinesubc);
348:   PetscFree((*nep)->sc);
349:   /* just in case the initial vectors have not been used */
350:   SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS);
351:   if ((*nep)->convergeddestroy) (*(*nep)->convergeddestroy)((*nep)->convergedctx);
352:   NEPMonitorCancel(*nep);
353:   PetscHeaderDestroy(nep);
354:   PetscFunctionReturn(0);
355: }

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

360:    Collective on nep

362:    Input Parameters:
363: +  nep - eigensolver context obtained from NEPCreate()
364: -  bv  - the basis vectors object

366:    Note:
367:    Use NEPGetBV() to retrieve the basis vectors context (for example,
368:    to free it at the end of the computations).

370:    Level: advanced

372: .seealso: NEPGetBV()
373: @*/
374: PetscErrorCode NEPSetBV(NEP nep,BV bv)
375: {
379:   PetscObjectReference((PetscObject)bv);
380:   BVDestroy(&nep->V);
381:   nep->V = bv;
382:   PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->V);
383:   PetscFunctionReturn(0);
384: }

386: /*@
387:    NEPGetBV - Obtain the basis vectors object associated to the nonlinear
388:    eigensolver object.

390:    Not Collective

392:    Input Parameters:
393: .  nep - eigensolver context obtained from NEPCreate()

395:    Output Parameter:
396: .  bv - basis vectors context

398:    Level: advanced

400: .seealso: NEPSetBV()
401: @*/
402: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
403: {
406:   if (!nep->V) {
407:     BVCreate(PetscObjectComm((PetscObject)nep),&nep->V);
408:     PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0);
409:     PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->V);
410:     PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options);
411:   }
412:   *bv = nep->V;
413:   PetscFunctionReturn(0);
414: }

416: /*@
417:    NEPSetRG - Associates a region object to the nonlinear eigensolver.

419:    Collective on nep

421:    Input Parameters:
422: +  nep - eigensolver context obtained from NEPCreate()
423: -  rg  - the region object

425:    Note:
426:    Use NEPGetRG() to retrieve the region context (for example,
427:    to free it at the end of the computations).

429:    Level: advanced

431: .seealso: NEPGetRG()
432: @*/
433: PetscErrorCode NEPSetRG(NEP nep,RG rg)
434: {
436:   if (rg) {
439:   }
440:   PetscObjectReference((PetscObject)rg);
441:   RGDestroy(&nep->rg);
442:   nep->rg = rg;
443:   PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->rg);
444:   PetscFunctionReturn(0);
445: }

447: /*@
448:    NEPGetRG - Obtain the region object associated to the
449:    nonlinear eigensolver object.

451:    Not Collective

453:    Input Parameters:
454: .  nep - eigensolver context obtained from NEPCreate()

456:    Output Parameter:
457: .  rg - region context

459:    Level: advanced

461: .seealso: NEPSetRG()
462: @*/
463: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
464: {
467:   if (!nep->rg) {
468:     RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg);
469:     PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0);
470:     PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->rg);
471:     PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options);
472:   }
473:   *rg = nep->rg;
474:   PetscFunctionReturn(0);
475: }

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

480:    Collective on nep

482:    Input Parameters:
483: +  nep - eigensolver context obtained from NEPCreate()
484: -  ds  - the direct solver object

486:    Note:
487:    Use NEPGetDS() to retrieve the direct solver context (for example,
488:    to free it at the end of the computations).

490:    Level: advanced

492: .seealso: NEPGetDS()
493: @*/
494: PetscErrorCode NEPSetDS(NEP nep,DS ds)
495: {
499:   PetscObjectReference((PetscObject)ds);
500:   DSDestroy(&nep->ds);
501:   nep->ds = ds;
502:   PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->ds);
503:   PetscFunctionReturn(0);
504: }

506: /*@
507:    NEPGetDS - Obtain the direct solver object associated to the
508:    nonlinear eigensolver object.

510:    Not Collective

512:    Input Parameters:
513: .  nep - eigensolver context obtained from NEPCreate()

515:    Output Parameter:
516: .  ds - direct solver context

518:    Level: advanced

520: .seealso: NEPSetDS()
521: @*/
522: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
523: {
526:   if (!nep->ds) {
527:     DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds);
528:     PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0);
529:     PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->ds);
530:     PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options);
531:   }
532:   *ds = nep->ds;
533:   PetscFunctionReturn(0);
534: }

536: /*@
537:    NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
538:    object in the refinement phase.

540:    Not Collective

542:    Input Parameters:
543: .  nep - eigensolver context obtained from NEPCreate()

545:    Output Parameter:
546: .  ksp - ksp context

548:    Level: advanced

550: .seealso: NEPSetRefine()
551: @*/
552: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
553: {
554:   MPI_Comm       comm;

558:   if (!nep->refineksp) {
559:     if (nep->npart>1) {
560:       /* Split in subcomunicators */
561:       PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc);
562:       PetscSubcommSetNumber(nep->refinesubc,nep->npart);
563:       PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS);
564:       PetscLogObjectMemory((PetscObject)nep,sizeof(PetscSubcomm));
565:       PetscSubcommGetChild(nep->refinesubc,&comm);
566:     } else PetscObjectGetComm((PetscObject)nep,&comm);
567:     KSPCreate(comm,&nep->refineksp);
568:     PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0);
569:     PetscLogObjectParent((PetscObject)nep,(PetscObject)nep->refineksp);
570:     PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options);
571:     KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix);
572:     KSPAppendOptionsPrefix(*ksp,"nep_refine_");
573:     KSPSetTolerances(nep->refineksp,SlepcDefaultTol(nep->rtol),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
574:   }
575:   *ksp = nep->refineksp;
576:   PetscFunctionReturn(0);
577: }

579: /*@
580:    NEPSetTarget - Sets the value of the target.

582:    Logically Collective on nep

584:    Input Parameters:
585: +  nep    - eigensolver context
586: -  target - the value of the target

588:    Options Database Key:
589: .  -nep_target <scalar> - the value of the target

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

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

599:    Level: intermediate

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

611: /*@
612:    NEPGetTarget - Gets the value of the target.

614:    Not Collective

616:    Input Parameter:
617: .  nep - eigensolver context

619:    Output Parameter:
620: .  target - the value of the target

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

625:    Level: intermediate

627: .seealso: NEPSetTarget()
628: @*/
629: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
630: {
633:   *target = nep->target;
634:   PetscFunctionReturn(0);
635: }

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

641:    Logically Collective on nep

643:    Input Parameters:
644: +  nep - the NEP context
645: .  A   - Function matrix
646: .  B   - preconditioner matrix (usually same as A)
647: .  fun - Function evaluation routine (if NULL then NEP retains any
648:          previously set value)
649: -  ctx - [optional] user-defined context for private data for the Function
650:          evaluation routine (may be NULL) (if NULL then NEP retains any
651:          previously set value)

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

656: +  nep    - the NEP context
657: .  lambda - the scalar argument where T(.) must be evaluated
658: .  T      - matrix that will contain T(lambda)
659: .  P      - (optional) different matrix to build the preconditioner
660: -  ctx    - (optional) user-defined context, as set by NEPSetFunction()

662:    Level: beginner

664: .seealso: NEPGetFunction(), NEPSetJacobian()
665: @*/
666: PetscErrorCode NEPSetFunction(NEP nep,Mat A,Mat B,PetscErrorCode (*fun)(NEP,PetscScalar,Mat,Mat,void*),void *ctx)
667: {

674:   if (nep->state) NEPReset(nep);
675:   else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) NEPReset_Problem(nep);

677:   if (fun) nep->computefunction = fun;
678:   if (ctx) nep->functionctx     = ctx;
679:   if (A) {
680:     PetscObjectReference((PetscObject)A);
681:     MatDestroy(&nep->function);
682:     nep->function = A;
683:   }
684:   if (B) {
685:     PetscObjectReference((PetscObject)B);
686:     MatDestroy(&nep->function_pre);
687:     nep->function_pre = B;
688:   }
689:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
690:   nep->state = NEP_STATE_INITIAL;
691:   PetscFunctionReturn(0);
692: }

694: /*@C
695:    NEPGetFunction - Returns the Function matrix and optionally the user
696:    provided context for evaluating the Function.

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

700:    Input Parameter:
701: .  nep - the nonlinear eigensolver context

703:    Output Parameters:
704: +  A   - location to stash Function matrix (or NULL)
705: .  B   - location to stash preconditioner matrix (or NULL)
706: .  fun - location to put Function function (or NULL)
707: -  ctx - location to stash Function context (or NULL)

709:    Level: advanced

711: .seealso: NEPSetFunction()
712: @*/
713: PetscErrorCode NEPGetFunction(NEP nep,Mat *A,Mat *B,PetscErrorCode (**fun)(NEP,PetscScalar,Mat,Mat,void*),void **ctx)
714: {
716:   NEPCheckCallback(nep,1);
717:   if (A)   *A   = nep->function;
718:   if (B)   *B   = nep->function_pre;
719:   if (fun) *fun = nep->computefunction;
720:   if (ctx) *ctx = nep->functionctx;
721:   PetscFunctionReturn(0);
722: }

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

728:    Logically Collective on nep

730:    Input Parameters:
731: +  nep - the NEP context
732: .  A   - Jacobian matrix
733: .  jac - Jacobian evaluation routine (if NULL then NEP retains any
734:          previously set value)
735: -  ctx - [optional] user-defined context for private data for the Jacobian
736:          evaluation routine (may be NULL) (if NULL then NEP retains any
737:          previously set value)

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

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

747:    Level: beginner

749: .seealso: NEPSetFunction(), NEPGetJacobian()
750: @*/
751: PetscErrorCode NEPSetJacobian(NEP nep,Mat A,PetscErrorCode (*jac)(NEP,PetscScalar,Mat,void*),void *ctx)
752: {

757:   if (nep->state) NEPReset(nep);
758:   else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) NEPReset_Problem(nep);

760:   if (jac) nep->computejacobian = jac;
761:   if (ctx) nep->jacobianctx     = ctx;
762:   if (A) {
763:     PetscObjectReference((PetscObject)A);
764:     MatDestroy(&nep->jacobian);
765:     nep->jacobian = A;
766:   }
767:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
768:   nep->state = NEP_STATE_INITIAL;
769:   PetscFunctionReturn(0);
770: }

772: /*@C
773:    NEPGetJacobian - Returns the Jacobian matrix and optionally the user
774:    provided routine and context for evaluating the Jacobian.

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

778:    Input Parameter:
779: .  nep - the nonlinear eigensolver context

781:    Output Parameters:
782: +  A   - location to stash Jacobian matrix (or NULL)
783: .  jac - location to put Jacobian function (or NULL)
784: -  ctx - location to stash Jacobian context (or NULL)

786:    Level: advanced

788: .seealso: NEPSetJacobian()
789: @*/
790: PetscErrorCode NEPGetJacobian(NEP nep,Mat *A,PetscErrorCode (**jac)(NEP,PetscScalar,Mat,void*),void **ctx)
791: {
793:   NEPCheckCallback(nep,1);
794:   if (A)   *A   = nep->jacobian;
795:   if (jac) *jac = nep->computejacobian;
796:   if (ctx) *ctx = nep->jacobianctx;
797:   PetscFunctionReturn(0);
798: }

800: /*@
801:    NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
802:    in split form.

804:    Collective on nep

806:    Input Parameters:
807: +  nep - the nonlinear eigensolver context
808: .  nt  - number of terms in the split form
809: .  A   - array of matrices
810: .  f   - array of functions
811: -  str - structure flag for matrices

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

818:    The structure flag provides information about A_i's nonzero pattern
819:    (see MatStructure enum). If all matrices have the same pattern, then
820:    use SAME_NONZERO_PATTERN. If the patterns are different but contained
821:    in the pattern of the first one, then use SUBSET_NONZERO_PATTERN. If
822:    patterns are known to be different, use DIFFERENT_NONZERO_PATTERN.
823:    If set to UNKNOWN_NONZERO_PATTERN, the patterns will be compared to
824:    determine if they are equal.

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

829:    Level: beginner

831: .seealso: NEPGetSplitOperatorTerm(), NEPGetSplitOperatorInfo(), NEPSetSplitPreconditioner()
832:  @*/
833: PetscErrorCode NEPSetSplitOperator(NEP nep,PetscInt nt,Mat A[],FN f[],MatStructure str)
834: {
835:   PetscInt       i,n=0,m,m0=0,mloc,nloc,mloc0=0;


844:   for (i=0;i<nt;i++) {
849:     MatGetSize(A[i],&m,&n);
850:     MatGetLocalSize(A[i],&mloc,&nloc);
853:     if (!i) { m0 = m; mloc0 = mloc; }
856:     PetscObjectReference((PetscObject)A[i]);
857:     PetscObjectReference((PetscObject)f[i]);
858:   }

860:   if (nep->state && (n!=nep->n || nloc!=nep->nloc)) NEPReset(nep);
861:   else NEPReset_Problem(nep);

863:   /* allocate space and copy matrices and functions */
864:   PetscMalloc1(nt,&nep->A);
865:   PetscLogObjectMemory((PetscObject)nep,nt*sizeof(Mat));
866:   for (i=0;i<nt;i++) nep->A[i] = A[i];
867:   PetscMalloc1(nt,&nep->f);
868:   PetscLogObjectMemory((PetscObject)nep,nt*sizeof(FN));
869:   for (i=0;i<nt;i++) nep->f[i] = f[i];
870:   PetscCalloc1(nt,&nep->nrma);
871:   PetscLogObjectMemory((PetscObject)nep,nt*sizeof(PetscReal));
872:   nep->nt    = nt;
873:   nep->mstr  = str;
874:   nep->fui   = NEP_USER_INTERFACE_SPLIT;
875:   nep->state = NEP_STATE_INITIAL;
876:   PetscFunctionReturn(0);
877: }

879: /*@
880:    NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
881:    the nonlinear operator in split form.

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

885:    Input Parameters:
886: +  nep - the nonlinear eigensolver context
887: -  k   - the index of the requested term (starting in 0)

889:    Output Parameters:
890: +  A - the matrix of the requested term
891: -  f - the function of the requested term

893:    Level: intermediate

895: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorInfo()
896: @*/
897: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
898: {
901:   NEPCheckSplit(nep,1);
903:   if (A) *A = nep->A[k];
904:   if (f) *f = nep->f[k];
905:   PetscFunctionReturn(0);
906: }

908: /*@
909:    NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
910:    the nonlinear operator, as well as the structure flag for matrices.

912:    Not collective

914:    Input Parameter:
915: .  nep - the nonlinear eigensolver context

917:    Output Parameters:
918: +  n   - the number of terms passed in NEPSetSplitOperator()
919: -  str - the matrix structure flag passed in NEPSetSplitOperator()

921:    Level: intermediate

923: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorTerm()
924: @*/
925: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
926: {
928:   NEPCheckSplit(nep,1);
929:   if (n)   *n = nep->nt;
930:   if (str) *str = nep->mstr;
931:   PetscFunctionReturn(0);
932: }

934: /*@
935:    NEPSetSplitPreconditioner - Sets an operator in split form from which
936:    to build the preconditioner to be used when solving the nonlinear
937:    eigenvalue problem in split form.

939:    Collective on nep

941:    Input Parameters:
942: +  nep  - the nonlinear eigensolver context
943: .  ntp  - number of terms in the split preconditioner
944: .  P    - array of matrices
945: -  strp - structure flag for matrices

947:    Notes:
948:    The matrix for the preconditioner is expressed as P(lambda) =
949:    sum_i P_i*f_i(lambda), for i=1,...,n, where the f_i functions
950:    are the same as in NEPSetSplitOperator(). It is not necessary to call
951:    this function. If it is not invoked, then the preconditioner is
952:    built from T(lambda), i.e., both matrices and functions passed in
953:    NEPSetSplitOperator().

955:    The structure flag provides information about P_i's nonzero pattern
956:    in the same way as in NEPSetSplitOperator().

958:    If the functions defining the preconditioner operator were different
959:    from the ones given in NEPSetSplitOperator(), then the split form
960:    cannot be used. Use the callback interface instead.

962:    Use ntp=0 to reset a previously set split preconditioner.

964:    Level: advanced

966: .seealso: NEPGetSplitPreconditionerTerm(), NEPGetSplitPreconditionerInfo(), NEPSetSplitOperator()
967:  @*/
968: PetscErrorCode NEPSetSplitPreconditioner(NEP nep,PetscInt ntp,Mat P[],MatStructure strp)
969: {
970:   PetscInt       i,n=0,m,m0=0,mloc,nloc,mloc0=0;


980:   for (i=0;i<ntp;i++) {
983:     MatGetSize(P[i],&m,&n);
984:     MatGetLocalSize(P[i],&mloc,&nloc);
987:     if (!i) { m0 = m; mloc0 = mloc; }
990:     PetscObjectReference((PetscObject)P[i]);
991:   }

994:   if (nep->P) MatDestroyMatrices(nep->nt,&nep->P);

996:   /* allocate space and copy matrices */
997:   if (ntp) {
998:     PetscMalloc1(ntp,&nep->P);
999:     PetscLogObjectMemory((PetscObject)nep,ntp*sizeof(Mat));
1000:     for (i=0;i<ntp;i++) nep->P[i] = P[i];
1001:   }
1002:   nep->mstrp = strp;
1003:   nep->state = NEP_STATE_INITIAL;
1004:   PetscFunctionReturn(0);
1005: }

1007: /*@
1008:    NEPGetSplitPreconditionerTerm - Gets the matrices associated with
1009:    the split preconditioner.

1011:    Not collective, though parallel Mats are returned if the NEP is parallel

1013:    Input Parameters:
1014: +  nep - the nonlinear eigensolver context
1015: -  k   - the index of the requested term (starting in 0)

1017:    Output Parameter:
1018: .  P  - the matrix of the requested term

1020:    Level: advanced

1022: .seealso: NEPSetSplitPreconditioner(), NEPGetSplitPreconditionerInfo()
1023: @*/
1024: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
1025: {
1029:   NEPCheckSplit(nep,1);
1032:   *P = nep->P[k];
1033:   PetscFunctionReturn(0);
1034: }

1036: /*@
1037:    NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
1038:    preconditioner, as well as the structure flag for matrices.

1040:    Not collective

1042:    Input Parameter:
1043: .  nep - the nonlinear eigensolver context

1045:    Output Parameters:
1046: +  n    - the number of terms passed in NEPSetSplitPreconditioner()
1047: -  strp - the matrix structure flag passed in NEPSetSplitPreconditioner()

1049:    Level: advanced

1051: .seealso: NEPSetSplitPreconditioner(), NEPGetSplitPreconditionerTerm()
1052: @*/
1053: PetscErrorCode NEPGetSplitPreconditionerInfo(NEP nep,PetscInt *n,MatStructure *strp)
1054: {
1056:   NEPCheckSplit(nep,1);
1057:   if (n)    *n    = nep->P? nep->nt: 0;
1058:   if (strp) *strp = nep->mstrp;
1059:   PetscFunctionReturn(0);
1060: }