Actual source code: nepbasic.c

slepc-main 2024-11-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 = NULL;
 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;

 49:   PetscFunctionBegin;
 50:   PetscAssertPointer(outnep,2);
 51:   PetscCall(NEPInitializePackage());
 52:   PetscCall(SlepcHeaderCreate(nep,NEP_CLASSID,"NEP","Nonlinear Eigenvalue Problem","NEP",comm,NEPDestroy,NEPView));

 54:   nep->max_it          = PETSC_DETERMINE;
 55:   nep->nev             = 1;
 56:   nep->ncv             = PETSC_DETERMINE;
 57:   nep->mpd             = PETSC_DETERMINE;
 58:   nep->nini            = 0;
 59:   nep->target          = 0.0;
 60:   nep->tol             = PETSC_DETERMINE;
 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_DETERMINE;
 68:   nep->rits            = PETSC_DETERMINE;
 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:   PetscCall(PetscNew(&nep->sc));
121:   *outnep = nep;
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

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

128:    Logically Collective

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;

158:   PetscFunctionBegin;
160:   PetscAssertPointer(type,2);

162:   PetscCall(PetscObjectTypeCompare((PetscObject)nep,type,&match));
163:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

165:   PetscCall(PetscFunctionListFind(NEPList,type,&r));
166:   PetscCheck(r,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown NEP type given: %s",type);

168:   PetscTryTypeMethod(nep,destroy);
169:   PetscCall(PetscMemzero(nep->ops,sizeof(struct _NEPOps)));

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

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

180:    Not Collective

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

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

188:    Level: intermediate

190: .seealso: NEPSetType()
191: @*/
192: PetscErrorCode NEPGetType(NEP nep,NEPType *type)
193: {
194:   PetscFunctionBegin;
196:   PetscAssertPointer(type,2);
197:   *type = ((PetscObject)nep)->type_name;
198:   PetscFunctionReturn(PETSC_SUCCESS);
199: }

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

204:    Not Collective

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

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

213:    Example Usage:
214: .vb
215:     NEPRegister("my_solver",MySolverCreate);
216: .ve

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

223:    Level: advanced

225: .seealso: NEPRegisterAll()
226: @*/
227: PetscErrorCode NEPRegister(const char *name,PetscErrorCode (*function)(NEP))
228: {
229:   PetscFunctionBegin;
230:   PetscCall(NEPInitializePackage());
231:   PetscCall(PetscFunctionListAdd(&NEPList,name,function));
232:   PetscFunctionReturn(PETSC_SUCCESS);
233: }

235: /*@C
236:    NEPMonitorRegister - Adds NEP monitor routine.

238:    Not Collective

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

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

251:    Example Usage:
252: .vb
253:    NEPMonitorRegister("my_monitor",PETSCVIEWERASCII,PETSC_VIEWER_ASCII_INFO_DETAIL,MyMonitor,NULL,NULL);
254: .ve

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

261:    Level: advanced

263: .seealso: NEPMonitorRegisterAll()
264: @*/
265: 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**))
266: {
267:   char           key[PETSC_MAX_PATH_LEN];

269:   PetscFunctionBegin;
270:   PetscCall(NEPInitializePackage());
271:   PetscCall(SlepcMonitorMakeKey_Internal(name,vtype,format,key));
272:   PetscCall(PetscFunctionListAdd(&NEPMonitorList,key,monitor));
273:   if (create)  PetscCall(PetscFunctionListAdd(&NEPMonitorCreateList,key,create));
274:   if (destroy) PetscCall(PetscFunctionListAdd(&NEPMonitorDestroyList,key,destroy));
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

278: /*
279:    NEPReset_Problem - Destroys the problem matrices.
280: */
281: PetscErrorCode NEPReset_Problem(NEP nep)
282: {
283:   PetscInt       i;

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

304:    Collective

306:    Input Parameter:
307: .  nep - eigensolver context obtained from NEPCreate()

309:    Level: advanced

311: .seealso: NEPDestroy()
312: @*/
313: PetscErrorCode NEPReset(NEP nep)
314: {
315:   PetscFunctionBegin;
317:   if (!nep) PetscFunctionReturn(PETSC_SUCCESS);
318:   PetscTryTypeMethod(nep,reset);
319:   if (nep->refineksp) PetscCall(KSPReset(nep->refineksp));
320:   PetscCall(NEPReset_Problem(nep));
321:   PetscCall(BVDestroy(&nep->V));
322:   PetscCall(BVDestroy(&nep->W));
323:   PetscCall(VecDestroyVecs(nep->nwork,&nep->work));
324:   PetscCall(MatDestroy(&nep->resolvent));
325:   nep->nwork = 0;
326:   nep->state = NEP_STATE_INITIAL;
327:   PetscFunctionReturn(PETSC_SUCCESS);
328: }

330: /*@
331:    NEPDestroy - Destroys the NEP context.

333:    Collective

335:    Input Parameter:
336: .  nep - eigensolver context obtained from NEPCreate()

338:    Level: beginner

340: .seealso: NEPCreate(), NEPSetUp(), NEPSolve()
341: @*/
342: PetscErrorCode NEPDestroy(NEP *nep)
343: {
344:   PetscFunctionBegin;
345:   if (!*nep) PetscFunctionReturn(PETSC_SUCCESS);
347:   if (--((PetscObject)*nep)->refct > 0) { *nep = NULL; PetscFunctionReturn(PETSC_SUCCESS); }
348:   PetscCall(NEPReset(*nep));
349:   PetscTryTypeMethod(*nep,destroy);
350:   if ((*nep)->eigr) PetscCall(PetscFree4((*nep)->eigr,(*nep)->eigi,(*nep)->errest,(*nep)->perm));
351:   PetscCall(RGDestroy(&(*nep)->rg));
352:   PetscCall(DSDestroy(&(*nep)->ds));
353:   PetscCall(KSPDestroy(&(*nep)->refineksp));
354:   PetscCall(PetscSubcommDestroy(&(*nep)->refinesubc));
355:   PetscCall(PetscFree((*nep)->sc));
356:   /* just in case the initial vectors have not been used */
357:   PetscCall(SlepcBasisDestroy_Private(&(*nep)->nini,&(*nep)->IS));
358:   if ((*nep)->convergeddestroy) PetscCall((*(*nep)->convergeddestroy)(&(*nep)->convergedctx));
359:   if ((*nep)->stoppingdestroy) PetscCall((*(*nep)->stoppingdestroy)(&(*nep)->stoppingctx));
360:   PetscCall(NEPMonitorCancel(*nep));
361:   PetscCall(PetscHeaderDestroy(nep));
362:   PetscFunctionReturn(PETSC_SUCCESS);
363: }

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

368:    Collective

370:    Input Parameters:
371: +  nep - eigensolver context obtained from NEPCreate()
372: -  bv  - the basis vectors object

374:    Note:
375:    Use NEPGetBV() to retrieve the basis vectors context (for example,
376:    to free it at the end of the computations).

378:    Level: advanced

380: .seealso: NEPGetBV()
381: @*/
382: PetscErrorCode NEPSetBV(NEP nep,BV bv)
383: {
384:   PetscFunctionBegin;
387:   PetscCheckSameComm(nep,1,bv,2);
388:   PetscCall(PetscObjectReference((PetscObject)bv));
389:   PetscCall(BVDestroy(&nep->V));
390:   nep->V = bv;
391:   PetscFunctionReturn(PETSC_SUCCESS);
392: }

394: /*@
395:    NEPGetBV - Obtain the basis vectors object associated to the nonlinear
396:    eigensolver object.

398:    Not Collective

400:    Input Parameters:
401: .  nep - eigensolver context obtained from NEPCreate()

403:    Output Parameter:
404: .  bv - basis vectors context

406:    Level: advanced

408: .seealso: NEPSetBV()
409: @*/
410: PetscErrorCode NEPGetBV(NEP nep,BV *bv)
411: {
412:   PetscFunctionBegin;
414:   PetscAssertPointer(bv,2);
415:   if (!nep->V) {
416:     PetscCall(BVCreate(PetscObjectComm((PetscObject)nep),&nep->V));
417:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->V,(PetscObject)nep,0));
418:     PetscCall(PetscObjectSetOptions((PetscObject)nep->V,((PetscObject)nep)->options));
419:   }
420:   *bv = nep->V;
421:   PetscFunctionReturn(PETSC_SUCCESS);
422: }

424: /*@
425:    NEPSetRG - Associates a region object to the nonlinear eigensolver.

427:    Collective

429:    Input Parameters:
430: +  nep - eigensolver context obtained from NEPCreate()
431: -  rg  - the region object

433:    Note:
434:    Use NEPGetRG() to retrieve the region context (for example,
435:    to free it at the end of the computations).

437:    Level: advanced

439: .seealso: NEPGetRG()
440: @*/
441: PetscErrorCode NEPSetRG(NEP nep,RG rg)
442: {
443:   PetscFunctionBegin;
445:   if (rg) {
447:     PetscCheckSameComm(nep,1,rg,2);
448:   }
449:   PetscCall(PetscObjectReference((PetscObject)rg));
450:   PetscCall(RGDestroy(&nep->rg));
451:   nep->rg = rg;
452:   PetscFunctionReturn(PETSC_SUCCESS);
453: }

455: /*@
456:    NEPGetRG - Obtain the region object associated to the
457:    nonlinear eigensolver object.

459:    Not Collective

461:    Input Parameters:
462: .  nep - eigensolver context obtained from NEPCreate()

464:    Output Parameter:
465: .  rg - region context

467:    Level: advanced

469: .seealso: NEPSetRG()
470: @*/
471: PetscErrorCode NEPGetRG(NEP nep,RG *rg)
472: {
473:   PetscFunctionBegin;
475:   PetscAssertPointer(rg,2);
476:   if (!nep->rg) {
477:     PetscCall(RGCreate(PetscObjectComm((PetscObject)nep),&nep->rg));
478:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->rg,(PetscObject)nep,0));
479:     PetscCall(PetscObjectSetOptions((PetscObject)nep->rg,((PetscObject)nep)->options));
480:   }
481:   *rg = nep->rg;
482:   PetscFunctionReturn(PETSC_SUCCESS);
483: }

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

488:    Collective

490:    Input Parameters:
491: +  nep - eigensolver context obtained from NEPCreate()
492: -  ds  - the direct solver object

494:    Note:
495:    Use NEPGetDS() to retrieve the direct solver context (for example,
496:    to free it at the end of the computations).

498:    Level: advanced

500: .seealso: NEPGetDS()
501: @*/
502: PetscErrorCode NEPSetDS(NEP nep,DS ds)
503: {
504:   PetscFunctionBegin;
507:   PetscCheckSameComm(nep,1,ds,2);
508:   PetscCall(PetscObjectReference((PetscObject)ds));
509:   PetscCall(DSDestroy(&nep->ds));
510:   nep->ds = ds;
511:   PetscFunctionReturn(PETSC_SUCCESS);
512: }

514: /*@
515:    NEPGetDS - Obtain the direct solver object associated to the
516:    nonlinear eigensolver object.

518:    Not Collective

520:    Input Parameters:
521: .  nep - eigensolver context obtained from NEPCreate()

523:    Output Parameter:
524: .  ds - direct solver context

526:    Level: advanced

528: .seealso: NEPSetDS()
529: @*/
530: PetscErrorCode NEPGetDS(NEP nep,DS *ds)
531: {
532:   PetscFunctionBegin;
534:   PetscAssertPointer(ds,2);
535:   if (!nep->ds) {
536:     PetscCall(DSCreate(PetscObjectComm((PetscObject)nep),&nep->ds));
537:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->ds,(PetscObject)nep,0));
538:     PetscCall(PetscObjectSetOptions((PetscObject)nep->ds,((PetscObject)nep)->options));
539:   }
540:   *ds = nep->ds;
541:   PetscFunctionReturn(PETSC_SUCCESS);
542: }

544: /*@
545:    NEPRefineGetKSP - Obtain the ksp object used by the eigensolver
546:    object in the refinement phase.

548:    Collective

550:    Input Parameters:
551: .  nep - eigensolver context obtained from NEPCreate()

553:    Output Parameter:
554: .  ksp - ksp context

556:    Level: advanced

558: .seealso: NEPSetRefine()
559: @*/
560: PetscErrorCode NEPRefineGetKSP(NEP nep,KSP *ksp)
561: {
562:   MPI_Comm       comm;

564:   PetscFunctionBegin;
566:   PetscAssertPointer(ksp,2);
567:   if (!nep->refineksp) {
568:     if (nep->npart>1) {
569:       /* Split in subcomunicators */
570:       PetscCall(PetscSubcommCreate(PetscObjectComm((PetscObject)nep),&nep->refinesubc));
571:       PetscCall(PetscSubcommSetNumber(nep->refinesubc,nep->npart));
572:       PetscCall(PetscSubcommSetType(nep->refinesubc,PETSC_SUBCOMM_CONTIGUOUS));
573:       PetscCall(PetscSubcommGetChild(nep->refinesubc,&comm));
574:     } else PetscCall(PetscObjectGetComm((PetscObject)nep,&comm));
575:     PetscCall(KSPCreate(comm,&nep->refineksp));
576:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)nep->refineksp,(PetscObject)nep,0));
577:     PetscCall(PetscObjectSetOptions((PetscObject)nep->refineksp,((PetscObject)nep)->options));
578:     PetscCall(KSPSetOptionsPrefix(*ksp,((PetscObject)nep)->prefix));
579:     PetscCall(KSPAppendOptionsPrefix(*ksp,"nep_refine_"));
580:     PetscCall(KSPSetTolerances(nep->refineksp,SlepcDefaultTol(nep->rtol),PETSC_CURRENT,PETSC_CURRENT,PETSC_CURRENT));
581:   }
582:   *ksp = nep->refineksp;
583:   PetscFunctionReturn(PETSC_SUCCESS);
584: }

586: /*@
587:    NEPSetTarget - Sets the value of the target.

589:    Logically Collective

591:    Input Parameters:
592: +  nep    - eigensolver context
593: -  target - the value of the target

595:    Options Database Key:
596: .  -nep_target <scalar> - the value of the target

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

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

606:    Level: intermediate

608: .seealso: NEPGetTarget(), NEPSetWhichEigenpairs()
609: @*/
610: PetscErrorCode NEPSetTarget(NEP nep,PetscScalar target)
611: {
612:   PetscFunctionBegin;
615:   nep->target = target;
616:   PetscFunctionReturn(PETSC_SUCCESS);
617: }

619: /*@
620:    NEPGetTarget - Gets the value of the target.

622:    Not Collective

624:    Input Parameter:
625: .  nep - eigensolver context

627:    Output Parameter:
628: .  target - the value of the target

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

633:    Level: intermediate

635: .seealso: NEPSetTarget()
636: @*/
637: PetscErrorCode NEPGetTarget(NEP nep,PetscScalar* target)
638: {
639:   PetscFunctionBegin;
641:   PetscAssertPointer(target,2);
642:   *target = nep->target;
643:   PetscFunctionReturn(PETSC_SUCCESS);
644: }

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

650:    Collective

652:    Input Parameters:
653: +  nep - the NEP context
654: .  A   - Function matrix
655: .  B   - preconditioner matrix (usually same as A)
656: .  fun - Function evaluation routine (if NULL then NEP retains any
657:          previously set value), see NEPFunctionFn for the calling sequence
658: -  ctx - [optional] user-defined context for private data for the Function
659:          evaluation routine (may be NULL) (if NULL then NEP retains any
660:          previously set value)

662:    Level: beginner

664: .seealso: NEPGetFunction(), NEPSetJacobian()
665: @*/
666: PetscErrorCode NEPSetFunction(NEP nep,Mat A,Mat B,NEPFunctionFn *fun,void *ctx)
667: {
668:   PetscFunctionBegin;
672:   if (A) PetscCheckSameComm(nep,1,A,2);
673:   if (B) PetscCheckSameComm(nep,1,B,3);

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

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

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

699:    Not Collective

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

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

710:    Level: advanced

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

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

730:    Collective

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

741:    Level: beginner

743: .seealso: NEPSetFunction(), NEPGetJacobian()
744: @*/
745: PetscErrorCode NEPSetJacobian(NEP nep,Mat A,NEPJacobianFn *jac,void *ctx)
746: {
747:   PetscFunctionBegin;
750:   if (A) PetscCheckSameComm(nep,1,A,2);

752:   if (nep->state) PetscCall(NEPReset(nep));
753:   else if (nep->fui && nep->fui!=NEP_USER_INTERFACE_CALLBACK) PetscCall(NEPReset_Problem(nep));

755:   if (jac) nep->computejacobian = jac;
756:   if (ctx) nep->jacobianctx     = ctx;
757:   if (A) {
758:     PetscCall(PetscObjectReference((PetscObject)A));
759:     PetscCall(MatDestroy(&nep->jacobian));
760:     nep->jacobian = A;
761:   }
762:   nep->fui   = NEP_USER_INTERFACE_CALLBACK;
763:   nep->state = NEP_STATE_INITIAL;
764:   PetscFunctionReturn(PETSC_SUCCESS);
765: }

767: /*@C
768:    NEPGetJacobian - Returns the Jacobian matrix and optionally the user
769:    provided routine and context for evaluating the Jacobian.

771:    Not Collective

773:    Input Parameter:
774: .  nep - the nonlinear eigensolver context

776:    Output Parameters:
777: +  A   - location to stash Jacobian matrix (or NULL)
778: .  jac - location to put Jacobian function (or NULL)
779: -  ctx - location to stash Jacobian context (or NULL)

781:    Level: advanced

783: .seealso: NEPSetJacobian()
784: @*/
785: PetscErrorCode NEPGetJacobian(NEP nep,Mat *A,NEPJacobianFn **jac,void **ctx)
786: {
787:   PetscFunctionBegin;
789:   NEPCheckCallback(nep,1);
790:   if (A)   *A   = nep->jacobian;
791:   if (jac) *jac = nep->computejacobian;
792:   if (ctx) *ctx = nep->jacobianctx;
793:   PetscFunctionReturn(PETSC_SUCCESS);
794: }

796: /*@
797:    NEPSetSplitOperator - Sets the operator of the nonlinear eigenvalue problem
798:    in split form.

800:    Collective

802:    Input Parameters:
803: +  nep - the nonlinear eigensolver context
804: .  nt  - number of terms in the split form
805: .  A   - array of matrices
806: .  f   - array of functions
807: -  str - structure flag for matrices

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

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

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

825:    Level: beginner

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

833:   PetscFunctionBegin;
836:   PetscCheck(nt>0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Must have one or more terms, you have %" PetscInt_FMT,nt);
837:   PetscAssertPointer(A,3);
838:   PetscAssertPointer(f,4);

841:   for (i=0;i<nt;i++) {
843:     PetscCheckSameComm(nep,1,A[i],3);
845:     PetscCheckSameComm(nep,1,f[i],4);
846:     PetscCall(MatGetSize(A[i],&m,&n));
847:     PetscCall(MatGetLocalSize(A[i],&mloc,&nloc));
848:     PetscCheck(m==n,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,m,n);
849:     PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"A[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
850:     if (!i) { m0 = m; mloc0 = mloc; }
851:     PetscCheck(m==m0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,m,m0);
852:     PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Local dimensions of A[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
853:     PetscCall(PetscObjectReference((PetscObject)A[i]));
854:     PetscCall(PetscObjectReference((PetscObject)f[i]));
855:   }

857:   if (nep->state && (n!=nep->n || nloc!=nep->nloc)) PetscCall(NEPReset(nep));
858:   else PetscCall(NEPReset_Problem(nep));

860:   /* allocate space and copy matrices and functions */
861:   PetscCall(PetscMalloc1(nt,&nep->A));
862:   for (i=0;i<nt;i++) nep->A[i] = A[i];
863:   PetscCall(PetscMalloc1(nt,&nep->f));
864:   for (i=0;i<nt;i++) nep->f[i] = f[i];
865:   PetscCall(PetscCalloc1(nt,&nep->nrma));
866:   nep->nt    = nt;
867:   nep->mstr  = str;
868:   nep->fui   = NEP_USER_INTERFACE_SPLIT;
869:   nep->state = NEP_STATE_INITIAL;
870:   PetscFunctionReturn(PETSC_SUCCESS);
871: }

873: /*@
874:    NEPGetSplitOperatorTerm - Gets the matrices and functions associated with
875:    the nonlinear operator in split form.

877:    Not Collective

879:    Input Parameters:
880: +  nep - the nonlinear eigensolver context
881: -  k   - the index of the requested term (starting in 0)

883:    Output Parameters:
884: +  A - the matrix of the requested term
885: -  f - the function of the requested term

887:    Level: intermediate

889: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorInfo()
890: @*/
891: PetscErrorCode NEPGetSplitOperatorTerm(NEP nep,PetscInt k,Mat *A,FN *f)
892: {
893:   PetscFunctionBegin;
896:   NEPCheckSplit(nep,1);
897:   PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
898:   if (A) *A = nep->A[k];
899:   if (f) *f = nep->f[k];
900:   PetscFunctionReturn(PETSC_SUCCESS);
901: }

903: /*@
904:    NEPGetSplitOperatorInfo - Returns the number of terms of the split form of
905:    the nonlinear operator, as well as the structure flag for matrices.

907:    Not Collective

909:    Input Parameter:
910: .  nep - the nonlinear eigensolver context

912:    Output Parameters:
913: +  n   - the number of terms passed in NEPSetSplitOperator()
914: -  str - the matrix structure flag passed in NEPSetSplitOperator()

916:    Level: intermediate

918: .seealso: NEPSetSplitOperator(), NEPGetSplitOperatorTerm()
919: @*/
920: PetscErrorCode NEPGetSplitOperatorInfo(NEP nep,PetscInt *n,MatStructure *str)
921: {
922:   PetscFunctionBegin;
924:   NEPCheckSplit(nep,1);
925:   if (n)   *n = nep->nt;
926:   if (str) *str = nep->mstr;
927:   PetscFunctionReturn(PETSC_SUCCESS);
928: }

930: /*@
931:    NEPSetSplitPreconditioner - Sets an operator in split form from which
932:    to build the preconditioner to be used when solving the nonlinear
933:    eigenvalue problem in split form.

935:    Collective

937:    Input Parameters:
938: +  nep  - the nonlinear eigensolver context
939: .  ntp  - number of terms in the split preconditioner
940: .  P    - array of matrices
941: -  strp - structure flag for matrices

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

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

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

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

960:    Level: advanced

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

968:   PetscFunctionBegin;
971:   PetscCheck(ntp>=0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"Negative value of ntp = %" PetscInt_FMT,ntp);
972:   PetscCheck(nep->fui==NEP_USER_INTERFACE_SPLIT,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"Must call NEPSetSplitOperator first");
973:   PetscCheck(ntp==0 || nep->nt==ntp,PetscObjectComm((PetscObject)nep),PETSC_ERR_SUP,"The number of terms must be the same as in NEPSetSplitOperator()");
974:   if (ntp) PetscAssertPointer(P,3);

977:   for (i=0;i<ntp;i++) {
979:     PetscCheckSameComm(nep,1,P[i],3);
980:     PetscCall(MatGetSize(P[i],&m,&n));
981:     PetscCall(MatGetLocalSize(P[i],&mloc,&nloc));
982:     PetscCheck(m==n,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"P[%" PetscInt_FMT "] is a non-square matrix (%" PetscInt_FMT " rows, %" PetscInt_FMT " cols)",i,m,n);
983:     PetscCheck(mloc==nloc,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_WRONG,"P[%" PetscInt_FMT "] does not have equal row and column local sizes (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,nloc);
984:     if (!i) { m0 = m; mloc0 = mloc; }
985:     PetscCheck(m==m0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Dimensions of P[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,m,m0);
986:     PetscCheck(mloc==mloc0,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_INCOMP,"Local dimensions of P[%" PetscInt_FMT "] do not match with previous matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")",i,mloc,mloc0);
987:     PetscCall(PetscObjectReference((PetscObject)P[i]));
988:   }

990:   PetscCheck(!nep->state,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"To call this function after NEPSetUp(), you must call NEPSetSplitOperator() again");
991:   if (nep->P) PetscCall(MatDestroyMatrices(nep->nt,&nep->P));

993:   /* allocate space and copy matrices */
994:   if (ntp) {
995:     PetscCall(PetscMalloc1(ntp,&nep->P));
996:     for (i=0;i<ntp;i++) nep->P[i] = P[i];
997:   }
998:   nep->mstrp = strp;
999:   nep->state = NEP_STATE_INITIAL;
1000:   PetscFunctionReturn(PETSC_SUCCESS);
1001: }

1003: /*@
1004:    NEPGetSplitPreconditionerTerm - Gets the matrices associated with
1005:    the split preconditioner.

1007:    Not Collective

1009:    Input Parameters:
1010: +  nep - the nonlinear eigensolver context
1011: -  k   - the index of the requested term (starting in 0)

1013:    Output Parameter:
1014: .  P  - the matrix of the requested term

1016:    Level: advanced

1018: .seealso: NEPSetSplitPreconditioner(), NEPGetSplitPreconditionerInfo()
1019: @*/
1020: PetscErrorCode NEPGetSplitPreconditionerTerm(NEP nep,PetscInt k,Mat *P)
1021: {
1022:   PetscFunctionBegin;
1025:   PetscAssertPointer(P,3);
1026:   NEPCheckSplit(nep,1);
1027:   PetscCheck(k>=0 && k<nep->nt,PetscObjectComm((PetscObject)nep),PETSC_ERR_ARG_OUTOFRANGE,"k must be between 0 and %" PetscInt_FMT,nep->nt-1);
1028:   PetscCheck(nep->P,PetscObjectComm((PetscObject)nep),PETSC_ERR_ORDER,"You have not called NEPSetSplitPreconditioner()");
1029:   *P = nep->P[k];
1030:   PetscFunctionReturn(PETSC_SUCCESS);
1031: }

1033: /*@
1034:    NEPGetSplitPreconditionerInfo - Returns the number of terms of the split
1035:    preconditioner, as well as the structure flag for matrices.

1037:    Not Collective

1039:    Input Parameter:
1040: .  nep - the nonlinear eigensolver context

1042:    Output Parameters:
1043: +  n    - the number of terms passed in NEPSetSplitPreconditioner()
1044: -  strp - the matrix structure flag passed in NEPSetSplitPreconditioner()

1046:    Level: advanced

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